libMesh
overlapping_coupling_test.C
Go to the documentation of this file.
1 #include <libmesh/equation_systems.h>
2 #include <libmesh/replicated_mesh.h>
3 #include <libmesh/numeric_vector.h>
4 #include <libmesh/dof_map.h>
5 #include <libmesh/elem.h>
6 #include <libmesh/face_quad4.h>
7 #include <libmesh/ghosting_functor.h>
8 #include <libmesh/coupling_matrix.h>
9 #include <libmesh/quadrature_gauss.h>
10 #include <libmesh/linear_implicit_system.h>
11 #include <libmesh/mesh_refinement.h>
12 #include <libmesh/mesh_modification.h>
13 #include <libmesh/partitioner.h>
14 #include <libmesh/sparse_matrix.h>
15 
16 #include "test_comm.h"
17 #include "libmesh_cppunit.h"
18 
19 #include <memory>
20 
21 using namespace libMesh;
22 
23 
24 // This coupling functor is attempting to mimic in a basic way
25 // a particular use case. The idea is that the mesh has overlapping domains
26 // (we assume only two in this test) and that the elements are not
27 // topologically connected between the overlapping domains, but when the user
28 // solves their problem, some dofs on are coupled between the overlapping
29 // elements. So we want to exercise the GhostingFunctor for this scenario,
30 // in particular both algebraic ghosting (dof ghosting) and coupling
31 // ghosting (sparsity pattern). In the use case, the coupling happens at
32 // quadrature points, so we look for the overlap based on the quadrature points.
33 //
34 // Additionally, the use case also has an issue that the coupling functor
35 // cannot be added until after the System is initialized because the overlap
36 // depends on the solution. So we will exercise that case as well.
38 {
39  public:
40 
42  : GhostingFunctor(),
43  _system(system),
44  _mesh(system.get_mesh()),
45  _point_locator(system.get_mesh().sub_point_locator()),
46  _coupling_matrix(nullptr)
47  {
48  // This is a highly specalized coupling functor where we are assuming
49  // only two subdomains that are overlapping so make sure there's only two.
50  CPPUNIT_ASSERT_EQUAL(2, static_cast<int>(_mesh.n_subdomains()));
51 
52  // We're assuming a specific set of subdomain ids
53  std::set<subdomain_id_type> ids;
54  _mesh.subdomain_ids(ids);
55  CPPUNIT_ASSERT( ids.find(1) != ids.end() );
56  CPPUNIT_ASSERT( ids.find(2) != ids.end() );
57 
58  _point_locator->enable_out_of_mesh_mode();
59  }
60 
62 
63  void set_coupling_matrix (std::unique_ptr<CouplingMatrix> & coupling_matrix)
64  { _coupling_matrix = std::move(coupling_matrix); }
65 
66  virtual void operator()
67  (const MeshBase::const_element_iterator & range_begin,
68  const MeshBase::const_element_iterator & range_end,
70  map_type & coupled_elements) override
71  {
72  std::unique_ptr<PointLocatorBase> sub_point_locator = _mesh.sub_point_locator();
73 
74  for( const auto & elem : as_range(range_begin,range_end) )
75  {
76  std::set<subdomain_id_type> allowed_subdomains;
77  unsigned int var = 0;
78  if( elem->subdomain_id() == 1 )
79  {
80  var = _system.variable_number("V");
81  allowed_subdomains.insert(2);
82  }
83  else if( elem->subdomain_id() == 2 )
84  {
85  var = _system.variable_number("U");
86  allowed_subdomains.insert(1);
87  }
88  else
89  CPPUNIT_FAIL("Something bad happpend");
90 
91  unsigned int dim = 2;
92  FEType fe_type = _system.variable_type(var);
93  std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
94  QGauss qrule (dim, fe_type.default_quadrature_order());
95  fe->attach_quadrature_rule (&qrule);
96 
97  const std::vector<libMesh::Point> & qpoints = fe->get_xyz();
98 
99  fe->reinit(elem);
100 
101  for ( const auto & qp : qpoints )
102  {
103  const Elem * overlapping_elem = (*sub_point_locator)( qp, &allowed_subdomains );
104 
105  if( overlapping_elem && overlapping_elem->processor_id() != p )
106  coupled_elements.emplace(overlapping_elem, _coupling_matrix.get());
107  }
108  }
109  }
110 
111 private:
112 
114 
115  const MeshBase & _mesh;
116 
117  std::unique_ptr<PointLocatorBase> _point_locator;
118 
119  std::unique_ptr<CouplingMatrix> _coupling_matrix;
120 };
121 
122 
123 // We use a custom partioner for this test to help ensure we're
124 // testing the cases we want. In particular, we're going to
125 // ensure that elements in subdomain one are on different processors
126 // from elements subdomain two so we know ghosting will be required
127 // for algebraic and coupling ghosting between the two subdomains.
129 {
130 public:
131 
132  OverlappingTestPartitioner () = default;
135  OverlappingTestPartitioner & operator= (const OverlappingTestPartitioner &) = default;
136  OverlappingTestPartitioner & operator= (OverlappingTestPartitioner &&) = default;
137  virtual ~OverlappingTestPartitioner() = default;
138 
142  virtual std::unique_ptr<Partitioner> clone () const override
143  {
144  return std::make_unique<OverlappingTestPartitioner>(*this);
145  }
146 
147 protected:
148 
152  virtual void _do_partition (MeshBase & mesh,
153  const unsigned int n) override
154  {
155  // If we're on one partition, then everyone gets to be on that partition
156  if (n == 1)
157  this->single_partition_range (mesh.active_elements_begin(), mesh.active_elements_end());
158  else
159  {
160  libmesh_assert_greater (n, 0);
161 
162  // If there are more partitions than elements, then just
163  // assign one element to each partition
164  if (mesh.n_active_elem() <= n)
165  {
166  processor_id_type e = 0;
167  for (auto & elem : mesh.active_element_ptr_range() )
168  {
169  elem->processor_id() = e;
170  e++;
171  }
172  }
173  // Otherwise, we'll split up the partitions into two groups
174  // and then assign the elements of each subdomain into each of
175  // the two groups
176  else
177  {
178  unsigned int n_sub_two_elems = std::distance( mesh.active_subdomain_elements_begin(2),
179  mesh.active_subdomain_elements_end(2) );
180 
181  unsigned int n_parts_sub_two = 0;
182  if( n_sub_two_elems < n/2 )
183  n_parts_sub_two = n_sub_two_elems;
184  else
185  n_parts_sub_two = n/2;
186 
187  const unsigned int n_parts_sub_one = n - n_parts_sub_two;
188 
189  const dof_id_type sub_two_blk_size = cast_int<dof_id_type>
190  (std::distance( mesh.active_subdomain_elements_begin(2),
191  mesh.active_subdomain_elements_end(2) )/n_parts_sub_two );
192 
193  this->assign_proc_id_subdomain( mesh, 2, sub_two_blk_size, n_parts_sub_two, 0 );
194 
195 
196  const dof_id_type sub_one_blk_size = cast_int<dof_id_type>
197  (std::distance( mesh.active_subdomain_elements_begin(1),
198  mesh.active_subdomain_elements_end(1) )/n_parts_sub_one );
199 
200  this->assign_proc_id_subdomain( mesh, 1, sub_one_blk_size, n_parts_sub_one, n_parts_sub_two );
201  }
202  }
203  }
204 
206  const subdomain_id_type sid,
207  const dof_id_type blksize,
208  const unsigned int n_parts,
209  const processor_id_type offset )
210  {
211  dof_id_type e = 0;
212  for (auto & elem : mesh.active_subdomain_elements_ptr_range(sid))
213  {
214  if ((e/blksize) < n_parts)
215  elem->processor_id() = offset + cast_int<processor_id_type>(e/blksize);
216  else
217  elem->processor_id() = offset;
218 
219  e++;
220  }
221  }
222 
223 };
224 
225 
226 // Common functionality for all the different tests we'll be running
227 // We're creating two elements that will live in subdomain 1 and one
228 // element that is lives in subdomain 2 that overlaps both elements in
229 // subdomain 1, but is not topologically connected.
231 {
232 protected:
233 
234  std::unique_ptr<MeshBase> _mesh;
235 
236  std::unique_ptr<EquationSystems> _es;
237 
238  void build_quad_mesh(unsigned int n_refinements = 0)
239  {
240  // We are making assumptions in various places about the presence
241  // of the elements on the current processor so we're restricting to
242  // ReplicatedMesh for now.
243  _mesh = std::make_unique<ReplicatedMesh>(*TestCommWorld);
244 
245  _mesh->set_mesh_dimension(2);
246 
247  _mesh->add_point( Point(0.0,0.0),0 );
248  _mesh->add_point( Point(1.0,0.0),1 );
249  _mesh->add_point( Point(1.0,1.0),2 );
250  _mesh->add_point( Point(0.0,1.0),3 );
251 
252  {
253  Elem * elem = _mesh->add_elem(Elem::build_with_id(QUAD4, 0));
254  elem->subdomain_id() = 1;
255 
256  for (unsigned int n=0; n<4; n++)
257  elem->set_node(n, _mesh->node_ptr(n));
258  }
259 
260  _mesh->add_point( Point(1.0,2.0),4 );
261  _mesh->add_point( Point(0.0,2.0),5 );
262 
263  {
264  Elem * elem = _mesh->add_elem(Elem::build_with_id(QUAD4, 1));
265  elem->subdomain_id() = 1;
266 
267  elem->set_node(0, _mesh->node_ptr(3));
268  elem->set_node(1, _mesh->node_ptr(2));
269  elem->set_node(2, _mesh->node_ptr(4));
270  elem->set_node(3, _mesh->node_ptr(5));
271  }
272 
273  _mesh->add_point( Point(0.0,0.0),6 );
274  _mesh->add_point( Point(1.0,0.0),7 );
275  _mesh->add_point( Point(1.0,2.0),8 );
276  _mesh->add_point( Point(0.0,2.0),9 );
277 
278  {
279  Elem* elem = _mesh->add_elem(Elem::build_with_id(QUAD4, 2));
280  elem->subdomain_id() = 2;
281 
282  elem->set_node(0, _mesh->node_ptr(6));
283  elem->set_node(1, _mesh->node_ptr(7));
284  elem->set_node(2, _mesh->node_ptr(8));
285  elem->set_node(3, _mesh->node_ptr(9));
286  }
287 
288  _mesh->partitioner() = std::make_unique<OverlappingTestPartitioner>();
289 
290  _mesh->prepare_for_use();
291 
292 #ifdef LIBMESH_ENABLE_AMR
293  if (n_refinements > 0)
294  {
295  MeshRefinement refine(*_mesh);
296  refine.uniformly_refine(n_refinements);
297  }
298 #endif // LIBMESH_ENABLE_AMR
299  }
300 
301  void init()
302  {
303  _es = std::make_unique<EquationSystems>(*_mesh);
304  LinearImplicitSystem & sys = _es->add_system<LinearImplicitSystem> ("SimpleSystem");
305 
306  std::set<subdomain_id_type> sub_one;
307  sub_one.insert(1);
308 
309  std::set<subdomain_id_type> sub_two;
310  sub_two.insert(2);
311 
312  sys.add_variable("U", FIRST, LAGRANGE, &sub_two);
313  sys.add_variable("L", FIRST, LAGRANGE, &sub_two);
314 
315  sys.add_variable("V", FIRST, LAGRANGE, &sub_one);
316  sys.add_variable("p", FIRST, LAGRANGE, &sub_one);
317 
318  _es->init();
319  }
320 
321  void clear()
322  {
323  _es.reset();
324  _mesh.reset();
325  }
326 
327  void setup_coupling_matrix( std::unique_ptr<CouplingMatrix> & coupling )
328  {
329  System & system = _es->get_system("SimpleSystem");
330 
331  coupling = std::make_unique<CouplingMatrix>(system.n_vars());
332 
333  const unsigned int u_var = system.variable_number("U");
334  const unsigned int l_var = system.variable_number("L");
335  const unsigned int v_var = system.variable_number("V");
336  const unsigned int p_var = system.variable_number("p");
337 
338  // Only adding the overlapping couplings since the primary ones should
339  // be there by default.
340  (*coupling)(u_var,v_var) = true;
341  (*coupling)(l_var,v_var) = true;
342  (*coupling)(l_var,p_var) = true;
343  (*coupling)(v_var,u_var) = true;
344  (*coupling)(v_var,l_var) = true;
345  }
346 
347 };
348 
349 // Mainly just to sanity check the mesh construction and
350 // the assumptions made in the OverlappingCouplingFunctor
351 // as well as the custom Partitioner.
352 class OverlappingFunctorTest : public CppUnit::TestCase,
353  public OverlappingTestBase
354 {
355 public:
356 
357  LIBMESH_CPPUNIT_TEST_SUITE( OverlappingFunctorTest );
358 
359 #ifdef LIBMESH_HAVE_SOLVER
360  CPPUNIT_TEST( checkCouplingFunctorQuad );
361  CPPUNIT_TEST( checkCouplingFunctorQuadUnifRef );
362  CPPUNIT_TEST( checkCouplingFunctorTri );
363  CPPUNIT_TEST( checkCouplingFunctorTriUnifRef );
364  CPPUNIT_TEST( checkOverlappingPartitioner );
365  CPPUNIT_TEST( checkOverlappingPartitionerUnifRef );
366 #endif
367 
368  CPPUNIT_TEST_SUITE_END();
369 
370 public:
371 
372  public:
373  void setUp()
374  {
375  this->build_quad_mesh();
376  this->init();
377  }
378 
379  void tearDown()
380  { this->clear(); }
381 
383  { this->run_coupling_functor_test(0); }
384 
386  { this->run_coupling_functor_test(1); }
387 
389  {
391  _es->reinit();
392  this->run_coupling_functor_test(0);
393  }
394 
396  {
398  _es->reinit();
399  this->run_coupling_functor_test(1);
400  }
401 
403  {
404  this->run_partitioner_test(0);
405  }
406 
408  {
409  this->run_partitioner_test(1);
410  }
411 
412 private:
413 
414  // This is basically to sanity check the coupling functor
415  // with the supplied mesh to make sure all the assumptions
416  // are kosher. This test requires AMR and some kind of a
417  // linear solver, so let's only run it if we have PETSc.
418  void run_coupling_functor_test(unsigned int n_refinements)
419  {
420 #if defined(LIBMESH_ENABLE_AMR) && defined(LIBMESH_HAVE_SOLVER)
421  if( n_refinements > 0 )
422  {
423  MeshRefinement refine(*_mesh);
424  refine.uniformly_refine(n_refinements);
425  _es->reinit();
426  }
427 #endif
428 
429  System & system = _es->get_system("SimpleSystem");
430 
431  OverlappingCouplingFunctor coupling_functor(system);
432 
433  GhostingFunctor::map_type subdomain_one_couplings;
434  GhostingFunctor::map_type subdomain_two_couplings;
435 
436  coupling_functor( _mesh->active_subdomain_elements_begin(1),
437  _mesh->active_subdomain_elements_end(1),
439  subdomain_one_couplings );
440 
441  coupling_functor( _mesh->active_subdomain_elements_begin(2),
442  _mesh->active_subdomain_elements_end(2),
444  subdomain_two_couplings );
445 
446  dof_id_type n_elems_subdomain_two = std::distance( _mesh->active_subdomain_elements_begin(2),
447  _mesh->active_subdomain_elements_end(2) );
448 
449  CPPUNIT_ASSERT_EQUAL(n_elems_subdomain_two, static_cast<dof_id_type>(subdomain_one_couplings.size()));
450  CPPUNIT_ASSERT_EQUAL(static_cast<dof_id_type>(2*n_elems_subdomain_two), static_cast<dof_id_type>(subdomain_two_couplings.size()));
451  }
452 
453  void run_partitioner_test(unsigned int n_refinements)
454  {
455 #ifdef LIBMESH_ENABLE_AMR
456  if( n_refinements > 0 )
457  {
458  MeshRefinement refine(*_mesh);
459  refine.uniformly_refine(n_refinements);
460  _es->reinit();
461  }
462 #endif // LIBMESH_ENABLE_AMR
463 
464  System & system = _es->get_system("SimpleSystem");
465 
466  std::set<processor_id_type> sub_one_proc_ids, sub_two_proc_ids;
467  for (auto & elem : _mesh->active_subdomain_elements_ptr_range(1))
468  sub_one_proc_ids.insert(elem->processor_id());
469 
470  for (auto & elem : _mesh->active_subdomain_elements_ptr_range(2))
471  sub_two_proc_ids.insert(elem->processor_id());
472 
473 
474  // Everyone should be on the same processor if only 1 processor
475  if( system.n_processors() == 1 )
476  {
477  CPPUNIT_ASSERT_EQUAL(1, static_cast<int>(sub_one_proc_ids.size()));
478  CPPUNIT_ASSERT_EQUAL(1, static_cast<int>(sub_two_proc_ids.size()));
479  CPPUNIT_ASSERT( sub_one_proc_ids == sub_two_proc_ids );
480  }
481  // Otherwise these sets should be disjoint
482  else
483  {
484  // Make sure no subdomain one ids in the subdomain two ids
485  for (auto & id : sub_one_proc_ids )
486  CPPUNIT_ASSERT( sub_two_proc_ids.find(id) == sub_two_proc_ids.end() );
487 
488  // Vice-versa
489  for (auto & id : sub_two_proc_ids )
490  CPPUNIT_ASSERT( sub_one_proc_ids.find(id) == sub_one_proc_ids.end() );
491  }
492  }
493 
494 };
495 
496 // In this testing, we're relying on the presence of PETSc
497 #ifdef LIBMESH_HAVE_PETSC
498 
499 
500 // This testing class the algebraic ghosting of the
501 // OverlappingCouplingFunctor.
502 class OverlappingAlgebraicGhostingTest : public CppUnit::TestCase,
503  public OverlappingTestBase
504 {
505 public:
506  LIBMESH_CPPUNIT_TEST_SUITE( OverlappingAlgebraicGhostingTest );
507 
508  CPPUNIT_TEST( testGhostingCouplingMatrix );
509  CPPUNIT_TEST( testGhostingNullCouplingMatrix );
510  CPPUNIT_TEST( testGhostingNullCouplingMatrixUnifRef );
511 
512  CPPUNIT_TEST_SUITE_END();
513 
514 public:
515 
516  void setUp()
517  {}
518 
519  void tearDown()
520  { this->clear(); }
521 
523  {
524  LOG_UNIT_TEST;
525  this->run_ghosting_test(0, true);
526  }
527 
529  {
530  LOG_UNIT_TEST;
531  this->run_ghosting_test(0, false);
532  }
533 
535  {
536  LOG_UNIT_TEST;
537  std::unique_ptr<CouplingMatrix> coupling_matrix;
538 
539  this->run_ghosting_test(2, false);
540  }
541 
542 private:
543 
544  void run_ghosting_test(const unsigned int n_refinements, bool build_coupling_matrix)
545  {
546  this->build_quad_mesh(n_refinements);
547  this->init();
548 
549  std::unique_ptr<CouplingMatrix> coupling_matrix;
550  if (build_coupling_matrix)
551  this->setup_coupling_matrix(coupling_matrix);
552 
553  LinearImplicitSystem & system = _es->get_system<LinearImplicitSystem>("SimpleSystem");
554 
555  // If we don't add this coupling functor and properly recompute the
556  // sparsity pattern, then PETSc will throw a malloc error when we
557  // try to assemble into the global matrix
558  OverlappingCouplingFunctor functor(system);
559  functor.set_coupling_matrix(coupling_matrix);
560 
561  DofMap & dof_map = system.get_dof_map();
562  dof_map.add_algebraic_ghosting_functor(functor);
563  dof_map.reinit_send_list(system.get_mesh());
564 
565  // Update current local solution
567 
568  system.current_local_solution->init(system.n_dofs(), system.n_local_dofs(),
569  dof_map.get_send_list(), false,
571 
572  system.solution->localize(*(system.current_local_solution),dof_map.get_send_list());
573 
574  std::unique_ptr<PointLocatorBase> point_locator = _mesh->sub_point_locator();
575 
576  const unsigned int u_var = system.variable_number("U");
577 
579 
580  FEMContext subdomain_one_context(system);
581  FEMContext subdomain_two_context(system);
582 
583  // The use case on which this test is based, we only add terms to the residual
584  // corresponding to the dofs in the second subdomain, but that have couplings
585  // to dofs in the first subdomain.
586  for (const auto & elem : _mesh->active_local_subdomain_elements_ptr_range(2))
587  {
588  // A little extra unit testing on the range iterator
589  CPPUNIT_ASSERT_EQUAL(2, static_cast<int>(elem->subdomain_id()));
590 
591  subdomain_one_context.get_element_fe(u_var)->get_nothing(); // for this unit test
592  const std::vector<libMesh::Point> & qpoints = subdomain_two_context.get_element_fe(u_var)->get_xyz();
593 
594  // Setup the context for the current element
595  subdomain_two_context.pre_fe_reinit(system,elem);
596  subdomain_two_context.elem_fe_reinit();
597 
598  std::set<subdomain_id_type> allowed_subdomains;
599  allowed_subdomains.insert(1);
600 
601  // Now loop over the quadrature points and find the subdomain-one element that overlaps
602  // with the current subdomain-two element and then initialize the FEMContext for the
603  // subdomain-one element. If the algebraic ghosting has not been done properly,
604  // this will error.
605  for ( const auto & qp : qpoints )
606  {
607  const Elem * overlapping_elem = (*point_locator)( qp, &allowed_subdomains );
608  CPPUNIT_ASSERT(overlapping_elem);
609 
610  // Setup the context for the overlapping element
611  subdomain_one_context.pre_fe_reinit(system,overlapping_elem);
612  subdomain_one_context.elem_fe_reinit();
613  }
614  }
615  }
616 
617 };
618 
619 
620 // This testing class now exercises testing the sparsity
621 // pattern augmented by the OverlappingCouplingFunctor
622 class OverlappingCouplingGhostingTest : public CppUnit::TestCase,
623  public OverlappingTestBase
624 {
625 public:
626  LIBMESH_CPPUNIT_TEST_SUITE( OverlappingCouplingGhostingTest );
627 
628  CPPUNIT_TEST( testSparsityCouplingMatrix );
629  CPPUNIT_TEST( testSparsityNullCouplingMatrix );
630  CPPUNIT_TEST( testSparsityNullCouplingMatrixUnifRef );
631 
632  CPPUNIT_TEST_SUITE_END();
633 
634 public:
635 
636  void setUp()
637  {}
638 
639  void tearDown()
640  { this->clear(); }
641 
643  {
644  LOG_UNIT_TEST;
645  this->run_sparsity_pattern_test(0, true);
646  }
647 
649  {
650  LOG_UNIT_TEST;
651  this->run_sparsity_pattern_test(0, false);
652  }
653 
655  {
656  LOG_UNIT_TEST;
657  this->run_sparsity_pattern_test(1, false);
658  }
659 
660 private:
661 
662  void run_sparsity_pattern_test(const unsigned int n_refinements, bool build_coupling_matrix)
663  {
664  this->build_quad_mesh(n_refinements);
665  this->init();
666 
667  std::unique_ptr<CouplingMatrix> coupling_matrix;
668  if (build_coupling_matrix)
669  this->setup_coupling_matrix(coupling_matrix);
670 
671  LinearImplicitSystem & system = _es->get_system<LinearImplicitSystem>("SimpleSystem");
672 
673  // If we don't add this coupling functor and properly recompute the
674  // sparsity pattern, then PETSc will throw a malloc error when we
675  // try to assemble into the global matrix
676  OverlappingCouplingFunctor coupling_functor(system);
677  coupling_functor.set_coupling_matrix(coupling_matrix);
678 
679  DofMap & dof_map = system.get_dof_map();
680  dof_map.add_coupling_functor(coupling_functor);
681  dof_map.reinit_send_list(system.get_mesh());
682 
683  // Update current local solution
685 
686  system.current_local_solution->init(system.n_dofs(), system.n_local_dofs(),
687  dof_map.get_send_list(), false,
689 
690  system.solution->localize(*(system.current_local_solution),dof_map.get_send_list());
691 
692  // Now that we've added the coupling functor, we need
693  // to recompute the sparsity
694  dof_map.clear_sparsity();
695  dof_map.compute_sparsity(system.get_mesh());
696 
697  // Now that we've recomputed the sparsity pattern, we need
698  // to reinitialize the system matrix.
699  SparseMatrix<Number> & matrix = system.get_system_matrix();
700  libmesh_assert(dof_map.is_attached(matrix));
701  matrix.init();
702 
703  std::unique_ptr<PointLocatorBase> point_locator = _mesh->sub_point_locator();
704 
705  const unsigned int u_var = system.variable_number("U");
706  const unsigned int v_var = system.variable_number("V");
707 
708  DenseMatrix<Number> K12, K21;
709 
710  FEMContext subdomain_one_context(system);
711  subdomain_one_context.get_element_fe(u_var)->get_nothing();
712  subdomain_one_context.get_element_fe(v_var)->get_nothing();
713 
714  FEMContext subdomain_two_context(system);
715  subdomain_two_context.get_element_fe(u_var)->get_xyz();
716  subdomain_two_context.get_element_fe(v_var)->get_nothing();
717 
718  // Add normally coupled parts of the matrix
719  for (const auto & elem : _mesh->active_local_subdomain_elements_ptr_range(1))
720  {
721  subdomain_one_context.pre_fe_reinit(system,elem);
722  subdomain_one_context.elem_fe_reinit();
723 
724  std::vector<dof_id_type> & rows = subdomain_one_context.get_dof_indices();
725 
726  // Fill with ones in case PETSc ignores the zeros at some point
727  std::fill( subdomain_one_context.get_elem_jacobian().get_values().begin(),
728  subdomain_one_context.get_elem_jacobian().get_values().end(),
729  1);
730 
731  // Insert the Jacobian for the dofs for this element
732  matrix.add_matrix( subdomain_one_context.get_elem_jacobian(), rows );
733  }
734 
735  for (const auto & elem : _mesh->active_local_subdomain_elements_ptr_range(2))
736  {
737  // A little extra unit testing on the range iterator
738  CPPUNIT_ASSERT_EQUAL(2, static_cast<int>(elem->subdomain_id()));
739 
740  const std::vector<libMesh::Point> & qpoints = subdomain_two_context.get_element_fe(u_var)->get_xyz();
741 
742  // Setup the context for the current element
743  subdomain_two_context.pre_fe_reinit(system,elem);
744  subdomain_two_context.elem_fe_reinit();
745 
746  // We're only assembling rows for the dofs on subdomain 2 (U,L), so
747  // the current element will have all those dof_indices.
748  std::vector<dof_id_type> & rows = subdomain_two_context.get_dof_indices();
749 
750  std::fill( subdomain_two_context.get_elem_jacobian().get_values().begin(),
751  subdomain_two_context.get_elem_jacobian().get_values().end(),
752  1);
753 
754  // Insert the Jacobian for the normally coupled dofs for this element
755  matrix.add_matrix( subdomain_two_context.get_elem_jacobian(), rows );
756 
757  std::set<subdomain_id_type> allowed_subdomains;
758  allowed_subdomains.insert(1);
759 
760  // Now loop over the quadrature points and find the subdomain-one element that overlaps
761  // with the current subdomain-two element and then add a local element matrix with
762  // the coupling to the global matrix to try and trip any issues with sparsity pattern
763  // construction
764  for ( const auto & qp : qpoints )
765  {
766  const Elem * overlapping_elem = (*point_locator)( qp, &allowed_subdomains );
767  CPPUNIT_ASSERT(overlapping_elem);
768 
769  // Setup the context for the overlapping element
770  subdomain_one_context.pre_fe_reinit(system,overlapping_elem);
771  subdomain_one_context.elem_fe_reinit();
772 
773  // We're only coupling to the "V" variable so only need those dof indices
774  std::vector<dof_id_type> & v_indices = subdomain_one_context.get_dof_indices(v_var);
775  std::vector<dof_id_type> columns(rows);
776  columns.insert( columns.end(), v_indices.begin(), v_indices.end() );
777 
778  // This will also zero the matrix so we can just insert zeros for this test
779  K21.resize( rows.size(), columns.size() );
780 
781  std::fill(K21.get_values().begin(), K21.get_values().end(), 1);
782 
783  // Now adding this local matrix to the global would trip a PETSc
784  // malloc error if the sparsity pattern hasn't been correctly
785  // built to include the overlapping coupling.
786  matrix.add_matrix (K21, rows, columns);
787 
788  // Now add the other part of the overlapping coupling
789  K12.resize(v_indices.size(), rows.size());
790  std::fill(K12.get_values().begin(), K12.get_values().end(), 1);
791  matrix.add_matrix(K12,v_indices,rows);
792  }
793  } // end element loop
794 
795  // We need to make sure to close the matrix for this test. There could still
796  // be PETSc malloc errors tripped here if we didn't allocate the off-processor
797  // part of the sparsity pattern correctly.
798  matrix.close();
799  }
800 
801 };
802 #endif // LIBMESH_HAVE_PETSC
803 
804 
806 
807 #ifdef LIBMESH_HAVE_PETSC
810 #endif
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
void run_coupling_functor_test(unsigned int n_refinements)
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:274
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:2558
virtual dof_id_type n_active_elem() const =0
This abstract base class defines the interface by which library code and user code can report associa...
virtual void pre_fe_reinit(const System &, const Elem *e)
Reinitializes local data vectors/matrices on the current geometric element.
Definition: fem_context.C:1683
void setup_coupling_matrix(std::unique_ptr< CouplingMatrix > &coupling)
The definition of the const_element_iterator struct.
Definition: mesh_base.h:2216
unsigned int dim
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and linear solvers ...
virtual std::unique_ptr< Partitioner > clone() const override
bool is_attached(SparseMatrix< Number > &matrix)
Matrices should not be attached more than once.
Definition: dof_map.C:329
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
const Parallel::Communicator & comm() const
Order default_quadrature_order() const
Definition: fe_type.h:371
std::unique_ptr< CouplingMatrix > _coupling_matrix
std::map< const Elem *, const CouplingMatrix *, CompareDofObjectsByPIDAndThenID > map_type
What elements do we care about and what variables do we care about on each element?
The libMesh namespace provides an interface to certain functionality in the library.
dof_id_type n_local_dofs() const
Definition: system.C:158
void run_partitioner_test(unsigned int n_refinements)
const MeshBase & get_mesh() const
Definition: system.h:2358
Real distance(const Point &p)
void add_coupling_functor(GhostingFunctor &coupling_functor, bool to_mesh=true)
Adds a functor which can specify coupling requirements for creation of sparse matrices.
Definition: dof_map.C:2033
dof_id_type n_dofs() const
Definition: system.C:121
This is the MeshBase class.
Definition: mesh_base.h:75
virtual void elem_fe_reinit(const std::vector< Point > *const pts=nullptr)
Reinitializes interior FE objects on the current geometric element.
Definition: fem_context.C:1477
The Partitioner class provides a uniform interface for partitioning algorithms.
Definition: partitioner.h:51
unsigned int variable_number(std::string_view var) const
Definition: system.C:1609
std::unique_ptr< MeshBase > _mesh
std::unique_ptr< EquationSystems > _es
virtual void _do_partition(MeshBase &mesh, const unsigned int n) override
Partition the MeshBase into n subdomains.
std::unique_ptr< PointLocatorBase > _point_locator
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
void build_quad_mesh(unsigned int n_refinements=0)
uint8_t processor_id_type
void all_tri(MeshBase &mesh)
Subdivides any non-simplex elements in a Mesh to produce simplex (triangular in 2D, tetrahedral in 3D) elements.
void reinit_send_list(MeshBase &mesh)
Clears the _send_list vector and then rebuilds it.
Definition: dof_map.C:1906
processor_id_type n_processors() const
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
Add the full matrix dm to the SparseMatrix.
Implements (adaptive) mesh refinement algorithms for a MeshBase.
CPPUNIT_TEST_SUITE_REGISTRATION(OverlappingFunctorTest)
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:493
void assign_proc_id_subdomain(MeshBase &mesh, const subdomain_id_type sid, const dof_id_type blksize, const unsigned int n_parts, const processor_id_type offset)
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:96
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
Helper function that allows us to treat a homogenous pair as a range.
Definition: simple_range.h:57
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1593
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
libmesh_assert(ctx)
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:62
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:363
unsigned int add_variable(std::string_view var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1357
void run_ghosting_test(const unsigned int n_refinements, bool build_coupling_matrix)
std::vector< T > & get_values()
Definition: dense_matrix.h:382
virtual void close()=0
Calls the SparseMatrix&#39;s internal assembly routines, ensuring that the values are consistent across p...
void set_coupling_matrix(std::unique_ptr< CouplingMatrix > &coupling_matrix)
subdomain_id_type subdomain_id() const
Definition: elem.h:2582
static std::unique_ptr< Elem > build_with_id(const ElemType type, dof_id_type id)
Calls the build() method above with a nullptr parent, and additionally sets the newly-created Elem&#39;s ...
Definition: elem.C:558
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition: dense_matrix.h:895
static std::unique_ptr< NumericVector< T > > build(const Parallel::Communicator &comm, SolverPackage solver_package=libMesh::default_solver_package(), ParallelType parallel_type=AUTOMATIC)
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
void clear_sparsity()
Clears the sparsity pattern.
Definition: dof_map.C:2009
This class implements specific orders of Gauss quadrature.
unsigned int level ElemType type std::set< subdomain_id_type > ss processor_id_type pid unsigned int level std::set< subdomain_id_type > virtual ss SimpleRange< element_iterator > active_subdomain_elements_ptr_range(subdomain_id_type sid)=0
std::unique_ptr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1605
void get_element_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for interior finite element object for variable var for the largest dimension in the mesh...
Definition: fem_context.h:277
void run_sparsity_pattern_test(const unsigned int n_refinements, bool build_coupling_matrix)
unsigned int n_vars() const
Definition: system.h:2430
void compute_sparsity(const MeshBase &)
Computes the sparsity pattern for the matrices corresponding to proc_id and sends that data to Linear...
Definition: dof_map.C:1988
void add_algebraic_ghosting_functor(GhostingFunctor &evaluable_functor, bool to_mesh=true)
Adds a functor which can specify algebraic ghosting requirements for use with distributed vectors...
Definition: dof_map.C:2058
const DofMap & get_dof_map() const
Definition: system.h:2374
processor_id_type processor_id() const
Definition: dof_object.h:905
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
const SparseMatrix< Number > & get_system_matrix() const
const std::vector< dof_id_type > & get_send_list() const
Definition: dof_map.h:526
uint8_t dof_id_type
Definition: id_types.h:67
virtual void init(const numeric_index_type m, const numeric_index_type n, const numeric_index_type m_l, const numeric_index_type n_l, const numeric_index_type nnz=30, const numeric_index_type noz=10, const numeric_index_type blocksize=1)=0
Initialize SparseMatrix with the specified sizes.
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.