libMesh
sparsity_pattern.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // Local Includes
21 #include "libmesh/coupling_matrix.h"
22 #include "libmesh/dof_map.h"
23 #include "libmesh/elem.h"
24 #include "libmesh/ghosting_functor.h"
25 #include "libmesh/sparsity_pattern.h"
26 
27 
28 namespace libMesh
29 {
30 namespace SparsityPattern
31 {
32 
33 //-------------------------------------------------------
34 // we need to implement these constructors here so that
35 // a full DofMap definition is available.
36 Build::Build (const MeshBase & mesh_in,
37  const DofMap & dof_map_in,
38  const CouplingMatrix * dof_coupling_in,
39  std::set<GhostingFunctor *> coupling_functors_in,
40  const bool implicit_neighbor_dofs_in,
41  const bool need_full_sparsity_pattern_in) :
42  ParallelObject(dof_map_in),
43  mesh(mesh_in),
44  dof_map(dof_map_in),
45  dof_coupling(dof_coupling_in),
46  coupling_functors(coupling_functors_in),
47  implicit_neighbor_dofs(implicit_neighbor_dofs_in),
48  need_full_sparsity_pattern(need_full_sparsity_pattern_in),
49  sparsity_pattern(),
50  nonlocal_pattern(),
51  n_nz(),
52  n_oz()
53 {}
54 
55 
56 
58  ParallelObject(other),
59  mesh(other.mesh),
60  dof_map(other.dof_map),
67  n_nz(),
68  n_oz()
69 {}
70 
71 
72 
73 #if defined(__GNUC__) && (__GNUC__ < 4) && !defined(__INTEL_COMPILER)
74 
75 void _dummy_function(void) {}
76 
77 #endif
78 
79 
80 
81 void Build::handle_vi_vj(const Elem * partner,
82  const std::vector<dof_id_type> & element_dofs_i,
83  unsigned int vj)
84 {
85  const unsigned int n_dofs_on_element_i =
86  cast_int<unsigned int>(element_dofs_i.size());
87 
88  const processor_id_type proc_id = mesh.processor_id();
89  const dof_id_type first_dof_on_proc = dof_map.first_dof(proc_id);
90  const dof_id_type end_dof_on_proc = dof_map.end_dof(proc_id);
91 
92  std::vector<dof_id_type>
93  element_dofs_j,
94  dofs_to_add;
95 
96  // Find element dofs for variable vj
97  dof_map.dof_indices (partner, element_dofs_j, vj);
98 #ifdef LIBMESH_ENABLE_CONSTRAINTS
99  dof_map.find_connected_dofs (element_dofs_j);
100 #endif
101 
102  // We can be more efficient if we sort the element DOFs
103  // into increasing order
104  std::sort (element_dofs_j.begin(), element_dofs_j.end());
105 
106  const unsigned int n_dofs_on_element_j =
107  cast_int<unsigned int>(element_dofs_j.size());
108 
109  // there might be 0 dofs for the other variable on the same element
110  // (when subdomain variables do not overlap) and that's when we do
111  // not do anything
112  if (n_dofs_on_element_j > 0)
113  {
114  for (unsigned int i=0; i<n_dofs_on_element_i; i++)
115  {
116  const dof_id_type ig = element_dofs_i[i];
117 
118  SparsityPattern::Row * row;
119 
120  // We save non-local row components for now so we can
121  // communicate them to other processors later.
122 
123  if ((ig >= first_dof_on_proc) &&
124  (ig < end_dof_on_proc))
125  {
126  // This is what I mean
127  // libmesh_assert_greater_equal ((ig - first_dof_on_proc), 0);
128  // but do the test like this because ig and
129  // first_dof_on_proc are unsigned ints
130  libmesh_assert_greater_equal (ig, first_dof_on_proc);
131  libmesh_assert_less (ig, (sparsity_pattern.size() +
132  first_dof_on_proc));
133 
134  row = &sparsity_pattern[ig - first_dof_on_proc];
135  }
136  else
137  {
138  row = &nonlocal_pattern[ig];
139  }
140 
141  // If the row is empty we will add *all*
142  // the element j DOFs, so just do that.
143  if (row->empty())
144  {
145  row->insert(row->end(),
146  element_dofs_j.begin(),
147  element_dofs_j.end());
148  }
149  else
150  {
151  // Build a list of the DOF indices not found in the
152  // sparsity pattern
153  dofs_to_add.clear();
154 
155  // Cache iterators. Low will move forward, subsequent
156  // searches will be on smaller ranges
157  SparsityPattern::Row::iterator
158  low = std::lower_bound
159  (row->begin(), row->end(), element_dofs_j.front()),
160  high = std::upper_bound
161  (low, row->end(), element_dofs_j.back());
162 
163  for (unsigned int j=0; j<n_dofs_on_element_j; j++)
164  {
165  const dof_id_type jg = element_dofs_j[j];
166 
167  // See if jg is in the sorted range
168  std::pair<SparsityPattern::Row::iterator,
169  SparsityPattern::Row::iterator>
170  pos = std::equal_range (low, high, jg);
171 
172  // Must add jg if it wasn't found
173  if (pos.first == pos.second)
174  dofs_to_add.push_back(jg);
175 
176  // pos.first is now a valid lower bound for any
177  // remaining element j DOFs. (That's why we sorted them.)
178  // Use it for the next search
179  low = pos.first;
180  }
181 
182  // Add to the sparsity pattern
183  if (!dofs_to_add.empty())
184  {
185  const std::size_t old_size = row->size();
186 
187  row->insert (row->end(),
188  dofs_to_add.begin(),
189  dofs_to_add.end());
190 
192  (row->begin(), row->begin()+old_size,
193  row->end());
194  }
195  }
196  } // End dofs-of-var-i loop
197  } // End if-dofs-of-var-j
198 }
199 
200 
201 
203 {
204  // Compute the sparsity structure of the global matrix. This can be
205  // fed into a PetscMatrix to allocate exactly the number of nonzeros
206  // necessary to store the matrix. This algorithm should be linear
207  // in the (# of elements)*(# nodes per element)
208  const processor_id_type proc_id = mesh.processor_id();
209  const dof_id_type n_dofs_on_proc = dof_map.n_dofs_on_processor(proc_id);
210  const dof_id_type first_dof_on_proc = dof_map.first_dof(proc_id);
211  const dof_id_type end_dof_on_proc = dof_map.end_dof(proc_id);
212 
213  sparsity_pattern.resize(n_dofs_on_proc);
214 
215  // Handle dof coupling specified by library and user coupling functors
216  {
217  const unsigned int n_var = dof_map.n_variables();
218 
219  std::vector<dof_id_type> element_dofs_i;
220 
221  std::vector<const Elem *> coupled_neighbors;
222  for (ConstElemRange::const_iterator elem_it = range.begin() ; elem_it != range.end(); ++elem_it)
223  {
224  const Elem * const elem = *elem_it;
225 
226  // Make some fake element iterators defining a range
227  // pointing to only this element.
228  Elem * const * elempp = const_cast<Elem * const *>(&elem);
229  Elem * const * elemend = elempp+1;
230 
231  const MeshBase::const_element_iterator fake_elem_it =
233  elemend,
235 
236  const MeshBase::const_element_iterator fake_elem_end =
238  elemend,
240 
241  GhostingFunctor::map_type elements_to_couple;
242 
243  // Man, I wish we had guaranteed unique_ptr availability...
244  std::set<CouplingMatrix *> temporary_coupling_matrices;
245 
246  dof_map.merge_ghost_functor_outputs(elements_to_couple,
247  temporary_coupling_matrices,
250  fake_elem_it,
251  fake_elem_end,
253 
254  for (unsigned int vi=0; vi<n_var; vi++)
255  {
256  // Find element dofs for variable vi
257  dof_map.dof_indices (elem, element_dofs_i, vi);
258 #ifdef LIBMESH_ENABLE_CONSTRAINTS
259  dof_map.find_connected_dofs (element_dofs_i);
260 #endif
261 
262  // We can be more efficient if we sort the element DOFs
263  // into increasing order
264  std::sort(element_dofs_i.begin(), element_dofs_i.end());
265 
266  GhostingFunctor::map_type::iterator etg_it = elements_to_couple.begin();
267  const GhostingFunctor::map_type::iterator etg_end = elements_to_couple.end();
268  for (; etg_it != etg_end; ++etg_it)
269  {
270  const Elem * const partner = etg_it->first;
271  const CouplingMatrix * ghost_coupling = etg_it->second;
272 
273  // Loop over coupling matrix row variables if we have a
274  // coupling matrix, or all variables if not.
275  if (ghost_coupling)
276  {
277  libmesh_assert_equal_to (ghost_coupling->size(), n_var);
278  ConstCouplingRow ccr(vi, *ghost_coupling);
279 
280  for (ConstCouplingRow::const_iterator it = ccr.begin(),
281  end = ccr.end();
282  it != end; ++it)
283  this->handle_vi_vj(partner, element_dofs_i, *it);
284  }
285  else
286  {
287  for (unsigned int vj = 0; vj != n_var; ++vj)
288  this->handle_vi_vj(partner, element_dofs_i, vj);
289  }
290  } // End ghosted element loop
291  } // End vi loop
292 
293  for (std::set<CouplingMatrix *>::iterator
294  it = temporary_coupling_matrices.begin(),
295  end = temporary_coupling_matrices.begin();
296  it != end; ++it)
297  delete *it;
298 
299  } // End range element loop
300  } // End ghosting functor section
301 
302  // Now a new chunk of sparsity structure is built for all of the
303  // DOFs connected to our rows of the matrix.
304 
305  // If we're building a full sparsity pattern, then we've got
306  // complete rows to work with, so we can just count them from
307  // scratch.
309  {
310  n_nz.clear();
311  n_oz.clear();
312  }
313 
314  n_nz.resize (n_dofs_on_proc, 0);
315  n_oz.resize (n_dofs_on_proc, 0);
316 
317  for (dof_id_type i=0; i<n_dofs_on_proc; i++)
318  {
319  // Get the row of the sparsity pattern
321 
322  for (std::size_t j=0; j<row.size(); j++)
323  if ((row[j] < first_dof_on_proc) || (row[j] >= end_dof_on_proc))
324  n_oz[i]++;
325  else
326  n_nz[i]++;
327 
328  // If we're not building a full sparsity pattern, then we want
329  // to avoid overcounting these entries as much as possible.
331  row.clear();
332  }
333 }
334 
335 
336 
338 {
339  const processor_id_type proc_id = mesh.processor_id();
340  const dof_id_type n_global_dofs = dof_map.n_dofs();
341  const dof_id_type n_dofs_on_proc = dof_map.n_dofs_on_processor(proc_id);
342  const dof_id_type first_dof_on_proc = dof_map.first_dof(proc_id);
343  const dof_id_type end_dof_on_proc = dof_map.end_dof(proc_id);
344 
345  libmesh_assert_equal_to (sparsity_pattern.size(), other.sparsity_pattern.size());
346  libmesh_assert_equal_to (n_nz.size(), sparsity_pattern.size());
347  libmesh_assert_equal_to (n_oz.size(), sparsity_pattern.size());
348 
349  for (dof_id_type r=0; r<n_dofs_on_proc; r++)
350  {
351  // increment the number of on and off-processor nonzeros in this row
352  // (note this will be an upper bound unless we need the full sparsity pattern)
354  {
356  const SparsityPattern::Row & their_row = other.sparsity_pattern[r];
357 
358  // simple copy if I have no dofs
359  if (my_row.empty())
360  my_row = their_row;
361 
362  // otherwise add their DOFs to mine, resort, and re-unique the row
363  else if (!their_row.empty()) // do nothing for the trivial case where
364  { // their row is empty
365  my_row.insert (my_row.end(),
366  their_row.begin(),
367  their_row.end());
368 
369  // We cannot use SparsityPattern::sort_row() here because it expects
370  // the [begin,middle) [middle,end) to be non-overlapping. This is not
371  // necessarily the case here, so use std::sort()
372  std::sort (my_row.begin(), my_row.end());
373 
374  my_row.erase(std::unique (my_row.begin(), my_row.end()), my_row.end());
375  }
376 
377  // fix the number of on and off-processor nonzeros in this row
378  n_nz[r] = n_oz[r] = 0;
379 
380  for (std::size_t j=0; j<my_row.size(); j++)
381  if ((my_row[j] < first_dof_on_proc) || (my_row[j] >= end_dof_on_proc))
382  n_oz[r]++;
383  else
384  n_nz[r]++;
385  }
386  else
387  {
388  n_nz[r] += other.n_nz[r];
389  n_nz[r] = std::min(n_nz[r], n_dofs_on_proc);
390  n_oz[r] += other.n_oz[r];
391  n_oz[r] =std::min(n_oz[r], static_cast<dof_id_type>(n_global_dofs-n_nz[r]));
392  }
393  }
394 
395  // Move nonlocal row information to ourselves; the other thread
396  // won't need it in the map after that.
397  NonlocalGraph::const_iterator it = other.nonlocal_pattern.begin();
398  for (; it != other.nonlocal_pattern.end(); ++it)
399  {
400 #ifndef NDEBUG
401  const dof_id_type dof_id = it->first;
402 
403  processor_id_type dbg_proc_id = 0;
404  while (dof_id >= dof_map.end_dof(dbg_proc_id))
405  dbg_proc_id++;
406  libmesh_assert (dbg_proc_id != this->processor_id());
407 #endif
408 
409  const SparsityPattern::Row & their_row = it->second;
410 
411  // We should have no empty values in a map
412  libmesh_assert (!their_row.empty());
413 
414  NonlocalGraph::iterator my_it = nonlocal_pattern.find(it->first);
415  if (my_it == nonlocal_pattern.end())
416  {
417  // nonlocal_pattern[it->first].swap(their_row);
418  nonlocal_pattern[it->first] = their_row;
419  }
420  else
421  {
422  SparsityPattern::Row & my_row = my_it->second;
423 
424  my_row.insert (my_row.end(),
425  their_row.begin(),
426  their_row.end());
427 
428  // We cannot use SparsityPattern::sort_row() here because it expects
429  // the [begin,middle) [middle,end) to be non-overlapping. This is not
430  // necessarily the case here, so use std::sort()
431  std::sort (my_row.begin(), my_row.end());
432 
433  my_row.erase(std::unique (my_row.begin(), my_row.end()), my_row.end());
434  }
435  }
436 }
437 
438 
439 
441 {
442  parallel_object_only();
444 
445  const dof_id_type n_global_dofs = dof_map.n_dofs();
446  const dof_id_type n_dofs_on_proc = dof_map.n_dofs_on_processor(this->processor_id());
447  const dof_id_type local_first_dof = dof_map.first_dof();
448  const dof_id_type local_end_dof = dof_map.end_dof();
449 
450  // Trade sparsity rows with other processors
451  for (processor_id_type p=1; p != this->n_processors(); ++p)
452  {
453  // Push to processor procup while receiving from procdown
454  processor_id_type procup =
455  cast_int<processor_id_type>((this->processor_id() + p) %
456  this->n_processors());
457  processor_id_type procdown =
458  cast_int<processor_id_type>((this->n_processors() + this->processor_id() - p) %
459  this->n_processors());
460 
461  // Pack the sparsity pattern rows to push to procup
462  std::vector<dof_id_type> pushed_row_ids,
463  pushed_row_ids_to_me;
464  std::vector<std::vector<dof_id_type>> pushed_rows,
465  pushed_rows_to_me;
466 
467  // Move nonlocal row information to a structure to send it from;
468  // we don't need it in the map after that.
469  NonlocalGraph::iterator it = nonlocal_pattern.begin();
470  while (it != nonlocal_pattern.end())
471  {
472  const dof_id_type dof_id = it->first;
473  processor_id_type proc_id = 0;
474  while (dof_id >= dof_map.end_dof(proc_id))
475  proc_id++;
476 
477  libmesh_assert (proc_id != this->processor_id());
478 
479  if (proc_id == procup)
480  {
481  pushed_row_ids.push_back(dof_id);
482 
483  // We can't just do the swap trick here, thanks to the
484  // differing vector allocators?
485  pushed_rows.push_back(std::vector<dof_id_type>());
486  pushed_rows.back().assign
487  (it->second.begin(), it->second.end());
488 
489  nonlocal_pattern.erase(it++);
490  }
491  else
492  ++it;
493  }
494 
495  this->comm().send_receive(procup, pushed_row_ids,
496  procdown, pushed_row_ids_to_me);
497  this->comm().send_receive(procup, pushed_rows,
498  procdown, pushed_rows_to_me);
499  pushed_row_ids.clear();
500  pushed_rows.clear();
501 
502  const std::size_t n_rows = pushed_row_ids_to_me.size();
503  for (std::size_t i=0; i != n_rows; ++i)
504  {
505  const dof_id_type r = pushed_row_ids_to_me[i];
506  const dof_id_type my_r = r - local_first_dof;
507 
508  std::vector<dof_id_type> & their_row = pushed_rows_to_me[i];
509 
511  {
512  SparsityPattern::Row & my_row =
513  sparsity_pattern[my_r];
514 
515  // They wouldn't have sent an empty row
516  libmesh_assert(!their_row.empty());
517 
518  // We can end up with an empty row on a dof that touches our
519  // inactive elements but not our active ones
520  if (my_row.empty())
521  {
522  my_row.assign (their_row.begin(),
523  their_row.end());
524  }
525  else
526  {
527  my_row.insert (my_row.end(),
528  their_row.begin(),
529  their_row.end());
530 
531  // We cannot use SparsityPattern::sort_row() here because it expects
532  // the [begin,middle) [middle,end) to be non-overlapping. This is not
533  // necessarily the case here, so use std::sort()
534  std::sort (my_row.begin(), my_row.end());
535 
536  my_row.erase(std::unique (my_row.begin(), my_row.end()), my_row.end());
537  }
538 
539  // fix the number of on and off-processor nonzeros in this row
540  n_nz[my_r] = n_oz[my_r] = 0;
541 
542  for (std::size_t j=0; j<my_row.size(); j++)
543  if ((my_row[j] < local_first_dof) || (my_row[j] >= local_end_dof))
544  n_oz[my_r]++;
545  else
546  n_nz[my_r]++;
547  }
548  else
549  {
550  for (std::size_t j=0; j<their_row.size(); j++)
551  if ((their_row[j] < local_first_dof) || (their_row[j] >= local_end_dof))
552  n_oz[my_r]++;
553  else
554  n_nz[my_r]++;
555 
556  n_nz[my_r] = std::min(n_nz[my_r], n_dofs_on_proc);
557  n_oz[my_r] = std::min(n_oz[my_r],
558  static_cast<dof_id_type>(n_global_dofs-n_nz[my_r]));
559  }
560  }
561  }
562 
563  // We should have sent everything at this point.
565 }
566 
567 
568 } // namespace SparsityPattern
569 } // namespace libMesh
SparsityPattern::Graph sparsity_pattern
This helper class can be called on multiple threads to compute the sparsity pattern (or graph) of the...
const std::set< GhostingFunctor * > & coupling_functors
The definition of the const_element_iterator struct.
Definition: mesh_base.h:1494
Dummy "splitting object" used to distinguish splitting constructors from copy constructors.
Definition: threads_none.h:63
std::set< GhostingFunctor * >::const_iterator coupling_functors_begin() const
Beginning of range of coupling functors.
Definition: dof_map.h:275
processor_id_type n_processors() const
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
MeshBase & mesh
uint8_t processor_id_type
Definition: id_types.h:99
This proxy class acts like a container of indices from a single coupling row.
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end...
The StoredRange class defines a contiguous, divisible set of objects.
Definition: stored_range.h:52
The libMesh namespace provides an interface to certain functionality in the library.
dof_id_type dof_id
Definition: xdr_io.C:48
const_iterator begin() const
Beginning of the range.
Definition: stored_range.h:261
dof_id_type n_dofs() const
Definition: dof_map.h:510
void handle_vi_vj(const Elem *partner, const std::vector< dof_id_type > &element_dofs_i, unsigned int vj)
Used to iterate over non-NULL entries in a container.
This is the MeshBase class.
Definition: mesh_base.h:68
dof_id_type n_dofs_on_processor(const processor_id_type proc) const
Definition: dof_map.h:526
libmesh_assert(j)
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
std::vector< dof_id_type, Threads::scalable_allocator< dof_id_type > > Row
dof_id_type first_dof(const processor_id_type proc) const
Definition: dof_map.h:535
dof_id_type end_dof(const processor_id_type proc) const
Definition: dof_map.h:577
void send_receive(const unsigned int dest_processor_id, const T1 &send, const unsigned int source_processor_id, T2 &recv, const MessageTag &send_tag=no_tag, const MessageTag &recv_tag=any_tag) const
Send data send to one processor while simultaneously receiving other data recv from a (potentially di...
void join(const Build &other)
const CouplingMatrix * dof_coupling
static const processor_id_type invalid_processor_id
An invalid processor_id to distinguish DoFs that have not been assigned to a processor.
Definition: dof_object.h:335
std::unordered_map< const Elem *, const CouplingMatrix * > map_type
What elements do we care about and what variables do we care about on each element?
void _dummy_function(void)
Dummy function that does nothing but can be used to prohibit compiler optimization in some situations...
std::vector< dof_id_type > n_nz
static void merge_ghost_functor_outputs(GhostingFunctor::map_type &elements_to_ghost, std::set< CouplingMatrix * > &temporary_coupling_matrices, const std::set< GhostingFunctor * >::iterator &gf_begin, const std::set< GhostingFunctor * >::iterator &gf_end, const MeshBase::const_element_iterator &elems_begin, const MeshBase::const_element_iterator &elems_end, processor_id_type p)
Definition: dof_map.C:1408
const_iterator end() const
End of the range.
Definition: stored_range.h:266
This class forms the base class for all other classes that are expected to be implemented in parallel...
std::vector< object_type >::const_iterator const_iterator
Allows an StoredRange to behave like an STL container.
Definition: stored_range.h:58
SparsityPattern::NonlocalGraph nonlocal_pattern
unsigned int n_variables() const
Definition: dof_map.h:477
bool verify(const T &r, const Communicator &comm=Communicator_World)
const Parallel::Communicator & comm() const
void find_connected_dofs(std::vector< dof_id_type > &elem_dofs) const
Finds all the DOFS associated with the element DOFs elem_dofs.
Definition: dof_map.C:2561
void operator()(const ConstElemRange &range)
std::set< GhostingFunctor * >::const_iterator coupling_functors_end() const
End of range of coupling functors.
Definition: dof_map.h:281
unsigned int size() const
long double min(long double a, double b)
Build(const MeshBase &mesh_in, const DofMap &dof_map_in, const CouplingMatrix *dof_coupling_in, std::set< GhostingFunctor * > coupling_functors_in, const bool implicit_neighbor_dofs_in, const bool need_full_sparsity_pattern_in)
uint8_t dof_id_type
Definition: id_types.h:64
processor_id_type processor_id() const
This class defines a coupling matrix.
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917
static void sort_row(const BidirectionalIterator begin, BidirectionalIterator middle, const BidirectionalIterator end)
Splices the two sorted ranges [begin,middle) and [middle,end) into one sorted range [begin...
std::vector< dof_id_type > n_oz