Reference documentation for deal.II version 8.1.0
dof_info.templates.h
1 // ---------------------------------------------------------------------
2 // @f$Id: dof_info.templates.h 31932 2013-12-08 02:15:54Z heister @f$
3 //
4 // Copyright (C) 2011 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 
18 #include <deal.II/base/memory_consumption.h>
19 #include <deal.II/base/multithread_info.h>
20 #include <deal.II/lac/compressed_simple_sparsity_pattern.h>
21 #include <deal.II/lac/sparsity_pattern.h>
22 #include <deal.II/matrix_free/dof_info.h>
23 #include <deal.II/matrix_free/helper_functions.h>
24 
26 
27 namespace internal
28 {
29  namespace MatrixFreeFunctions
30  {
31 
33  {
34  bool operator()(const std::pair<types::global_dof_index,double> &p1,
35  const std::pair<types::global_dof_index,double> &p2) const
36  {
37  return p1.second < p2.second;
38  }
39  };
40 
45  template <typename Number>
46  struct ConstraintValues
47  {
49 
57  unsigned short
58  insert_entries (const std::vector<std::pair<types::global_dof_index,double> > &entries);
59 
60  std::vector<std::pair<types::global_dof_index, double> > constraint_entries;
61  std::vector<types::global_dof_index> constraint_indices;
62 
63  std::pair<std::vector<Number>, types::global_dof_index> next_constraint;
64  std::map<std::vector<Number>, types::global_dof_index, FPArrayComparator<double> > constraints;
65  };
66 
67 
68  template <typename Number>
70  :
71  constraints(FPArrayComparator<double>(1.))
72  {}
73 
74  template <typename Number>
75  unsigned short
77  insert_entries (const std::vector<std::pair<types::global_dof_index,double> > &entries)
78  {
79  next_constraint.first.resize(entries.size());
80  if (entries.size() > 0)
81  {
82  constraint_indices.resize(entries.size());
83  constraint_entries = entries;
84  std::sort(constraint_entries.begin(), constraint_entries.end(),
86  for (types::global_dof_index j=0; j<constraint_entries.size(); j++)
87  {
88  // copy the indices of the constraint entries after sorting.
89  constraint_indices[j] = constraint_entries[j].first;
90 
91  // one_constraint takes the weights of the constraint
92  next_constraint.first[j] = constraint_entries[j].second;
93  }
94  }
95  next_constraint.second = constraints.size();
96 
97  // check whether or not constraint is already in pool. the initial
98  // implementation computed a hash value based on the truncated array (to
99  // given accuracy around 1e-13) in order to easily detect different
100  // arrays and then made a fine-grained check when the hash values were
101  // equal. this was quite lenghty and now we use a std::map with a
102  // user-defined comparator to compare floating point arrays to a
103  // tolerance 1e-13.
104  std::pair<typename std::map<std::vector<double>, types::global_dof_index,
105  FPArrayComparator<double> >::iterator,
106  bool> it = constraints.insert(next_constraint);
107 
108  types::global_dof_index insert_position = deal_II_numbers::invalid_dof_index;
109  if (it.second == false)
110  insert_position = it.first->second;
111  else
112  insert_position = next_constraint.second;
113 
114  // we want to store the result as a short variable, so we have to make
115  // sure that the result does not exceed the limits when casting.
116  Assert(insert_position < (1<<(8*sizeof(unsigned short))),
117  ExcInternalError());
118  return static_cast<unsigned short>(insert_position);
119  }
120 
121 
122 
123  // ----------------- actual DoFInfo functions -----------------------------
124 
125  DoFInfo::DoFInfo ()
126  {
127  clear();
128  }
129 
130 
131  DoFInfo::DoFInfo (const DoFInfo &dof_info_in)
132  :
133  row_starts (dof_info_in.row_starts),
134  dof_indices (dof_info_in.dof_indices),
135  constraint_indicator (dof_info_in.constraint_indicator),
136  vector_partitioner (dof_info_in.vector_partitioner),
137  constrained_dofs (dof_info_in.constrained_dofs),
138  row_starts_plain_indices (dof_info_in.row_starts_plain_indices),
139  plain_dof_indices (dof_info_in.plain_dof_indices),
140  dimension (dof_info_in.dimension),
141  n_components (dof_info_in.n_components),
142  dofs_per_cell (dof_info_in.dofs_per_cell),
143  dofs_per_face (dof_info_in.dofs_per_face),
144  store_plain_indices (dof_info_in.store_plain_indices),
145  cell_active_fe_index (dof_info_in.cell_active_fe_index),
146  max_fe_index (dof_info_in.max_fe_index),
147  fe_index_conversion (dof_info_in.fe_index_conversion),
148  ghost_dofs (dof_info_in.ghost_dofs)
149  {}
150 
151 
152 
153  void
154  DoFInfo::clear ()
155  {
156  row_starts.clear();
157  dof_indices.clear();
158  constraint_indicator.clear();
159  vector_partitioner.reset();
160  ghost_dofs.clear();
161  dofs_per_cell.clear();
162  dofs_per_face.clear();
163  dimension = 2;
164  n_components = 0;
165  row_starts_plain_indices.clear();
166  plain_dof_indices.clear();
167  store_plain_indices = false;
168  cell_active_fe_index.clear();
169  max_fe_index = 0;
170  fe_index_conversion.clear();
171  }
172 
173 
174 
175  void
176  DoFInfo::read_dof_indices (const std::vector<types::global_dof_index> &local_indices,
177  const std::vector<unsigned int> &lexicographic_inv,
178  const ConstraintMatrix &constraints,
179  const unsigned int cell_number,
180  ConstraintValues<double> &constraint_values,
181  bool &cell_at_boundary)
182  {
183  Assert (vector_partitioner.get() !=0, ExcInternalError());
184  const unsigned int n_mpi_procs = vector_partitioner->n_mpi_processes();
185  const types::global_dof_index first_owned = vector_partitioner->local_range().first;
186  const types::global_dof_index last_owned = vector_partitioner->local_range().second;
187  Assert (last_owned-first_owned < std::numeric_limits<unsigned int>::max(),
188  ExcMessage("The size local range of owned indices must not "
189  "exceed the size of unsigned int"));
190  const unsigned int n_owned = last_owned - first_owned;
191  std::pair<unsigned short,unsigned short> constraint_iterator (0,0);
192 
193  unsigned int dofs_this_cell = (cell_active_fe_index.empty()) ?
194  dofs_per_cell[0] : dofs_per_cell[cell_active_fe_index[cell_number]];
195  for (unsigned int i=0; i<dofs_this_cell; i++)
196  {
197  types::global_dof_index current_dof =
198  local_indices[lexicographic_inv[i]];
199  const std::vector<std::pair<types::global_dof_index,double> >
200  *entries_ptr =
201  constraints.get_constraint_entries(current_dof);
202 
203  // dof is constrained
204  if (entries_ptr != 0)
205  {
206  // in case we want to access plain indices, we need to know
207  // about the location of constrained indices as well (all the
208  // other indices are collected by the cases below)
209  if (current_dof < first_owned || current_dof >= last_owned)
210  {
211  ghost_dofs.push_back (current_dof);
212  cell_at_boundary = true;
213  }
214 
215  // check whether this dof is identity constrained to another
216  // dof. then we can simply insert that dof and there is no need
217  // to actually resolve the constraint entries
218  const std::vector<std::pair<types::global_dof_index,double> >
219  &entries = *entries_ptr;
220  const types::global_dof_index n_entries = entries.size();
221  if (n_entries == 1 && std::fabs(entries[0].second-1.)<1e-14)
222  {
223  current_dof = entries[0].first;
224  goto no_constraint;
225  }
226 
227  // append a new index to the indicators
228  constraint_indicator.push_back (constraint_iterator);
229  constraint_indicator.back().second =
230  constraint_values.insert_entries (entries);
231 
232  // reset constraint iterator for next round
233  constraint_iterator.first = 0;
234 
235  // add the local_to_global indices computed in the
236  // insert_entries function. transform the index to local index
237  // space or mark it as ghost if necessary
238  if (n_entries > 0)
239  {
240  const std::vector<types::global_dof_index> &constraint_indices =
241  constraint_values.constraint_indices;
242  for (unsigned int j=0; j<n_entries; ++j)
243  {
244  if (n_mpi_procs > 1 &&
245  (constraint_indices[j] < first_owned ||
246  constraint_indices[j] >= last_owned))
247  {
248  dof_indices.push_back (n_owned + ghost_dofs.size());
249 
250  // collect ghosts so that we can later construct an
251  // IndexSet for them. also store whether the current
252  // cell is on the boundary
253  ghost_dofs.push_back(constraint_indices[j]);
254  cell_at_boundary = true;
255  }
256  else
257  // not ghost, so transform to the local index space
258  // directly
259  dof_indices.push_back
260  (static_cast<unsigned int>(constraint_indices[j] -
261  first_owned));
262  }
263  }
264  }
265  else
266  {
267 no_constraint:
268  // Not constrained, we simply have to add the local index to the
269  // indices_local_to_global list and increment constraint
270  // iterator. transform to local index space/mark as ghost
271  if (n_mpi_procs > 1 &&
272  (current_dof < first_owned ||
273  current_dof >= last_owned))
274  {
275  ghost_dofs.push_back(current_dof);
276  current_dof = n_owned + ghost_dofs.size()-1;
277  cell_at_boundary = true;
278  }
279  else
280  current_dof -= first_owned;
281 
282  dof_indices.push_back (static_cast<unsigned int>(current_dof));
283 
284  // make sure constraint_iterator.first is always within the
285  // bounds of unsigned short
286  Assert (constraint_iterator.first <
287  (1<<(8*sizeof(unsigned short)))-1,
288  ExcInternalError());
289  constraint_iterator.first++;
290  }
291  }
292  row_starts[cell_number+1][0] = dof_indices.size();
293  row_starts[cell_number+1][1] = constraint_indicator.size();
294  row_starts[cell_number+1][2] = 0;
295 
296  // now to the plain indices: in case we have constraints on this cell,
297  // store the indices without the constraints resolve once again
298  if (store_plain_indices == true)
299  {
300  if (cell_number == 0)
301  row_starts_plain_indices.resize (row_starts.size());
302  row_starts_plain_indices[cell_number] = plain_dof_indices.size();
303  bool cell_has_constraints = (row_starts[cell_number+1][1] >
304  row_starts[cell_number][1]);
305  if (cell_has_constraints == true)
306  {
307  for (unsigned int i=0; i<dofs_this_cell; ++i)
308  {
309  types::global_dof_index current_dof =
310  local_indices[lexicographic_inv[i]];
311  if (n_mpi_procs > 1 &&
312  (current_dof < first_owned ||
313  current_dof >= last_owned))
314  {
315  ghost_dofs.push_back(current_dof);
316  current_dof = n_owned + ghost_dofs.size()-1;
317  cell_at_boundary = true;
318  }
319  else
320  current_dof -= first_owned;
321  plain_dof_indices.push_back (static_cast<unsigned int>
322  (current_dof));
323  }
324  }
325  }
326  }
327 
328 
329 
330  void
331  DoFInfo::assign_ghosts (const std::vector<unsigned int> &boundary_cells)
332  {
333  Assert (boundary_cells.size() < row_starts.size(), ExcInternalError());
334 
335  // sort ghost dofs and compress out duplicates
336  const unsigned int n_owned = (vector_partitioner->local_range().second-
337  vector_partitioner->local_range().first);
338  const std::size_t n_ghosts = ghost_dofs.size();
339  unsigned int n_unique_ghosts= 0;
340 #ifdef DEBUG
341  for (std::vector<unsigned int>::iterator dof = dof_indices.begin();
342  dof!=dof_indices.end(); ++dof)
343  AssertIndexRange (*dof, n_owned+n_ghosts);
344 #endif
345 
346  std::vector<unsigned int> ghost_numbering (n_ghosts);
347  IndexSet ghost_indices (vector_partitioner->size());
348  if (n_ghosts > 0)
349  {
350  // since we need to go back to the local_to_global indices and
351  // replace the temporary numbering of ghosts by the real number in
352  // the index set, we need to store these values
353  std::vector<std::pair<types::global_dof_index,unsigned int> > ghost_origin(n_ghosts);
354  for (std::size_t i=0; i<n_ghosts; ++i)
355  {
356  ghost_origin[i].first = ghost_dofs[i];
357  ghost_origin[i].second = i;
358  }
359  std::sort (ghost_origin.begin(), ghost_origin.end());
360 
361  types::global_dof_index last_contiguous_start = ghost_origin[0].first;
362  ghost_numbering[ghost_origin[0].second] = 0;
363  for (std::size_t i=1; i<n_ghosts; i++)
364  {
365  if (ghost_origin[i].first > ghost_origin[i-1].first+1)
366  {
367  ghost_indices.add_range (last_contiguous_start,
368  ghost_origin[i-1].first+1);
369  last_contiguous_start = ghost_origin[i].first;
370  }
371  if (ghost_origin[i].first>ghost_origin[i-1].first)
372  ++n_unique_ghosts;
373  ghost_numbering[ghost_origin[i].second] = n_unique_ghosts;
374  }
375  ++n_unique_ghosts;
376  ghost_indices.add_range (last_contiguous_start,
377  ghost_origin.back().first+1);
378  ghost_indices.compress();
379 
380  // make sure that we got the correct local numbering of the ghost
381  // dofs. the ghost index set should store the same number
382  {
383  AssertDimension (n_unique_ghosts, ghost_indices.n_elements());
384  for (std::size_t i=0; i<n_ghosts; ++i)
385  Assert (ghost_numbering[i] ==
386  ghost_indices.index_within_set(ghost_dofs[i]),
387  ExcInternalError());
388  }
389 
390  // apply correct numbering for ghost indices: We previously just
391  // enumerated them according to their appearance in the
392  // local_to_global structure. Above, we derived a relation between
393  // this enumeration and the actual number
394  const unsigned int n_boundary_cells = boundary_cells.size();
395  for (unsigned int i=0; i<n_boundary_cells; ++i)
396  {
397  unsigned int *data_ptr = const_cast<unsigned int *> (begin_indices(boundary_cells[i]));
398 
399  const unsigned int *row_end = end_indices(boundary_cells[i]);
400  for ( ; data_ptr != row_end; ++data_ptr)
401  *data_ptr = ((*data_ptr < n_owned)
402  ?
403  *data_ptr
404  :
405  n_owned +
406  ghost_numbering[*data_ptr - n_owned]);
407 
408  // now the same procedure for plain indices
409  if (store_plain_indices == true)
410  {
411  if (row_length_indicators(boundary_cells[i]) > 0)
412  {
413  unsigned int *data_ptr = const_cast<unsigned int *> (begin_indices_plain(boundary_cells[i]));
414  const unsigned int *row_end = end_indices_plain(boundary_cells[i]);
415  for ( ; data_ptr != row_end; ++data_ptr)
416  *data_ptr = ((*data_ptr < n_owned)
417  ?
418  *data_ptr
419  :
420  n_owned +
421  ghost_numbering[*data_ptr - n_owned]);
422  }
423  }
424  }
425  }
426 
427  std::vector<types::global_dof_index> empty;
428  ghost_dofs.swap(empty);
429 
430  // set the ghost indices now. need to cast away constness here, but that
431  // is uncritical since we reset the Partitioner in the same initialize
432  // call as this call here.
433  Utilities::MPI::Partitioner *vec_part =
434  const_cast<Utilities::MPI::Partitioner *>(vector_partitioner.get());
435  vec_part->set_ghost_indices (ghost_indices);
436  }
437 
438 
439 
440  void
441  DoFInfo::compute_renumber_serial (const std::vector<unsigned int> &boundary_cells,
442  const SizeInfo &size_info,
443  std::vector<unsigned int> &renumbering)
444  {
445  std::vector<unsigned int> reverse_numbering (size_info.n_active_cells,
447  const unsigned int n_boundary_cells = boundary_cells.size();
448  for (unsigned int j=0; j<n_boundary_cells; ++j)
449  reverse_numbering[boundary_cells[j]] =
450  j + size_info.vectorization_length*size_info.boundary_cells_start;
451  unsigned int counter = 0;
452  unsigned int j = 0;
453  while (counter < size_info.n_active_cells &&
454  counter < size_info.vectorization_length * size_info.boundary_cells_start)
455  {
456  if (reverse_numbering[j] == numbers::invalid_unsigned_int)
457  reverse_numbering[j] = counter++;
458  j++;
459  }
460  counter = std::min (size_info.vectorization_length*
461  size_info.boundary_cells_start+n_boundary_cells,
462  size_info.n_active_cells);
463  if (counter < size_info.n_active_cells)
464  {
465  for ( ; j<size_info.n_active_cells; ++j)
466  if (reverse_numbering[j] == numbers::invalid_unsigned_int)
467  reverse_numbering[j] = counter++;
468  }
469  AssertDimension (counter, size_info.n_active_cells);
470  renumbering = Utilities::invert_permutation (reverse_numbering);
471  }
472 
473 
474 
475  void
476  DoFInfo::compute_renumber_hp_serial (SizeInfo &size_info,
477  std::vector<unsigned int> &renumbering,
478  std::vector<unsigned int> &irregular_cells)
479  {
480  if (max_fe_index < 2)
481  return;
482  const unsigned int n_active_cells = size_info.n_active_cells;
483  const unsigned int vectorization_length = size_info.vectorization_length;
484  irregular_cells.resize (0);
485  irregular_cells.resize (size_info.n_macro_cells+3*max_fe_index);
486  std::vector<std::vector<unsigned int> > renumbering_fe_index;
487  renumbering_fe_index.resize(max_fe_index);
488  unsigned int counter,n_macro_cells_before = 0;
489  const unsigned int
490  start_bound = std::min (size_info.n_active_cells,
491  size_info.boundary_cells_start*vectorization_length),
492  end_bound = std::min (size_info.n_active_cells,
493  size_info.boundary_cells_end*vectorization_length);
494  for (counter=0; counter<start_bound; counter++)
495  {
496  renumbering_fe_index[cell_active_fe_index[renumbering[counter]]].
497  push_back(renumbering[counter]);
498  }
499  counter = 0;
500  for (unsigned int j=0; j<max_fe_index; j++)
501  {
502  for (unsigned int jj=0; jj<renumbering_fe_index[j].size(); jj++)
503  renumbering[counter++] = renumbering_fe_index[j][jj];
504  irregular_cells[renumbering_fe_index[j].size()/vectorization_length+
505  n_macro_cells_before] =
506  renumbering_fe_index[j].size()%vectorization_length;
507  n_macro_cells_before += (renumbering_fe_index[j].size()+vectorization_length-1)/
508  vectorization_length;
509  renumbering_fe_index[j].resize(0);
510  }
511  unsigned int new_boundary_start = n_macro_cells_before;
512  for (counter = start_bound; counter < end_bound; counter++)
513  {
514  renumbering_fe_index[cell_active_fe_index[renumbering[counter]]].
515  push_back(renumbering[counter]);
516  }
517  counter = start_bound;
518  for (unsigned int j=0; j<max_fe_index; j++)
519  {
520  for (unsigned int jj=0; jj<renumbering_fe_index[j].size(); jj++)
521  renumbering[counter++] = renumbering_fe_index[j][jj];
522  irregular_cells[renumbering_fe_index[j].size()/vectorization_length+
523  n_macro_cells_before] =
524  renumbering_fe_index[j].size()%vectorization_length;
525  n_macro_cells_before += (renumbering_fe_index[j].size()+vectorization_length-1)/
526  vectorization_length;
527  renumbering_fe_index[j].resize(0);
528  }
529  unsigned int new_boundary_end = n_macro_cells_before;
530  for (counter=end_bound; counter<n_active_cells; counter++)
531  {
532  renumbering_fe_index[cell_active_fe_index[renumbering[counter]]].
533  push_back(renumbering[counter]);
534  }
535  counter = end_bound;
536  for (unsigned int j=0; j<max_fe_index; j++)
537  {
538  for (unsigned int jj=0; jj<renumbering_fe_index[j].size(); jj++)
539  renumbering[counter++] = renumbering_fe_index[j][jj];
540  irregular_cells[renumbering_fe_index[j].size()/vectorization_length+
541  n_macro_cells_before] =
542  renumbering_fe_index[j].size()%vectorization_length;
543  n_macro_cells_before += (renumbering_fe_index[j].size()+vectorization_length-1)/
544  vectorization_length;
545  }
546  AssertIndexRange (n_macro_cells_before,
547  size_info.n_macro_cells + 3*max_fe_index+1);
548  irregular_cells.resize (n_macro_cells_before);
549  size_info.n_macro_cells = n_macro_cells_before;
550  size_info.boundary_cells_start = new_boundary_start;
551  size_info.boundary_cells_end = new_boundary_end;
552  }
553 
554 
555 
556  void
557  DoFInfo::compute_renumber_parallel (const std::vector<unsigned int> &boundary_cells,
558  SizeInfo &size_info,
559  std::vector<unsigned int> &renumbering)
560  {
561  std::vector<unsigned int> reverse_numbering (size_info.n_active_cells,
563  const unsigned int n_boundary_cells = boundary_cells.size();
564  for (unsigned int j=0; j<n_boundary_cells; ++j)
565  reverse_numbering[boundary_cells[j]] = j;
566  unsigned int counter = n_boundary_cells;
567  for (unsigned int j=0; j<size_info.n_active_cells; ++j)
568  if (reverse_numbering[j] == numbers::invalid_unsigned_int)
569  reverse_numbering[j] = counter++;
570 
571  size_info.boundary_cells_end = (size_info.boundary_cells_end -
572  size_info.boundary_cells_start);
573  size_info.boundary_cells_start = 0;
574 
575  AssertDimension (counter, size_info.n_active_cells);
576  renumbering = Utilities::invert_permutation (reverse_numbering);
577  }
578 
579 
580 
581  void
582  DoFInfo::reorder_cells (const SizeInfo &size_info,
583  const std::vector<unsigned int> &renumbering,
584  const std::vector<unsigned int> &constraint_pool_row_index,
585  const std::vector<unsigned int> &irregular_cells,
586  const unsigned int vectorization_length)
587  {
588  // first reorder the active fe index.
589  if (cell_active_fe_index.size() > 0)
590  {
591  std::vector<unsigned int> new_active_fe_index;
592  new_active_fe_index.reserve (size_info.n_macro_cells);
593  std::vector<unsigned int> fe_indices(vectorization_length);
594  unsigned int position_cell = 0;
595  for (unsigned int cell=0; cell<size_info.n_macro_cells; ++cell)
596  {
597  const unsigned int n_comp = (irregular_cells[cell] > 0 ?
598  irregular_cells[cell] : vectorization_length);
599  for (unsigned int j=0; j<n_comp; ++j)
600  fe_indices[j]=cell_active_fe_index[renumbering[position_cell+j]];
601 
602  // by construction, all cells should have the same fe index.
603  for (unsigned int j=1; j<n_comp; ++j)
604  Assert (fe_indices[j] == fe_indices[0], ExcInternalError());
605 
606  new_active_fe_index.push_back(fe_indices[0]);
607  position_cell += n_comp;
608  }
609  std::swap (new_active_fe_index, cell_active_fe_index);
610  }
611 
612  std::vector<std_cxx1x::array<unsigned int, 3> > new_row_starts;
613  std::vector<unsigned int> new_dof_indices;
614  std::vector<std::pair<unsigned short,unsigned short> >
615  new_constraint_indicator;
616  std::vector<unsigned int> new_plain_indices, new_rowstart_plain;
617  unsigned int position_cell = 0;
618  new_row_starts.resize (size_info.n_macro_cells + 1);
619  new_dof_indices.reserve (dof_indices.size());
620  new_constraint_indicator.reserve (constraint_indicator.size());
621  if (store_plain_indices == true)
622  {
623  new_rowstart_plain.resize (size_info.n_macro_cells + 1,
625  new_plain_indices.reserve (plain_dof_indices.size());
626  }
627 
628  // copy the indices and the constraint indicators to the new data field:
629  // Store the indices in a way so that adjacent data fields in local
630  // vectors are adjacent, i.e., first dof index 0 for all vectors, then
631  // dof index 1 for all vectors, and so on. This involves some extra
632  // resorting.
633  std::vector<const unsigned int *> glob_indices (vectorization_length);
634  std::vector<const unsigned int *> plain_glob_indices (vectorization_length);
635  std::vector<const std::pair<unsigned short,unsigned short>*>
636  constr_ind(vectorization_length), constr_end(vectorization_length);
637  std::vector<unsigned int> index(vectorization_length);
638  for (unsigned int i=0; i<size_info.n_macro_cells; ++i)
639  {
640  const unsigned int dofs_mcell =
641  dofs_per_cell[cell_active_fe_index.size() == 0 ? 0 :
642  cell_active_fe_index[i]] * vectorization_length;
643  new_row_starts[i][0] = new_dof_indices.size();
644  new_row_starts[i][1] = new_constraint_indicator.size();
645  new_row_starts[i][2] = irregular_cells[i];
646 
647  const unsigned int n_comp = (irregular_cells[i]>0 ?
648  irregular_cells[i] : vectorization_length);
649 
650  for (unsigned int j=0; j<n_comp; ++j)
651  {
652  glob_indices[j] = begin_indices(renumbering[position_cell+j]);
653  constr_ind[j] = begin_indicators(renumbering[position_cell+j]);
654  constr_end[j] = end_indicators(renumbering[position_cell+j]);
655  index[j] = 0;
656  }
657 
658  bool has_constraints = false;
659  if (store_plain_indices == true)
660  {
661  for (unsigned int j=0; j<n_comp; ++j)
662  if (begin_indicators(renumbering[position_cell+j]) <
663  end_indicators(renumbering[position_cell+j]))
664  {
665  plain_glob_indices[j] =
666  begin_indices_plain (renumbering[position_cell+j]);
667  has_constraints = true;
668  }
669  else
670  plain_glob_indices[j] =
671  begin_indices (renumbering[position_cell+j]);
672  if (has_constraints == true)
673  new_rowstart_plain[i] = new_plain_indices.size();
674  }
675 
676  unsigned int m_ind_local = 0, m_index = 0;
677  while (m_ind_local < dofs_mcell)
678  for (unsigned int j=0; j<vectorization_length; ++j)
679  {
680  // last cell: nothing to do
681  if (j >= n_comp)
682  {
683  ++m_ind_local;
684  continue;
685  }
686 
687  // otherwise, check if we are a constrained dof. The dof is
688  // not constrained if we are at the end of the row for the
689  // constraints (indi[j] == n_indi[j]) or if the local index[j]
690  // is smaller than the next position for a constraint. Then,
691  // just copy it. otherwise, copy all the entries that come
692  // with this dof
693  if (constr_ind[j] == constr_end[j] ||
694  index[j] < constr_ind[j]->first)
695  {
696  new_dof_indices.push_back (*glob_indices[j]);
697  ++m_index;
698  ++index[j];
699  ++glob_indices[j];
700  }
701  else
702  {
703  const unsigned short constraint_loc = constr_ind[j]->second;
704  new_constraint_indicator.push_back
705  (std::pair<unsigned short,unsigned short> (m_index, constraint_loc));
706  for (unsigned int k=constraint_pool_row_index[constraint_loc];
707  k<constraint_pool_row_index[constraint_loc+1];
708  ++k, ++glob_indices[j])
709  new_dof_indices.push_back (*glob_indices[j]);
710  ++constr_ind[j];
711  m_index = 0;
712  index[j] = 0;
713  }
714  if (store_plain_indices==true && has_constraints==true)
715  new_plain_indices.push_back (*plain_glob_indices[j]++);
716  ++m_ind_local;
717  }
718 
719  for (unsigned int j=0; j<n_comp; ++j)
720  Assert (glob_indices[j]==end_indices(renumbering[position_cell+j]),
721  ExcInternalError());
722  position_cell += n_comp;
723  }
724  AssertDimension (position_cell+1, row_starts.size());
725 
726  new_row_starts[size_info.n_macro_cells][0] = new_dof_indices.size();
727  new_row_starts[size_info.n_macro_cells][1] = new_constraint_indicator.size();
728  new_row_starts[size_info.n_macro_cells][2] = 0;
729 
730  AssertDimension(dof_indices.size(), new_dof_indices.size());
731  AssertDimension(constraint_indicator.size(),
732  new_constraint_indicator.size());
733 
734  new_row_starts.swap (row_starts);
735  new_dof_indices.swap (dof_indices);
736  new_constraint_indicator.swap (constraint_indicator);
737  new_plain_indices.swap (plain_dof_indices);
738  new_rowstart_plain.swap (row_starts_plain_indices);
739 
740 #ifdef DEBUG
741  // sanity check 1: all indices should be smaller than the number of dofs
742  // locally owned plus the number of ghosts
743  const unsigned int index_range = (vector_partitioner->local_range().second-
744  vector_partitioner->local_range().first)
745  + vector_partitioner->ghost_indices().n_elements();
746  for (std::size_t i=0; i<dof_indices.size(); ++i)
747  AssertIndexRange (dof_indices[i], index_range);
748 
749  // sanity check 2: for the constraint indicators, the first index should
750  // be smaller than the number of indices in the row, and the second
751  // index should be smaller than the number of constraints in the
752  // constraint pool.
753  for (unsigned int row=0; row<size_info.n_macro_cells; ++row)
754  {
755  const unsigned int row_length_ind = row_length_indices(row);
756  const std::pair<unsigned short,unsigned short>
757  *con_it = begin_indicators(row), * end_con = end_indicators(row);
758  for ( ; con_it != end_con; ++con_it)
759  {
760  AssertIndexRange (con_it->first, row_length_ind+1);
761  AssertIndexRange (con_it->second,
762  constraint_pool_row_index.size()-1);
763  }
764  }
765 
766  // sanity check 3: all non-boundary cells should have indices that only
767  // refer to the locally owned range
768  const unsigned int local_size = (vector_partitioner->local_range().second-
769  vector_partitioner->local_range().first);
770  for (unsigned int row=0; row<size_info.boundary_cells_start; ++row)
771  {
772  const unsigned int *ptr = begin_indices(row);
773  const unsigned int *end_ptr = end_indices (row);
774  for ( ; ptr != end_ptr; ++ptr)
775  AssertIndexRange (*ptr, local_size);
776  }
777  for (unsigned int row=size_info.boundary_cells_end;
778  row<size_info.n_macro_cells; ++row)
779  {
780  const unsigned int *ptr = begin_indices(row);
781  const unsigned int *end_ptr = end_indices (row);
782  for ( ; ptr != end_ptr; ++ptr)
783  AssertIndexRange (*ptr, local_size);
784  }
785 #endif
786  }
787 
788 
789 
790  void DoFInfo::guess_block_size (const SizeInfo &size_info,
791  TaskInfo &task_info)
792  {
793  // user did not say a positive number, so we have to guess
794  if (task_info.block_size == 0)
795  {
796  // we would like to have enough work to do, so as first guess, try
797  // to get 50 times as many chunks as we have threads on the system.
798  task_info.block_size =
799  size_info.n_macro_cells / (multithread_info.n_threads() * 50);
800 
801  // if there are too few degrees of freedom per cell, need to
802  // increase the block size
803  const unsigned int minimum_parallel_grain_size = 500;
804  if (dofs_per_cell[0] * task_info.block_size <
805  minimum_parallel_grain_size)
806  task_info.block_size = (minimum_parallel_grain_size /
807  dofs_per_cell[0] + 1);
808  }
809  if (task_info.block_size > size_info.n_macro_cells)
810  task_info.block_size = size_info.n_macro_cells;
811  }
812 
813 
814 
815  void DoFInfo::make_thread_graph_partition_color
816  (SizeInfo &size_info,
817  TaskInfo &task_info,
818  std::vector<unsigned int> &renumbering,
819  std::vector<unsigned int> &irregular_cells,
820  const bool hp_bool)
821  {
822  if (size_info.n_macro_cells == 0)
823  return;
824 
825  const std::size_t vectorization_length = size_info.vectorization_length;
826  Assert (vectorization_length > 0, ExcInternalError());
827 
828  guess_block_size (size_info, task_info);
829 
830  // set up partitions. if we just use coloring without partitions, do
831  // nothing here, assume all cells to belong to the zero partition (that
832  // we otherwise use for MPI boundary cells)
833  unsigned int partition = 0, start_up = 0, counter = 0;
834  unsigned int start_nonboundary = numbers::invalid_unsigned_int;
835  bool work = true;
836  if (task_info.use_coloring_only == false)
837  {
838  start_nonboundary =
839  std::min(((size_info.boundary_cells_end+task_info.block_size-1)/
840  task_info.block_size)*task_info.block_size,
841  size_info.n_macro_cells);
842  start_up = start_nonboundary;
843  size_info.boundary_cells_end = start_nonboundary;
844  }
845  else
846  {
847  start_nonboundary = size_info.n_macro_cells;
848  start_up = size_info.n_macro_cells;
849  size_info.boundary_cells_start = 0;
850  size_info.boundary_cells_end = size_info.n_macro_cells;
851  }
852  if (hp_bool == true)
853  {
854  irregular_cells.resize (0);
855  irregular_cells.resize (size_info.n_macro_cells+2*max_fe_index);
856  std::vector<std::vector<unsigned int> > renumbering_fe_index;
857  renumbering_fe_index.resize(max_fe_index);
858  unsigned int counter,n_macro_cells_before = 0;
859  for (counter=0; counter<start_nonboundary*vectorization_length;
860  counter++)
861  {
862  renumbering_fe_index[cell_active_fe_index[renumbering[counter]]].
863  push_back(renumbering[counter]);
864  }
865  counter = 0;
866  for (unsigned int j=0; j<max_fe_index; j++)
867  {
868  for (unsigned int jj=0; jj<renumbering_fe_index[j].size(); jj++)
869  renumbering[counter++] = renumbering_fe_index[j][jj];
870  irregular_cells[renumbering_fe_index[j].size()/vectorization_length+
871  n_macro_cells_before] =
872  renumbering_fe_index[j].size()%vectorization_length;
873  n_macro_cells_before += (renumbering_fe_index[j].size()+vectorization_length-1)/
874  vectorization_length;
875  renumbering_fe_index[j].resize(0);
876  }
877 
878  unsigned int new_boundary_end = n_macro_cells_before;
879  for (counter=start_nonboundary*vectorization_length;
880  counter<size_info.n_active_cells; counter++)
881  {
882  renumbering_fe_index[cell_active_fe_index.empty() ? 0 :
883  cell_active_fe_index[renumbering[counter]]].
884  push_back(renumbering[counter]);
885  }
886  counter = start_nonboundary * vectorization_length;
887  for (unsigned int j=0; j<max_fe_index; j++)
888  {
889  for (unsigned int jj=0; jj<renumbering_fe_index[j].size(); jj++)
890  renumbering[counter++] = renumbering_fe_index[j][jj];
891  irregular_cells[renumbering_fe_index[j].size()/vectorization_length+
892  n_macro_cells_before] =
893  renumbering_fe_index[j].size()%vectorization_length;
894  n_macro_cells_before += (renumbering_fe_index[j].size()+vectorization_length-1)/
895  vectorization_length;
896  }
897  AssertIndexRange (n_macro_cells_before,
898  size_info.n_macro_cells + 2*max_fe_index+1);
899  irregular_cells.resize (n_macro_cells_before);
900  size_info.n_macro_cells = n_macro_cells_before;
901  size_info.boundary_cells_start = 0;
902  size_info.boundary_cells_end = new_boundary_end;
903  task_info.n_blocks = (size_info.n_macro_cells+task_info.block_size-1)
904  /task_info.block_size;
905  task_info.block_size_last = size_info.n_macro_cells%task_info.block_size;
906  if (task_info.block_size_last == 0)
907  task_info.block_size_last = task_info.block_size;
908  }
909 
910  // assume that all FEs have the same connectivity graph, so take the
911  // zeroth FE
912  task_info.n_blocks = (size_info.n_macro_cells+task_info.block_size-1)/
913  task_info.block_size;
914  task_info.block_size_last = size_info.n_macro_cells-
915  (task_info.block_size*(task_info.n_blocks-1));
916 
917  // create the connectivity graph with internal blocking
918  CompressedSimpleSparsityPattern connectivity;
919  make_connectivity_graph (size_info, task_info, renumbering,irregular_cells,
920  true, connectivity);
921 
922  // Create cell-block partitioning.
923 
924  // For each block of cells, this variable saves to which partitions the
925  // block belongs. Initialize all to n_macro_cells to mark them as not
926  // yet assigned a partition.
927  std::vector<unsigned int> cell_partition(task_info.n_blocks,
928  size_info.n_macro_cells);
929  std::vector<unsigned int> neighbor_list;
930  std::vector<unsigned int> neighbor_neighbor_list;
931 
932  // In element j of this variable, one puts the old number of the block
933  // that should be the jth block in the new numeration.
934  std::vector<unsigned int> partition_list (task_info.n_blocks,0);
935  std::vector<unsigned int> partition_color_list(task_info.n_blocks,0);
936 
937  // This vector points to the start of each partition.
938  std::vector<unsigned int> partition_blocks (2,0);
939  std::vector<unsigned int> cell_color(task_info.n_blocks,
940  size_info.n_macro_cells);
941  std::vector<bool> color_finder;
942 
943  // this performs a classical breath-first search in the connectivity
944  // graph of the cell chunks
945  while (work)
946  {
947  // put all cells up to begin_inner_cells into first partition. if
948  // the numbers do not add up exactly, assign an additional block
949  if (start_nonboundary>0 && start_up == start_nonboundary)
950  {
951  unsigned int n_blocks = ((start_nonboundary+task_info.block_size-1)
952  /task_info.block_size);
953  for (unsigned int cell=0; cell<n_blocks; ++cell)
954  {
955  cell_partition[cell] = partition;
956  neighbor_list.push_back(cell);
957  partition_list[counter++] = cell;
958  partition_blocks.back()++;
959  }
960  }
961  else
962  {
963  // To start up, set the start_up cell to partition and list all
964  // its neighbors.
965  AssertIndexRange(start_up, cell_partition.size());
966  cell_partition[start_up] = partition;
967  neighbor_list.push_back(start_up);
968  partition_list[counter++] = start_up;
969  partition_blocks.back()++;
970  }
971 
972  while (neighbor_list.size()>0)
973  {
974  partition++;
975  partition_blocks.push_back(partition_blocks.back());
976  for (unsigned int j=0; j<neighbor_list.size(); ++j)
977  {
978  Assert(cell_partition[neighbor_list[j]]==partition-1,
979  ExcInternalError());
981  connectivity.row_begin(neighbor_list[j]),
982  end = connectivity.row_end(neighbor_list[j]);
983  for (; neighbor!=end ; ++neighbor)
984  {
985  if (cell_partition[*neighbor]==size_info.n_macro_cells)
986  {
987  partition_blocks.back()++;
988  cell_partition[*neighbor] = partition;
989  neighbor_neighbor_list.push_back(*neighbor);
990  partition_list[counter++] = *neighbor;
991  }
992  }
993  }
994  neighbor_list = neighbor_neighbor_list;
995  neighbor_neighbor_list.resize(0);
996  }
997 
998  // One has to check if the graph is not connected so we have to find
999  // another partition.
1000  work = false;
1001  for (unsigned int j=start_up; j<task_info.n_blocks; ++j)
1002  if (cell_partition[j] == size_info.n_macro_cells)
1003  {
1004  start_up = j;
1005  work = true;
1006  break;
1007  }
1008  }
1009  AssertDimension (partition_blocks[partition], task_info.n_blocks);
1010 
1011 
1012  // Color the cells within each partition
1013  task_info.partition_color_blocks_row_index.resize(partition+1);
1014  unsigned int color_counter = 0, index_counter = 0;
1015  for (unsigned int part=0; part<partition; part++)
1016  {
1017  task_info.partition_color_blocks_row_index[part] = index_counter;
1018  unsigned int max_color = 0;
1019  for (unsigned int k=partition_blocks[part]; k<partition_blocks[part+1];
1020  k++)
1021  {
1022  unsigned int cell = partition_list[k];
1023  unsigned int n_neighbors = connectivity.row_length(cell);
1024 
1025  // In the worst case, each neighbor has a different color. So we
1026  // find at least one available color between 0 and n_neighbors.
1027  color_finder.resize(n_neighbors+1);
1028  for (unsigned int j=0; j<=n_neighbors; ++j)
1029  color_finder[j]=true;
1031  neighbor = connectivity.row_begin(cell),
1032  end = connectivity.row_end(cell);
1033  for (; neighbor!=end ; ++neighbor)
1034  {
1035  // Mark the color that a neighbor within the partition has
1036  // as taken
1037  if (cell_partition[*neighbor] == part &&
1038  cell_color[*neighbor] <= n_neighbors)
1039  color_finder[cell_color[*neighbor]] = false;
1040  }
1041  // Choose the smallest color that is not taken for the block
1042  cell_color[cell]=0;
1043  while (color_finder[cell_color[cell]] == false)
1044  cell_color[cell]++;
1045  if (cell_color[cell] > max_color)
1046  max_color = cell_color[cell];
1047  }
1048  // Reorder within partition: First, all blocks that belong the 0 and
1049  // then so on until those with color max (Note that the smaller the
1050  // number the larger the partition)
1051  for (unsigned int color=0; color<=max_color; color++)
1052  {
1053  task_info.partition_color_blocks_data.push_back(color_counter);
1054  index_counter++;
1055  for (unsigned int k=partition_blocks[part];
1056  k<partition_blocks[part+1]; k++)
1057  {
1058  unsigned int cell=partition_list[k];
1059  if (cell_color[cell] == color)
1060  {
1061  partition_color_list[color_counter++] = cell;
1062  }
1063  }
1064  }
1065  }
1066  task_info.partition_color_blocks_data.push_back(task_info.n_blocks);
1067  task_info.partition_color_blocks_row_index[partition] = index_counter;
1068  AssertDimension (color_counter, task_info.n_blocks);
1069 
1070  partition_list = renumbering;
1071 
1072  // in debug mode, check that the partition color list is one-to-one
1073 #ifdef DEBUG
1074  {
1075  std::vector<unsigned int> sorted_pc_list (partition_color_list);
1076  std::sort(sorted_pc_list.begin(), sorted_pc_list.end());
1077  for (unsigned int i=0; i<sorted_pc_list.size(); ++i)
1078  Assert(sorted_pc_list[i] == i, ExcInternalError());
1079  }
1080 #endif
1081 
1082  // set the start list for each block and compute the renumbering of
1083  // cells
1084  std::vector<unsigned int> block_start(size_info.n_macro_cells+1);
1085  std::vector<unsigned int> irregular(size_info.n_macro_cells);
1086 
1087  unsigned int mcell_start=0;
1088  block_start[0] = 0;
1089  for (unsigned int block=0; block<task_info.n_blocks; block++)
1090  {
1091  block_start[block+1] = block_start[block];
1092  for (unsigned int mcell=mcell_start; mcell<
1093  std::min(mcell_start+task_info.block_size,
1094  size_info.n_macro_cells);
1095  ++mcell)
1096  {
1097  unsigned int n_comp = (irregular_cells[mcell]>0)
1098  ?irregular_cells[mcell]:size_info.vectorization_length;
1099  block_start[block+1] += n_comp;
1100  ++counter;
1101  }
1102  mcell_start += task_info.block_size;
1103  }
1104  counter = 0;
1105  unsigned int counter_macro = 0;
1106  for (unsigned int block=0; block<task_info.n_blocks; block++)
1107  {
1108  unsigned int present_block = partition_color_list[block];
1109  for (unsigned int cell = block_start[present_block];
1110  cell<block_start[present_block+1]; ++cell)
1111  renumbering[counter++] = partition_list[cell];
1112  unsigned int this_block_size = (present_block == task_info.n_blocks-1)?
1113  task_info.block_size_last:task_info.block_size;
1114  for (unsigned int j=0; j<this_block_size; j++)
1115  irregular[counter_macro++] =
1116  irregular_cells[present_block*task_info.block_size+j];
1117  if (present_block == task_info.n_blocks-1)
1118  task_info.position_short_block = block;
1119  }
1120  irregular_cells.swap(irregular);
1121  AssertDimension (counter, size_info.n_active_cells);
1122  AssertDimension (counter_macro, size_info.n_macro_cells);
1123 
1124  // check that the renumbering is one-to-one
1125 #ifdef DEBUG
1126  {
1127  std::vector<unsigned int> sorted_renumbering (renumbering);
1128  std::sort(sorted_renumbering.begin(), sorted_renumbering.end());
1129  for (unsigned int i=0; i<sorted_renumbering.size(); ++i)
1130  Assert(sorted_renumbering[i] == i, ExcInternalError());
1131  }
1132 #endif
1133  AssertDimension(counter,size_info.n_active_cells);
1134  task_info.evens = (partition+1)>>1;
1135  task_info.odds = (partition)>>1;
1136  task_info.n_blocked_workers = task_info.odds-
1137  (task_info.odds+task_info.evens+1)%2;
1138  task_info.n_workers = task_info.partition_color_blocks_data.size()-1-
1139  task_info.n_blocked_workers;
1140  }
1141 
1142 
1143 
1144  void
1145  DoFInfo::make_thread_graph_partition_partition
1146  (SizeInfo &size_info,
1147  TaskInfo &task_info,
1148  std::vector<unsigned int> &renumbering,
1149  std::vector<unsigned int> &irregular_cells,
1150  const bool hp_bool)
1151  {
1152  if (size_info.n_macro_cells == 0)
1153  return;
1154 
1155  const std::size_t vectorization_length = size_info.vectorization_length;
1156  Assert (vectorization_length > 0, ExcInternalError());
1157 
1158  guess_block_size (size_info, task_info);
1159 
1160  // assume that all FEs have the same connectivity graph, so take the
1161  // zeroth FE
1162  task_info.n_blocks = (size_info.n_macro_cells+task_info.block_size-1)/
1163  task_info.block_size;
1164  task_info.block_size_last = size_info.n_macro_cells-
1165  (task_info.block_size*(task_info.n_blocks-1));
1166  task_info.position_short_block = task_info.n_blocks-1;
1167  unsigned int cluster_size = task_info.block_size*vectorization_length;
1168 
1169  // create the connectivity graph without internal blocking
1170  CompressedSimpleSparsityPattern connectivity;
1171  make_connectivity_graph (size_info, task_info, renumbering,irregular_cells,
1172  false, connectivity);
1173 
1174  // Create cell-block partitioning.
1175 
1176  // For each block of cells, this variable saves to which partitions the
1177  // block belongs. Initialize all to n_macro_cells to mark them as not
1178  // yet assigned a partition.
1179  std::vector<unsigned int> cell_partition (size_info.n_active_cells,
1180  size_info.n_active_cells);
1181  std::vector<unsigned int> neighbor_list;
1182  std::vector<unsigned int> neighbor_neighbor_list;
1183 
1184  // In element j of this variable, one puts the old number of the block
1185  // that should be the jth block in the new numeration.
1186  std::vector<unsigned int> partition_list(size_info.n_active_cells,0);
1187  std::vector<unsigned int> partition_partition_list(size_info.n_active_cells,0);
1188 
1189  // This vector points to the start of each partition.
1190  std::vector<unsigned int> partition_size(2,0);
1191 
1192  unsigned int partition = 0,start_up=0,counter=0;
1193  unsigned int start_nonboundary = vectorization_length * size_info.boundary_cells_end;
1194  if (start_nonboundary > size_info.n_active_cells)
1195  start_nonboundary = size_info.n_active_cells;
1196  bool work = true;
1197  unsigned int remainder = cluster_size;
1198 
1199  // this performs a classical breath-first search in the connectivity
1200  // graph of the cells under the restriction that the size of the
1201  // partitions should be a multiple of the given block size
1202  while (work)
1203  {
1204  // put the cells with neighbors on remote MPI processes up front
1205  if (start_nonboundary>0)
1206  {
1207  for (unsigned int cell=0; cell<start_nonboundary; ++cell)
1208  {
1209  const unsigned int cell_nn = renumbering[cell];
1210  cell_partition[cell_nn] = partition;
1211  neighbor_list.push_back(cell_nn);
1212  partition_list[counter++] = cell_nn;
1213  partition_size.back()++;
1214  }
1215  remainder -= (start_nonboundary%cluster_size);
1216  if (remainder == cluster_size)
1217  remainder = 0;
1218 
1219  // adjust end of boundary cells to the remainder
1220  size_info.boundary_cells_end += (remainder+vectorization_length-1)/vectorization_length;
1221  }
1222  else
1223  {
1224  // To start up, set the start_up cell to partition and list all
1225  // its neighbors.
1226  cell_partition[start_up] = partition;
1227  neighbor_list.push_back(start_up);
1228  partition_list[counter++] = start_up;
1229  partition_size.back()++;
1230  start_up++;
1231  remainder--;
1232  if (remainder == cluster_size)
1233  remainder = 0;
1234  }
1235  int index_before = neighbor_list.size(), index = index_before,
1236  index_stop = 0;
1237  while (remainder>0)
1238  {
1239  if (index==index_stop)
1240  {
1241  index = neighbor_list.size();
1242  if (index == index_before)
1243  {
1244  neighbor_list.resize(0);
1245  goto not_connect;
1246  }
1247  index_stop = index_before;
1248  index_before = index;
1249  }
1250  index--;
1251  unsigned int additional = neighbor_list[index];
1253  connectivity.row_begin(additional),
1254  end = connectivity.row_end(additional);
1255  for (; neighbor!=end ; ++neighbor)
1256  {
1257  if (cell_partition[*neighbor]==size_info.n_active_cells)
1258  {
1259  partition_size.back()++;
1260  cell_partition[*neighbor] = partition;
1261  neighbor_list.push_back(*neighbor);
1262  partition_list[counter++] = *neighbor;
1263  remainder--;
1264  if (remainder == 0)
1265  break;
1266  }
1267  }
1268  }
1269 
1270  while (neighbor_list.size()>0)
1271  {
1272  partition++;
1273  unsigned int partition_counter = 0;
1274  partition_size.push_back(partition_size.back());
1275 
1276  for (unsigned int j=0; j<neighbor_list.size(); ++j)
1277  {
1278  Assert(cell_partition[neighbor_list[j]]==partition-1,
1279  ExcInternalError());
1281  connectivity.row_begin(neighbor_list[j]),
1282  end = connectivity.row_end(neighbor_list[j]);
1283  for (; neighbor!=end ; ++neighbor)
1284  {
1285  if (cell_partition[*neighbor]==size_info.n_active_cells)
1286  {
1287  partition_size.back()++;
1288  cell_partition[*neighbor] = partition;
1289  neighbor_neighbor_list.push_back(*neighbor);
1290  partition_list[counter++] = *neighbor;
1291  partition_counter++;
1292  }
1293  }
1294  }
1295  remainder = cluster_size-(partition_counter%cluster_size);
1296  if (remainder == cluster_size)
1297  remainder = 0;
1298  int index_stop = 0;
1299  int index_before = neighbor_neighbor_list.size(), index = index_before;
1300  while (remainder>0)
1301  {
1302  if (index==index_stop)
1303  {
1304  index = neighbor_neighbor_list.size();
1305  if (index == index_before)
1306  {
1307  neighbor_neighbor_list.resize(0);
1308  break;
1309  }
1310  index_stop = index_before;
1311  index_before = index;
1312  }
1313  index--;
1314  unsigned int additional = neighbor_neighbor_list[index];
1316  connectivity.row_begin(additional),
1317  end = connectivity.row_end(additional);
1318  for (; neighbor!=end ; ++neighbor)
1319  {
1320  if (cell_partition[*neighbor]==size_info.n_active_cells)
1321  {
1322  partition_size.back()++;
1323  cell_partition[*neighbor] = partition;
1324  neighbor_neighbor_list.push_back(*neighbor);
1325  partition_list[counter++] = *neighbor;
1326  remainder--;
1327  if (remainder == 0)
1328  break;
1329  }
1330  }
1331  }
1332 
1333  neighbor_list = neighbor_neighbor_list;
1334  neighbor_neighbor_list.resize(0);
1335  }
1336 not_connect:
1337  // One has to check if the graph is not connected so we have to find
1338  // another partition.
1339  work = false;
1340  for (unsigned int j=start_up; j<size_info.n_active_cells; ++j)
1341  if (cell_partition[j] == size_info.n_active_cells)
1342  {
1343  start_up = j;
1344  work = true;
1345  if (remainder == 0)
1346  remainder = cluster_size;
1347  break;
1348  }
1349  }
1350  if (remainder != 0)
1351  partition++;
1352 
1353  for (unsigned int j=0; j<renumbering.size(); j++)
1354  renumbering[j] = 0;
1355  irregular_cells.back() = 0;
1356  irregular_cells.resize(size_info.n_active_cells);
1357  unsigned int n_macro_cells_before = 0;
1358  {
1359  // Create partitioning within partitions.
1360 
1361  // For each block of cells, this variable saves to which partitions
1362  // the block belongs. Initialize all to n_macro_cells to mark them as
1363  // not yet assigned a partition.
1364  std::vector<unsigned int> cell_partition_l2(size_info.n_active_cells,
1365  size_info.n_active_cells);
1366  task_info.partition_color_blocks_row_index.resize(partition+1,0);
1367  task_info.partition_color_blocks_data.resize(1,0);
1368 
1369  start_up = 0;
1370  counter = 0;
1371  unsigned int missing_macros;
1372  for (unsigned int part=0; part<partition; ++part)
1373  {
1374  neighbor_neighbor_list.resize(0);
1375  neighbor_list.resize(0);
1376  bool work = true;
1377  unsigned int partition_l2 = 0;
1378  start_up = partition_size[part];
1379  unsigned int partition_counter = 0;
1380  while (work)
1381  {
1382  if (neighbor_list.size()==0)
1383  {
1384  work = false;
1385  partition_counter = 0;
1386  for (unsigned int j=start_up; j<partition_size[part+1]; ++j)
1387  if (cell_partition[partition_list[j]] == part &&
1388  cell_partition_l2[partition_list[j]] == size_info.n_active_cells)
1389  {
1390  start_up = j;
1391  work = true;
1392  partition_counter = 1;
1393  // To start up, set the start_up cell to partition
1394  // and list all its neighbors.
1395  AssertIndexRange (start_up, partition_size[part+1]);
1396  cell_partition_l2[partition_list[start_up]] =
1397  partition_l2;
1398  neighbor_neighbor_list.push_back
1399  (partition_list[start_up]);
1400  partition_partition_list[counter++] =
1401  partition_list[start_up];
1402  start_up++;
1403  break;
1404  }
1405  }
1406  else
1407  {
1408  partition_counter = 0;
1409  for (unsigned int j=0; j<neighbor_list.size(); ++j)
1410  {
1411  Assert(cell_partition[neighbor_list[j]]==part,
1412  ExcInternalError());
1413  Assert(cell_partition_l2[neighbor_list[j]]==partition_l2-1,
1414  ExcInternalError());
1416  connectivity.row_begin(neighbor_list[j]),
1417  end = connectivity.row_end(neighbor_list[j]);
1418  for (; neighbor!=end ; ++neighbor)
1419  {
1420  if (cell_partition[*neighbor] == part &&
1421  cell_partition_l2[*neighbor]==
1422  size_info.n_active_cells)
1423  {
1424  cell_partition_l2[*neighbor] = partition_l2;
1425  neighbor_neighbor_list.push_back(*neighbor);
1426  partition_partition_list[counter++] = *neighbor;
1427  partition_counter++;
1428  }
1429  }
1430  }
1431  }
1432  if (partition_counter>0)
1433  {
1434  int index_before = neighbor_neighbor_list.size(),
1435  index = index_before;
1436  {
1437  // put the cells into separate lists for each FE index
1438  // within one partition-partition
1439  missing_macros = 0;
1440  std::vector<unsigned int> remaining_per_macro_cell
1441  (max_fe_index);
1442  std::vector<std::vector<unsigned int> >
1443  renumbering_fe_index;
1444  unsigned int cell;
1445  bool filled = true;
1446  if (hp_bool == true)
1447  {
1448  renumbering_fe_index.resize(max_fe_index);
1449  for (cell=counter-partition_counter; cell<counter; ++cell)
1450  {
1451  renumbering_fe_index
1452  [cell_active_fe_index.empty() ? 0 :
1453  cell_active_fe_index[partition_partition_list
1454  [cell]]].
1455  push_back(partition_partition_list[cell]);
1456  }
1457  // check how many more cells are needed in the lists
1458  for (unsigned int j=0; j<max_fe_index; j++)
1459  {
1460  remaining_per_macro_cell[j] =
1461  renumbering_fe_index[j].size()%vectorization_length;
1462  if (remaining_per_macro_cell[j] != 0)
1463  filled = false;
1464  missing_macros += ((renumbering_fe_index[j].size()+
1465  vectorization_length-1)/vectorization_length);
1466  }
1467  }
1468  else
1469  {
1470  remaining_per_macro_cell.resize(1);
1471  remaining_per_macro_cell[0] = partition_counter%
1472  vectorization_length;
1473  missing_macros = partition_counter/vectorization_length;
1474  if (remaining_per_macro_cell[0] != 0)
1475  {
1476  filled = false;
1477  missing_macros++;
1478  }
1479  }
1480  missing_macros = task_info.block_size -
1481  (missing_macros%task_info.block_size);
1482 
1483  // now we realized that there are some cells missing.
1484  while (missing_macros>0 || filled == false)
1485  {
1486  if (index==0)
1487  {
1488  index = neighbor_neighbor_list.size();
1489  if (index == index_before)
1490  {
1491  if (missing_macros != 0)
1492  {
1493  neighbor_neighbor_list.resize(0);
1494  }
1495  start_up--;
1496  break;// not connected - start again
1497  }
1498  index_before = index;
1499  }
1500  index--;
1501  unsigned int additional = neighbor_neighbor_list
1502  [index];
1503 
1504  // go through the neighbors of the last cell in the
1505  // current partition and check if we find some to
1506  // fill up with.
1508  neighbor =
1509  connectivity.row_begin(additional),
1510  end = connectivity.row_end(additional);
1511  for (; neighbor!=end ; ++neighbor)
1512  {
1513  if (cell_partition[*neighbor] == part &&
1514  cell_partition_l2[*neighbor] ==
1515  size_info.n_active_cells)
1516  {
1517  unsigned int this_index = 0;
1518  if (hp_bool == true)
1519  this_index = cell_active_fe_index.empty() ? 0 :
1520  cell_active_fe_index[*neighbor];
1521 
1522  // Only add this cell if we need more macro
1523  // cells in the current block or if there is
1524  // a macro cell with the FE index that is
1525  // not yet fully populated
1526  if (missing_macros > 0 ||
1527  remaining_per_macro_cell[this_index] > 0)
1528  {
1529  cell_partition_l2[*neighbor] = partition_l2;
1530  neighbor_neighbor_list.push_back(*neighbor);
1531  if (hp_bool == true)
1532  renumbering_fe_index[this_index].
1533  push_back(*neighbor);
1534  partition_partition_list[counter] =
1535  *neighbor;
1536  counter++;
1537  partition_counter++;
1538  if (remaining_per_macro_cell[this_index]
1539  == 0 && missing_macros > 0)
1540  missing_macros--;
1541  remaining_per_macro_cell[this_index]++;
1542  if (remaining_per_macro_cell[this_index]
1543  == vectorization_length)
1544  {
1545  remaining_per_macro_cell[this_index] = 0;
1546  }
1547  if (missing_macros == 0)
1548  {
1549  filled = true;
1550  for (unsigned int fe_ind=0;
1551  fe_ind<max_fe_index; ++fe_ind)
1552  if (remaining_per_macro_cell[fe_ind]!=0)
1553  filled = false;
1554  }
1555  if (filled == true)
1556  break;
1557  }
1558  }
1559  }
1560  }
1561  if (hp_bool == true)
1562  {
1563  // set the renumbering according to their active FE
1564  // index within one partition-partition which was
1565  // implicitly assumed above
1566  cell = counter - partition_counter;
1567  for (unsigned int j=0; j<max_fe_index; j++)
1568  {
1569  for (unsigned int jj=0; jj<renumbering_fe_index[j].
1570  size(); jj++)
1571  renumbering[cell++] =
1572  renumbering_fe_index[j][jj];
1573  if (renumbering_fe_index[j].size()%vectorization_length != 0)
1574  irregular_cells[renumbering_fe_index[j].size()/
1575  vectorization_length+
1576  n_macro_cells_before] =
1577  renumbering_fe_index[j].size()%vectorization_length;
1578  n_macro_cells_before += (renumbering_fe_index[j].
1579  size()+vectorization_length-1)/
1580  vectorization_length;
1581  renumbering_fe_index[j].resize(0);
1582  }
1583  }
1584  else
1585  {
1586  n_macro_cells_before += partition_counter/vectorization_length;
1587  if (partition_counter%vectorization_length != 0)
1588  {
1589  irregular_cells[n_macro_cells_before] =
1590  partition_counter%vectorization_length;
1591  n_macro_cells_before++;
1592  }
1593  }
1594  }
1595  task_info.partition_color_blocks_data.
1596  push_back(n_macro_cells_before);
1597  partition_l2++;
1598  }
1599  neighbor_list = neighbor_neighbor_list;
1600  neighbor_neighbor_list.resize(0);
1601  }
1602  task_info.partition_color_blocks_row_index[part+1] =
1603  task_info.partition_color_blocks_row_index[part] + partition_l2;
1604  }
1605  }
1606 
1607  if (size_info.boundary_cells_end>0)
1608  size_info.boundary_cells_end = task_info.partition_color_blocks_data
1609  [task_info.partition_color_blocks_row_index[1]];
1610 
1611  if (hp_bool == false)
1612  renumbering.swap(partition_partition_list);
1613  irregular_cells.resize(n_macro_cells_before);
1614  size_info.n_macro_cells = n_macro_cells_before;
1615 
1616  task_info.evens = (partition+1)/2;
1617  task_info.odds = partition/2;
1618  task_info.n_blocked_workers =
1619  task_info.odds-(task_info.odds+task_info.evens+1)%2;
1620  task_info.n_workers = task_info.evens+task_info.odds-
1621  task_info.n_blocked_workers;
1622  task_info.partition_evens.resize(partition);
1623  task_info.partition_odds.resize(partition);
1624  task_info.partition_n_blocked_workers.resize(partition);
1625  task_info.partition_n_workers.resize(partition);
1626  for (unsigned int part=0; part<partition; part++)
1627  {
1628  task_info.partition_evens[part] =
1629  (task_info.partition_color_blocks_row_index[part+1]-
1630  task_info.partition_color_blocks_row_index[part]+1)/2;
1631  task_info.partition_odds[part] =
1632  (task_info.partition_color_blocks_row_index[part+1]-
1633  task_info.partition_color_blocks_row_index[part])/2;
1634  task_info.partition_n_blocked_workers[part] =
1635  task_info.partition_odds[part]-(task_info.partition_odds[part]+
1636  task_info.partition_evens[part]+1)%2;
1637  task_info.partition_n_workers[part] =
1638  task_info.partition_evens[part]+task_info.partition_odds[part]-
1639  task_info.partition_n_blocked_workers[part];
1640  }
1641  }
1642 
1643 
1644  namespace internal
1645  {
1646  // rudimentary version of a vector that keeps entries always ordered
1647  class ordered_vector : public std::vector<types::global_dof_index>
1648  {
1649  public:
1650  ordered_vector ()
1651  {
1652  reserve (2000);
1653  }
1654 
1655  void reserve (const std::size_t size)
1656  {
1657  if (size > 0)
1658  this->std::vector<types::global_dof_index>::reserve (size);
1659  }
1660 
1661 
1662  // insert a given entry. dat is a pointer within this vector (the user
1663  // needs to make sure that it really stays there)
1664  void insert (const unsigned int entry,
1665  std::vector<types::global_dof_index>::iterator &dat)
1666  {
1667  AssertIndexRange (static_cast<std::size_t>(dat - begin()), size()+1);
1668  AssertIndexRange (static_cast<std::size_t>(end() - dat), size()+1);
1669  AssertIndexRange (size(), capacity());
1670  while (dat != end() && *dat < entry)
1671  ++dat;
1672 
1673  if (dat == end())
1674  {
1675  push_back(entry);
1676  dat = end();
1677  }
1678  else if (*dat > entry)
1679  {
1680  dat = this->std::vector<types::global_dof_index>::insert (dat, entry);
1681  ++dat;
1682  }
1683  else
1684  ++dat;
1685  }
1686  };
1687  }
1688 
1689 
1690  void
1691  DoFInfo::make_connectivity_graph
1692  (const SizeInfo &size_info,
1693  const TaskInfo &task_info,
1694  const std::vector<unsigned int> &renumbering,
1695  const std::vector<unsigned int> &irregular_cells,
1696  const bool do_blocking,
1697  CompressedSimpleSparsityPattern &connectivity) const
1698  {
1699  AssertDimension (row_starts.size()-1, size_info.n_active_cells);
1700  const unsigned int n_rows =
1701  (vector_partitioner->local_range().second-
1702  vector_partitioner->local_range().first)
1703  + vector_partitioner->ghost_indices().n_elements();
1704  const unsigned int n_blocks = (do_blocking == true) ?
1705  task_info.n_blocks : size_info.n_active_cells;
1706 
1707  // first determine row lengths
1708  std::vector<unsigned int> row_lengths(n_rows);
1709  unsigned int cell_start = 0, mcell_start = 0;
1710  std::vector<unsigned int> scratch;
1711  for (unsigned int block = 0; block < n_blocks; ++block)
1712  {
1713  // if we have the blocking variant (used in the coloring scheme), we
1714  // want to build a graph with the blocks with interaction with
1715  // remote MPI processes up front. in the non-blocking variant, we do
1716  // not do this here. TODO: unify this approach!!!
1717  if (do_blocking == true)
1718  {
1719  scratch.clear();
1720  for (unsigned int mcell=mcell_start; mcell<
1721  std::min(mcell_start+task_info.block_size,
1722  size_info.n_macro_cells);
1723  ++mcell)
1724  {
1725  unsigned int n_comp = (irregular_cells[mcell]>0)
1726  ?irregular_cells[mcell]:size_info.vectorization_length;
1727  for (unsigned int cell = cell_start; cell < cell_start+n_comp;
1728  ++cell)
1729  scratch.insert(scratch.end(),
1730  begin_indices(renumbering[cell]),
1731  end_indices(renumbering[cell]));
1732  cell_start += n_comp;
1733  }
1734  std::sort(scratch.begin(), scratch.end());
1735  const unsigned int n_unique =
1736  std::unique(scratch.begin(), scratch.end())-scratch.begin();
1737  for (unsigned int i=0; i<n_unique; ++i)
1738  row_lengths[scratch[i]]++;
1739  mcell_start += task_info.block_size;
1740  }
1741  else
1742  {
1743  scratch.clear();
1744  scratch.insert(scratch.end(),
1745  begin_indices(block), end_indices(block));
1746  std::sort(scratch.begin(), scratch.end());
1747  const unsigned int n_unique =
1748  std::unique(scratch.begin(), scratch.end())-scratch.begin();
1749  for (unsigned int i=0; i<n_unique; ++i)
1750  row_lengths[scratch[i]]++;
1751  }
1752  }
1753 
1754  // disregard dofs that only sit on one cell
1755  for (unsigned int row=0; row<n_rows; ++row)
1756  if (row_lengths[row] == 1)
1757  row_lengths[row] = 0;
1758 
1759  SparsityPattern connectivity_dof (n_rows, n_blocks, row_lengths);
1760  cell_start = 0, mcell_start = 0;
1761  for (unsigned int block = 0; block < n_blocks; ++block)
1762  {
1763  // if we have the blocking variant (used in the coloring scheme), we
1764  // want to build a graph with the blocks with interaction with
1765  // remote MPI processes up front. in the non-blocking variant, we do
1766  // not do this here. TODO: unify this approach!!!
1767  if (do_blocking == true)
1768  {
1769  for (unsigned int mcell=mcell_start; mcell<
1770  std::min(mcell_start+task_info.block_size,
1771  size_info.n_macro_cells);
1772  ++mcell)
1773  {
1774  unsigned int n_comp = (irregular_cells[mcell]>0)
1775  ?irregular_cells[mcell]:size_info.vectorization_length;
1776  for (unsigned int cell = cell_start; cell < cell_start+n_comp;
1777  ++cell)
1778  {
1779  const unsigned int
1780  *it = begin_indices (renumbering[cell]),
1781  *end_cell = end_indices (renumbering[cell]);
1782  for ( ; it != end_cell; ++it)
1783  if (row_lengths[*it]>0)
1784  connectivity_dof.add(*it, block);
1785  }
1786  cell_start += n_comp;
1787  }
1788  mcell_start += task_info.block_size;
1789  }
1790  else
1791  {
1792  const unsigned int
1793  *it = begin_indices (block),
1794  *end_cell = end_indices (block);
1795  for ( ; it != end_cell; ++it)
1796  if (row_lengths[*it]>0)
1797  connectivity_dof.add(*it, block);
1798  }
1799  }
1800  connectivity_dof.compress();
1801 
1802  connectivity.reinit (n_blocks, n_blocks);
1803  internal::ordered_vector row_entries;
1804  cell_start = 0;
1805  mcell_start = 0;
1806  for (unsigned int block=0; block < n_blocks; ++block)
1807  {
1808  row_entries.clear();
1809 
1810  if (do_blocking==true)
1811  {
1812  for (unsigned int mcell=mcell_start; mcell<
1813  std::min(mcell_start+task_info.block_size,
1814  size_info.n_macro_cells);
1815  ++mcell)
1816  {
1817  unsigned int n_comp = (irregular_cells[mcell]>0)
1818  ?irregular_cells[mcell]:size_info.vectorization_length;
1819  for (unsigned int cell = cell_start; cell < cell_start+n_comp;
1820  ++cell)
1821  {
1822  // apply renumbering when we do blocking
1823  const unsigned int
1824  *it = begin_indices (renumbering[cell]),
1825  *end_cell = end_indices (renumbering[cell]);
1826  for ( ; it != end_cell; ++it)
1827  if (row_lengths[*it] > 0)
1828  {
1829  SparsityPattern::iterator sp = connectivity_dof.begin(*it);
1830  // jump over diagonal for square patterns
1831  if (connectivity_dof.n_rows()==connectivity_dof.n_cols())
1832  ++sp;
1833  row_entries.reserve (row_entries.size() + end_cell - it);
1834  std::vector<types::global_dof_index>::iterator insert_pos = row_entries.begin();
1835  for ( ; sp != connectivity_dof.end(*it); ++sp)
1836  if (sp->column() >= block)
1837  break;
1838  else
1839  row_entries.insert (sp->column(), insert_pos);
1840  }
1841  }
1842  cell_start +=n_comp;
1843  }
1844  mcell_start += task_info.block_size;
1845  }
1846  else
1847  {
1848  const unsigned int *it = begin_indices (block),
1849  * end_cell = end_indices (block);
1850  for ( ; it != end_cell; ++it)
1851  if (row_lengths[*it] > 0)
1852  {
1853  SparsityPattern::iterator sp = connectivity_dof.begin(*it);
1854  // jump over diagonal for square patterns
1855  if (connectivity_dof.n_rows()==connectivity_dof.n_cols())
1856  ++sp;
1857  row_entries.reserve (row_entries.size() + end_cell - it);
1858  std::vector<types::global_dof_index>::iterator insert_pos = row_entries.begin();
1859  for ( ; sp != connectivity_dof.end(*it); ++sp)
1860  if (sp->column() >= block)
1861  break;
1862  else
1863  row_entries.insert (sp->column(), insert_pos);
1864  }
1865  }
1866  connectivity.add_entries (block, row_entries.begin(), row_entries.end());
1867  }
1868  connectivity.symmetrize ();
1869  }
1870 
1871 
1872 
1873  void DoFInfo::renumber_dofs (std::vector<types::global_dof_index> &renumbering)
1874  {
1875  // first renumber all locally owned degrees of freedom
1876  AssertDimension (vector_partitioner->local_size(),
1877  vector_partitioner->size());
1878  const unsigned int local_size = vector_partitioner->local_size();
1879  renumbering.resize (0);
1880  renumbering.resize (local_size, numbers::invalid_dof_index);
1881 
1882  types::global_dof_index counter = 0;
1883  std::vector<unsigned int>::iterator dof_ind = dof_indices.begin(),
1884  end_ind = dof_indices.end();
1885  for ( ; dof_ind != end_ind; ++dof_ind)
1886  {
1887  if (*dof_ind < local_size)
1888  {
1889  if (renumbering[*dof_ind] == numbers::invalid_dof_index)
1890  renumbering[*dof_ind] = counter++;
1891  *dof_ind = renumbering[*dof_ind];
1892  }
1893  }
1894 
1895  AssertIndexRange (counter, local_size+1);
1896  for (std::size_t i=0; i<renumbering.size(); ++i)
1897  if (renumbering[i] == numbers::invalid_dof_index)
1898  renumbering[i] = counter++;
1899 
1900  // adjust the constrained DoFs
1901  std::vector<unsigned int> new_constrained_dofs (constrained_dofs.size());
1902  for (std::size_t i=0; i<constrained_dofs.size(); ++i)
1903  new_constrained_dofs[i] = renumbering[constrained_dofs[i]];
1904 
1905  // the new constrained DoFs should be sorted already as they are not
1906  // contained in dof_indices and then get contiguous numbers
1907 #ifdef DEBUG
1908  for (std::size_t i=1; i<new_constrained_dofs.size(); ++i)
1909  Assert (new_constrained_dofs[i] > new_constrained_dofs[i-1], ExcInternalError());
1910 #endif
1911  std::swap (constrained_dofs, new_constrained_dofs);
1912 
1913  // transform indices to global index space
1914  for (std::size_t i=0; i<renumbering.size(); ++i)
1915  renumbering[i] = vector_partitioner->local_to_global(renumbering[i]);
1916 
1917  AssertDimension (counter, renumbering.size());
1918  }
1919 
1920 
1921 
1922  std::size_t
1923  DoFInfo::memory_consumption () const
1924  {
1925  std::size_t memory = sizeof(*this);
1926  memory += (row_starts.capacity()*sizeof(std_cxx1x::array<unsigned int,3>));
1927  memory += MemoryConsumption::memory_consumption (dof_indices);
1928  memory += MemoryConsumption::memory_consumption (row_starts_plain_indices);
1929  memory += MemoryConsumption::memory_consumption (plain_dof_indices);
1930  memory += MemoryConsumption::memory_consumption (constraint_indicator);
1931  memory += MemoryConsumption::memory_consumption (*vector_partitioner);
1932  return memory;
1933  }
1934 
1935 
1936 
1937  template <typename STREAM>
1938  void
1939  DoFInfo::print_memory_consumption (STREAM &out,
1940  const SizeInfo &size_info) const
1941  {
1942  out << " Memory row starts indices: ";
1943  size_info.print_memory_statistics
1944  (out, (row_starts.capacity()*sizeof(std_cxx1x::array<unsigned int, 3>)));
1945  out << " Memory dof indices: ";
1946  size_info.print_memory_statistics
1947  (out, MemoryConsumption::memory_consumption (dof_indices));
1948  out << " Memory constraint indicators: ";
1949  size_info.print_memory_statistics
1950  (out, MemoryConsumption::memory_consumption (constraint_indicator));
1951  out << " Memory plain indices: ";
1952  size_info.print_memory_statistics
1953  (out, MemoryConsumption::memory_consumption (row_starts_plain_indices)+
1954  MemoryConsumption::memory_consumption (plain_dof_indices));
1955  out << " Memory vector partitioner: ";
1956  size_info.print_memory_statistics
1957  (out, MemoryConsumption::memory_consumption (*vector_partitioner));
1958  }
1959 
1960 
1961 
1962  template <typename Number>
1963  void
1964  DoFInfo::print (const std::vector<Number> &constraint_pool_data,
1965  const std::vector<unsigned int> &constraint_pool_row_index,
1966  std::ostream &out) const
1967  {
1968  const unsigned int n_rows = row_starts.size() - 1;
1969  for (unsigned int row=0 ; row<n_rows ; ++row)
1970  {
1971  out << "Entries row " << row << ": ";
1972  const unsigned int *glob_indices = begin_indices(row),
1973  *end_row = end_indices(row);
1974  unsigned int index = 0;
1975  const std::pair<unsigned short,unsigned short>
1976  *con_it = begin_indicators(row),
1977  * end_con = end_indicators(row);
1978  for ( ; con_it != end_con; ++con_it)
1979  {
1980  for ( ; index<con_it->first; index++)
1981  {
1982  Assert (glob_indices+index != end_row, ExcInternalError());
1983  out << glob_indices[index] << " ";
1984  }
1985 
1986  out << "[ ";
1987  for (unsigned int k=constraint_pool_row_index[con_it->second];
1988  k<constraint_pool_row_index[con_it->second+1];
1989  k++,index++)
1990  {
1991  Assert (glob_indices+index != end_row, ExcInternalError());
1992  out << glob_indices[index] << "/"
1993  << constraint_pool_data[k];
1994  if (k<constraint_pool_row_index[con_it->second+1]-1)
1995  out << " ";
1996  }
1997  out << "] ";
1998  }
1999  glob_indices += index;
2000  for (; glob_indices != end_row; ++glob_indices)
2001  out << *glob_indices << " ";
2002  out << std::endl;
2003  }
2004  }
2005 
2006 
2007  } // end of namespace MatrixFreeFunctions
2008 } // end of namespace internal
2009 
2010 DEAL_II_NAMESPACE_CLOSE
static const unsigned int invalid_unsigned_int
Definition: types.h:191
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
const std::vector< std::pair< size_type, double > > * get_constraint_entries(const size_type line) const
::ExceptionBase & ExcMessage(std::string arg1)
iterator begin() const
#define AssertIndexRange(index, range)
Definition: exceptions.h:888
void print_memory_statistics(STREAM &out, std::size_t data_length) const
void add(const size_type i, const size_type j)
unsigned int global_dof_index
Definition: types.h:100
#define Assert(cond, exc)
Definition: exceptions.h:299
std::size_t memory_consumption(const T &t)
size_type n_cols() const
unsigned short insert_entries(const std::vector< std::pair< types::global_dof_index, double > > &entries)
row_iterator row_begin(const size_type row) const
iterator end() const
std::vector< unsigned int > invert_permutation(const std::vector< unsigned int > &permutation)
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_unique_and_sorted=false)
size_type row_length(const size_type row) const
row_iterator row_end(const size_type row) const
size_type n_rows() const
void reinit(const size_type m, const size_type n, const IndexSet &rowset=IndexSet())
MultithreadInfo multithread_info
const types::global_dof_index invalid_dof_index
Definition: types.h:217
::ExceptionBase & ExcInternalError()
void set_ghost_indices(const IndexSet &ghost_indices)
std::vector< size_type >::const_iterator row_iterator