libMesh
boundary_info.C
Go to the documentation of this file.
1 #include <libmesh/mesh.h>
2 #include <libmesh/mesh_generation.h>
3 #include <libmesh/boundary_info.h>
4 #include <libmesh/elem.h>
5 #include <libmesh/face_quad4_shell.h>
6 #include <libmesh/equation_systems.h>
7 #include <libmesh/zero_function.h>
8 #include <libmesh/dirichlet_boundaries.h>
9 #include <libmesh/dof_map.h>
10 #include <libmesh/parallel.h>
11 #include <libmesh/mesh_refinement.h>
12 
13 #include "test_comm.h"
14 #include "libmesh_cppunit.h"
15 
16 #include <regex>
17 
18 using namespace libMesh;
19 
20 class BoundaryInfoTest : public CppUnit::TestCase {
24 public:
25  LIBMESH_CPPUNIT_TEST_SUITE( BoundaryInfoTest );
26 
27  CPPUNIT_TEST( testNameCopying );
28 
29 #if LIBMESH_DIM > 1
30  CPPUNIT_TEST( testMesh );
31  CPPUNIT_TEST( testRenumber );
32 # ifdef LIBMESH_ENABLE_AMR
33 # ifdef LIBMESH_ENABLE_EXCEPTIONS
34  CPPUNIT_TEST( testBoundaryOnChildrenErrors );
35 # endif
36  CPPUNIT_TEST( testBoundaryOnChildrenElementsRefineCoarsen );
37  CPPUNIT_TEST( testBoundaryOnChildrenBoundaryIDs );
38  CPPUNIT_TEST( testBoundaryOnChildrenBoundarySides );
39 # endif
40 # ifdef LIBMESH_ENABLE_DIRICHLET
41  CPPUNIT_TEST( testShellFaceConstraints );
42 # endif
43 #endif
44 #if LIBMESH_DIM > 2
45  CPPUNIT_TEST( testEdgeBoundaryConditions );
46 #endif
47 
48  CPPUNIT_TEST_SUITE_END();
49 
50 protected:
51 
52 public:
53  void setUp()
54  {
55  }
56 
57  void tearDown()
58  {
59  }
60 
61  void testMesh()
62  {
63  LOG_UNIT_TEST;
64 
66 
68  2, 2,
69  0., 1.,
70  0., 1.,
71  QUAD4);
72 
74 
75  // Side lists should be cleared and refilled by each call
76 #ifdef LIBMESH_ENABLE_DEPRECATED
77  std::vector<dof_id_type> element_id_list;
78  std::vector<unsigned short int> side_list;
79  std::vector<boundary_id_type> bc_id_list;
80 #endif
81 
82  // build_square adds boundary_ids 0,1,2,3 for the bottom, right,
83  // top, and left sides, respectively.
84 
85  // On a ReplicatedMesh, we should see all 4 ids on each processor
86  if (mesh.is_serial())
87  CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(4), bi.n_boundary_ids());
88 
89  // On any mesh, we should see each id on *some* processor
90  {
91  const std::set<boundary_id_type> & bc_ids = bi.get_boundary_ids();
92  for (boundary_id_type i = 0 ; i != 4; ++i)
93  {
94  bool has_bcid = bc_ids.count(i);
95  mesh.comm().max(has_bcid);
96  CPPUNIT_ASSERT(has_bcid);
97  }
98  }
99 
100  // Build the side list
101 #ifdef LIBMESH_ENABLE_DEPRECATED
102  bi.build_side_list (element_id_list, side_list, bc_id_list);
103 #endif
104 
105  // Test that the new vector-of-tuples API works equivalently.
106  auto bc_triples = bi.build_side_list();
107 
108  // Check that there are exactly 8 sides in the BoundaryInfo for a
109  // replicated mesh
110  if (mesh.is_serial())
111  {
112 #ifdef LIBMESH_ENABLE_DEPRECATED
113  CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(8), element_id_list.size());
114 #endif
115  CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(8), bc_triples.size());
116  }
117 
118  // Let's test that they are preserved (in a relative sense) when
119  // we clone a mesh.
120  std::unique_ptr<MeshBase> mesh_clone = mesh.clone();
121  CPPUNIT_ASSERT(mesh_clone->get_boundary_info() ==
123 
124  // Let's test that we can remove them successfully.
125  bi.remove_id(0);
126 
127  CPPUNIT_ASSERT(mesh_clone->get_boundary_info() !=
129 
130  // Check that there are now only 3 boundary ids total on the Mesh.
131  if (mesh.is_serial())
132  CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(3), bi.n_boundary_ids());
133 
134  {
135  const std::set<boundary_id_type> & bc_ids = bi.get_boundary_ids();
136  CPPUNIT_ASSERT(!bc_ids.count(0));
137  for (boundary_id_type i = 1 ; i != 4; ++i)
138  {
139  bool has_bcid = bc_ids.count(i);
140  mesh.comm().max(has_bcid);
141  CPPUNIT_ASSERT(has_bcid);
142  }
143  }
144 
145  // Build the side list again
146 #ifdef LIBMESH_ENABLE_DEPRECATED
147  bi.build_side_list (element_id_list, side_list, bc_id_list);
148 #endif
149  bc_triples = bi.build_side_list();
150 
151  // Check that there are now exactly 6 sides left in the
152  // BoundaryInfo on a replicated mesh
153  if (mesh.is_serial())
154  {
155 #ifdef LIBMESH_ENABLE_DEPRECATED
156  CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(6), element_id_list.size());
157 #endif
158  CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(6), bc_triples.size());
159  }
160 
161  // Check that the removed ID is really removed
162 #ifdef LIBMESH_ENABLE_DEPRECATED
163  CPPUNIT_ASSERT(std::find(bc_id_list.begin(), bc_id_list.end(), 0) == bc_id_list.end());
164 #endif
165  typedef std::tuple<dof_id_type, unsigned short int, boundary_id_type> Tuple;
166  CPPUNIT_ASSERT(std::find_if(bc_triples.begin(), bc_triples.end(),
167  [](const Tuple & t)->bool { return std::get<2>(t) == 0; }) == bc_triples.end());
168 
169  // Remove the same id again, make sure nothing changes.
170  bi.remove_id(0);
171  if (mesh.is_serial())
172  CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(3), bi.n_boundary_ids());
173 
174  // Remove the remaining IDs, verify that we have no sides left and
175  // that we can safely reuse the same vectors in the
176  // build_side_list() call.
177  bi.remove_id(1);
178  bi.remove_id(2);
179  bi.remove_id(3);
180 #ifdef LIBMESH_ENABLE_DEPRECATED
181  bi.build_side_list (element_id_list, side_list, bc_id_list);
182 #endif
183  bc_triples = bi.build_side_list();
184 
185  CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(0), bi.n_boundary_ids());
186 #ifdef LIBMESH_ENABLE_DEPRECATED
187  CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(0), element_id_list.size());
188 #endif
189  CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(0), bc_triples.size());
190  }
191 
192 
194  {
195  LOG_UNIT_TEST;
196 
198 
200  2, 2,
201  0., 1.,
202  0., 1.,
203  QUAD4);
204 
206 
207  // Side lists should be cleared and refilled by each call
208 #ifdef LIBMESH_ENABLE_DEPRECATED
209  std::vector<dof_id_type> element_id_list;
210  std::vector<unsigned short int> side_list;
211  std::vector<boundary_id_type> bc_id_list;
212 #endif
213 
214  // build_square adds boundary_ids 0,1,2,3 for the bottom, right,
215  // top, and left sides, respectively. Let's remap those, not 1-1.
216  bi.renumber_id(0, 4);
217  bi.renumber_id(1, 5);
218  bi.renumber_id(2, 6);
219  bi.renumber_id(3, 6);
220 
221  const std::map<boundary_id_type, std::string> expected_names =
222  {{4,"bottom"}, {5,"right"}, {6,"left"}};
223 
224  // On a ReplicatedMesh, we should see ids 4,5,6 on each processor
225  if (mesh.is_serial())
226  CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(3), bi.n_boundary_ids());
227 
228  // On any mesh, we should see each new id on *some* processor, and
229  // shouldn't see old ids on *any* processor
230  {
231  const std::set<boundary_id_type> & bc_ids = bi.get_boundary_ids();
232  for (boundary_id_type i = 0 ; i != 4; ++i)
233  {
234  bool has_bcid = bc_ids.count(i);
235  mesh.comm().max(has_bcid);
236  CPPUNIT_ASSERT(!has_bcid);
237  }
238  for (boundary_id_type i = 4 ; i != 7; ++i)
239  {
240  bool has_bcid = bc_ids.count(i);
241 
242  bool bad_name = false;
243  if (has_bcid)
244  {
245  const std::string & current_name = bi.sideset_name(i);
246 
247  bad_name = (current_name != libmesh_map_find(expected_names, i));
248  }
249 
250  // At least one proc should have each of these BCs
251  mesh.comm().max(has_bcid);
252  CPPUNIT_ASSERT(has_bcid);
253 
254  // No proc should have the wrong name for a BC it has
255  mesh.comm().max(bad_name);
256  CPPUNIT_ASSERT(!bad_name);
257  }
258  }
259 
260  // Check that there are still exactly 8 sides in the BoundaryInfo
261  // for a replicated mesh
262  auto bc_triples = bi.build_side_list();
263 
264  if (mesh.is_serial())
265  {
266  CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(8), bc_triples.size());
267  }
268 
269  // Remove the new IDs, verify that we have no sides left
270  bi.remove_id(4);
271  bi.remove_id(5);
272  bi.remove_id(6);
273  bc_triples = bi.build_side_list();
274 
275  CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(0), bi.n_boundary_ids());
276  CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(0), bc_triples.size());
277  }
278 
279 
281  {
282  LOG_UNIT_TEST;
283 
284  const unsigned int n_elem = 5;
285  const std::string mesh_filename = "cube_mesh.xda";
286 
287  {
290  n_elem, n_elem, n_elem,
291  0., 1.,
292  0., 1.,
293  0., 1.,
294  HEX8);
295 
297 
298  // build_cube does not add any edge boundary IDs
299  CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(0), bi.n_edge_conds());
300 
301  // Let's now add some edge boundary IDs.
302  // We loop over all elements (not just local elements) so that
303  // all processors know about the boundary IDs
304  const boundary_id_type BOUNDARY_ID_MAX_X = 2;
305  const boundary_id_type BOUNDARY_ID_MIN_Y = 1;
306  const boundary_id_type EDGE_BOUNDARY_ID = 20;
307 
308  for (const auto & elem : mesh.element_ptr_range())
309  {
310  unsigned short side_max_x = 0, side_min_y = 0;
311  bool found_side_max_x = false, found_side_min_y = false;
312 
313  for (unsigned short side=0; side<elem->n_sides(); side++)
314  {
315  if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MAX_X))
316  {
317  side_max_x = side;
318  found_side_max_x = true;
319  }
320 
321  if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MIN_Y))
322  {
323  side_min_y = side;
324  found_side_min_y = true;
325  }
326  }
327 
328  // If elem has sides on boundaries
329  // BOUNDARY_ID_MAX_X and BOUNDARY_ID_MIN_Y
330  // then let's set an edge boundary condition
331  if (found_side_max_x && found_side_min_y)
332  for (unsigned short e=0; e<elem->n_edges(); e++)
333  if (elem->is_edge_on_side(e, side_max_x) &&
334  elem->is_edge_on_side(e, side_min_y))
335  bi.add_edge(elem, e, EDGE_BOUNDARY_ID);
336  }
337 
338  // Check that we have the expected number of edge boundary IDs after
339  // updating bi
340  CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(n_elem), bi.n_edge_conds());
341 
342  // Let's test that edge BCIDs are preserved (in a relative
343  // sense) when we clone a mesh.
344  std::unique_ptr<MeshBase> mesh_clone = mesh.clone();
345  CPPUNIT_ASSERT(mesh_clone->get_boundary_info() ==
347 
348  mesh.write(mesh_filename);
349  }
350 
351  // Make sure all processors are done writing before we try to
352  // start reading
354 
356  mesh.read(mesh_filename);
357 
358  // Check that writing and reading preserves the edge boundary IDs
360  CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(n_elem), bi.n_edge_conds());
361  }
362 
364  {
365  LOG_UNIT_TEST;
366 
369  8,
370  0., 1.,
371  EDGE2);
372 
373  Mesh mesh2(mesh);
374 
376  bi.sideset_name(0) = "zero";
377  bi.sideset_name(1) = "one";
378  bi.sideset_name(2) = "two";
379  bi.sideset_name(3) = "three";
380  bi.nodeset_name(0) = "ZERO";
381  bi.nodeset_name(1) = "ONE";
382 
383  BoundaryInfo bi2 {bi};
384  CPPUNIT_ASSERT_EQUAL(bi2.get_sideset_name(0), std::string("zero"));
385  CPPUNIT_ASSERT_EQUAL(bi2.get_sideset_name(1), std::string("one"));
386  CPPUNIT_ASSERT_EQUAL(bi2.get_sideset_name(2), std::string("two"));
387  CPPUNIT_ASSERT_EQUAL(bi2.get_sideset_name(3), std::string("three"));
388  CPPUNIT_ASSERT_EQUAL(bi2.get_nodeset_name(0), std::string("ZERO"));
389  CPPUNIT_ASSERT_EQUAL(bi2.get_nodeset_name(1), std::string("ONE"));
390 
391  BoundaryInfo & bi3 = mesh2.get_boundary_info();
392  bi3 = bi;
393  CPPUNIT_ASSERT_EQUAL(bi3.get_sideset_name(0), std::string("zero"));
394  CPPUNIT_ASSERT_EQUAL(bi3.get_sideset_name(1), std::string("one"));
395  CPPUNIT_ASSERT_EQUAL(bi3.get_sideset_name(2), std::string("two"));
396  CPPUNIT_ASSERT_EQUAL(bi3.get_sideset_name(3), std::string("three"));
397  CPPUNIT_ASSERT_EQUAL(bi3.get_nodeset_name(0), std::string("ZERO"));
398  CPPUNIT_ASSERT_EQUAL(bi3.get_nodeset_name(1), std::string("ONE"));
399  }
400 
401 #ifdef LIBMESH_ENABLE_DIRICHLET
403  {
404  LOG_UNIT_TEST;
405 
406  // Make a simple two element mesh that we can use to test constraints
408 
409  // (0,1) (1,1)
410  // x---------------x
411  // | |
412  // | |
413  // | |
414  // | |
415  // | |
416  // x---------------x
417  // (0,0) (1,0)
418  // | |
419  // | |
420  // | |
421  // | |
422  // x---------------x
423  // (0,-1) (1,-1)
424 
425  mesh.add_point( Point(0.0,-1.0), 4 );
426  mesh.add_point( Point(1.0,-1.0), 5 );
427  mesh.add_point( Point(1.0, 0.0), 1 );
428  mesh.add_point( Point(1.0, 1.0), 2 );
429  mesh.add_point( Point(0.0, 1.0), 3 );
430  mesh.add_point( Point(0.0, 0.0), 0 );
431 
432  Elem * elem_top = mesh.add_elem(Elem::build(QUADSHELL4));
433  elem_top->set_node(0, mesh.node_ptr(0));
434  elem_top->set_node(1, mesh.node_ptr(1));
435  elem_top->set_node(2, mesh.node_ptr(2));
436  elem_top->set_node(3, mesh.node_ptr(3));
437 
438  Elem * elem_bottom = mesh.add_elem(Elem::build(QUADSHELL4));
439  elem_bottom->set_node(0, mesh.node_ptr(4));
440  elem_bottom->set_node(1, mesh.node_ptr(5));
441  elem_bottom->set_node(2, mesh.node_ptr(1));
442  elem_bottom->set_node(3, mesh.node_ptr(0));
443 
445  bi.add_shellface(elem_top, 0, 10);
446  bi.add_shellface(elem_bottom, 1, 20);
447 
448  mesh.allow_renumbering(true);
450 
451  // Let's test that shellface BCIDs are preserved (in a relative
452  // sense) when we clone a mesh.
453  std::unique_ptr<MeshBase> mesh_clone = mesh.clone();
454  CPPUNIT_ASSERT(mesh_clone->get_boundary_info() ==
456 
457  CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(2), bi.n_shellface_conds());
458 
459  EquationSystems es(mesh);
460  System & system = es.add_system<System> ("SimpleSystem");
461  system.add_variable("u", FIRST);
462 
463  // Add a Dirichlet constraint to check that we impose constraints
464  // correctly on shell faces.
465  std::vector<unsigned int> variables;
466  variables.push_back(0);
467  std::set<boundary_id_type> shellface_ids;
468  shellface_ids.insert(20);
469  ZeroFunction<> zf;
470  DirichletBoundary dirichlet_bc(shellface_ids,
471  variables,
472  &zf);
473  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
474  es.init();
475 
476  // Find elem_bottom again if we have it (it may have been deleted
477  // in a DistributedMesh or renumbered in theory)
478  elem_bottom = nullptr;
479  for (unsigned int e = 0; e != mesh.max_elem_id(); ++e)
480  {
481  Elem *elem = mesh.query_elem_ptr(e);
482  if (elem && elem->point(3) == Point(0,0))
483  elem_bottom = elem;
484  }
485 
486  // We expect to have a dof constraint on all four dofs of
487  // elem_bottom
488  CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(4), static_cast<std::size_t>(system.n_constrained_dofs()));
489 
490  // But we may only know the details of that
491  // constraint on the processor which owns elem_bottom.
492  if (elem_bottom &&
493  elem_bottom->processor_id() == mesh.processor_id())
494  {
495  std::vector<dof_id_type> dof_indices;
496  system.get_dof_map().dof_indices(elem_bottom, dof_indices);
497  CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(4), dof_indices.size());
498 
499  for(unsigned int i=0; i<dof_indices.size(); i++)
500  {
501  dof_id_type dof_id = dof_indices[i];
502  CPPUNIT_ASSERT( system.get_dof_map().is_constrained_dof(dof_id) );
503  }
504  }
505  }
506 #endif // LIBMESH_ENABLE_DIRICHLET
507 
508 #if LIBMESH_ENABLE_AMR
509 # if LIBMESH_ENABLE_EXCEPTIONS
511  {
512  LOG_UNIT_TEST;
513 
514  // We create one cell only. The default boundaries of the cell are below.
515  // ___2___
516  // 3 | | 1
517  // |_____|
518  // 0
519 
520  auto mesh = std::make_unique<Mesh>(*TestCommWorld);
522  1, 1,
523  0., 1.,
524  0., 1.,
525  QUAD4);
526 
528 
529  // We only have one element, but for easy access we use the iterator
530  for (auto & elem : mesh->active_element_ptr_range())
531  elem->set_refinement_flag(Elem::REFINE);
533 
536 
537  // Now we try to add boundary id 3 to a child on side 3. This should
538  // result in a "not implemented" error message
539  bool threw_desired_exception = false;
540  try {
541  for (auto & elem : mesh->active_element_ptr_range())
542  {
543  const Point c = elem->vertex_average();
544  if (c(0) < 0.5 && c(1) > 0.5)
545  bi.add_side(elem, 3, 3);
546  }
547  }
548  catch (libMesh::NotImplemented & e) {
549  std::regex msg_regex("Trying to add boundary ID 3 which already exists on the ancestors");
550  CPPUNIT_ASSERT(std::regex_search(e.what(), msg_regex));
551  threw_desired_exception = true;
552  }
553  // If we have more than 4 processors, or a poor partitioner, we
554  // might not get an exception on every processor
555  mesh->comm().max(threw_desired_exception);
556 
557  CPPUNIT_ASSERT(threw_desired_exception);
558 
559  threw_desired_exception = false;
560  try {
561  for (auto & elem : mesh->active_element_ptr_range())
562  {
563  const Point c = elem->vertex_average();
564  if (c(0) < 0.5 && c(1) > 0.5)
565  bi.add_side(elem, 3, {3,4});
566  }
567  }
568  catch (libMesh::NotImplemented & e) {
569  std::regex msg_regex("Trying to add boundary ID 3 which already exists on the ancestors");
570  CPPUNIT_ASSERT(std::regex_search(e.what(), msg_regex));
571  threw_desired_exception = true;
572  }
573 
574  // If we have more than 4 processors, or a poor partitioner, we
575  // might not get an exception on every processor
576  mesh->comm().max(threw_desired_exception);
577 
578  CPPUNIT_ASSERT(threw_desired_exception);
579 
580  // We tested the side addition errors, now we move to the removal parts.
581  // We will attempt the removal of boundary 3 through the child
582  threw_desired_exception = false;
584  try {
585  for (auto & elem : mesh->active_element_ptr_range())
586  {
587  const Point c = elem->vertex_average();
588  if (c(0) < 0.5 && c(1) > 0.5)
589  bi.remove_side(elem, 3, 3);
590  }
591  }
592  catch (libMesh::NotImplemented & e) {
593  std::regex msg_regex("We cannot delete boundary ID 3 using a child because it is inherited from an ancestor");
594  CPPUNIT_ASSERT(std::regex_search(e.what(), msg_regex));
595  threw_desired_exception = true;
596  }
597 
598  // If we have more than 4 processors, or a poor partitioner, we
599  // might not get an exception on every processor
600  mesh->comm().max(threw_desired_exception);
601 
602  CPPUNIT_ASSERT(threw_desired_exception);
603  }
604 # endif // LIBMESH_ENABLE_EXCEPTIONS
605 
607  {
608  LOG_UNIT_TEST;
609 
610  // Set subdomain ids for specific elements, we will refine/coarsen
611  // the cell on subdomain 1
612  // _____________
613  // | 1 | 2 |
614  // |_____|_____|
615 
616  auto mesh = std::make_unique<Mesh>(*TestCommWorld);
618  2, 1,
619  0., 2.,
620  0., 1.,
621  QUAD4);
622 
624 
625  for (auto & elem : mesh->active_element_ptr_range())
626  {
627  const Point c = elem->vertex_average();
628  if (c(0) < 1)
629  {
630  elem->subdomain_id() = 1;
631  elem->set_refinement_flag(Elem::REFINE);
632  }
633  else
634  elem->subdomain_id() = 2;
635  }
637 
638  // Refine the elements once in subdomain 1, and
639  // add the right side subdomain 1 as boundary 5
642 
643  for (auto & elem : mesh->active_element_ptr_range())
644  {
645  const Point c = elem->vertex_average();
646  if (c(0) < 1 && c(0) > 0.5)
647  bi.add_side(elem, 1, 5);
648  }
650 
651  // Check the middle boundary, we expect to have two sides in boundary 5
652  unsigned int count = 0;
653  for (auto & elem : mesh->active_element_ptr_range())
654  if (bi.has_boundary_id(elem, 1, 5))
655  count++;
656 
657  if (mesh->n_active_local_elem())
658  {
659  CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(2), count);
660  CPPUNIT_ASSERT(bi.is_children_on_boundary_side());
661  }
662 
663  // First, we will coarsen the the elements on subdomain 1. This
664  // is to check if the boundary information propagates upward upon
665  // coarsening.
666  for (auto & elem : mesh->active_element_ptr_range())
667  {
668  const Point c = elem->vertex_average();
669  if (c(0) < 1)
670  elem->set_refinement_flag(Elem::COARSEN);
671  }
673 
674  // The coarsened element should have its side on boundary 5
675  // This is boundary info transferred from this child element
678 
679  for (auto & elem : mesh->active_element_ptr_range())
680  {
681  const Point c = elem->vertex_average();
682  if (c(0) < 1)
683  {
684  CPPUNIT_ASSERT(bi.has_boundary_id(elem, 1, 5));
685  // We clean up this boundary ID for the next round of tests
686  bi.remove_side(elem, 1, 5);
687  // we will refine this element again
688  elem->set_refinement_flag(Elem::REFINE);
689  }
690  }
691 
694 
695  // This time we remove boundary 5 from one of the children. We expect
696  // the boundary not to propagate to the next level. Furthermore we
697  // expect boundary 5 to be deleted from the parent's boundaries
698  for (auto & elem : mesh->active_element_ptr_range())
699  {
700  const Point c = elem->vertex_average();
701  if (c(0) < 1)
702  elem->set_refinement_flag(Elem::COARSEN);
703  if (c(0) > 0.5 && c(0) < 1 && c(1) < 0.5)
704  bi.add_side(elem, 1, 5);
705  }
707 
709 
711 
712  // The parent element should not have any side associated with boundary 5
713  for (auto & elem : mesh->active_element_ptr_range())
714  {
715  const Point c = elem->vertex_average();
716  if (c(0) < 1)
717  CPPUNIT_ASSERT(!bi.has_boundary_id(elem, 1, 5));
718  }
719  }
720 
722  {
723  LOG_UNIT_TEST;
724 
725  // We create one cell only. The default boundaries of the cell are below.
726  // We will refine the mesh and add a new boundary id to the left side (side 3).
727  // Then will query the available boundary ids on the added side. It should return
728  // both the parent's and the child's boundaries.
729  // ___2___
730  // 3 | | 1
731  // |_____|
732  // 0
733 
734  auto mesh = std::make_unique<Mesh>(*TestCommWorld);
736  1, 1,
737  0., 1.,
738  0., 1.,
739  QUAD4);
740 
742 
743  // We only have one element, but for easy access we use the iterator
744  for (auto & elem : mesh->active_element_ptr_range())
745  elem->set_refinement_flag(Elem::REFINE);
747 
749 
750  // Now we add the extra boundary ID (5) to the element in the top
751  // left corner
752  for (auto & elem : mesh->active_element_ptr_range())
753  {
754  const Point c = elem->vertex_average();
755  if (c(0) < 0.5 && c(1) > 0.5)
756  bi.add_side(elem, 3, 5);
757  }
759 
760  // Okay, now we query the boundary ids on side 3 of the child and check if it has
761  // the right elements
762  for (auto & elem : mesh->active_element_ptr_range())
763  {
764  const Point c = elem->vertex_average();
765  if (c(0) < 0.5 && c(1) > 0.5)
766  {
767  std::vector<boundary_id_type> container;
768  bi.boundary_ids(elem, 3, container);
769 
770  CPPUNIT_ASSERT_EQUAL(static_cast<std::size_t>(2), container.size());
771  CPPUNIT_ASSERT_EQUAL(static_cast<boundary_id_type>(5), container[0]);
772  CPPUNIT_ASSERT_EQUAL(static_cast<boundary_id_type>(3), container[1]);
773  }
774  }
775  }
776 
778  {
779  LOG_UNIT_TEST;
780 
781  // We create one cell only. The default boundaries of the cell are below.
782  // We will refine mesh and see if we can get back the correct sides
783  // for a given boundary id on an internal boundary.
784  // ___2___
785  // 3 | | 1
786  // |_____|
787  // 0
788 
789  auto mesh = std::make_unique<Mesh>(*TestCommWorld);
791  1, 1,
792  0., 1.,
793  0., 1.,
794  QUAD4);
795 
797 
798 
799  // We only have one element, but for easy access we use the iterator
800  for (auto & elem : mesh->active_element_ptr_range())
801  elem->set_refinement_flag(Elem::REFINE);
804 
805  // Now we add the extra boundary ID (5) to two sides of
806  // the element in the bottom left corner. then we refine again
807  for (auto & elem : mesh->active_element_ptr_range())
808  {
809  const Point c = elem->vertex_average();
810  if (c(0) < 0.5 && c(1) < 0.5)
811  {
812  bi.add_side(elem, 1, 5);
813  bi.add_side(elem, 2, 5);
814  elem->set_refinement_flag(Elem::REFINE);
815  }
816  }
819 
820  // Okay, now we add another boundary id (6) to the cell which is in the bottom
821  // right corner of the refined element
822  for (auto & elem : mesh->active_element_ptr_range())
823  {
824  const Point c = elem->vertex_average();
825  if (c(0) < 0.5 && c(0) > 0.25 && c(1) < 0.25)
826  bi.add_side(elem, 1, 6);
827  }
828 
829  // Time to test if we can get back the boundary sides, first we
830  // check if we can get back boundary from the ancestors of (5) on
831  // the cell which only has boundary (6) registered. We also check
832  // if we can get boundary (6) back.
833 
834  for (auto & elem : mesh->active_element_ptr_range())
835  {
836  const Point c = elem->vertex_average();
837  if (c(0) < 0.5 && c(0) > 0.25 && c(1) < 0.25)
838  {
839  const auto side_5 = bi.side_with_boundary_id(elem, 5);
840  const auto side_6 = bi.side_with_boundary_id(elem, 6);
841  CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(1), side_5);
842  CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(1), side_6);
843  }
844  }
845 
846  // Now we go and try to query the sides with boundary id (5) using
847  // the element which is at the top right corner of the bottom
848  // right parent.
849  for (auto & elem : mesh->active_element_ptr_range())
850  {
851  const Point c = elem->vertex_average();
852  if (c(0) < 0.5 && c(0) > 0.25 && c(1) > 0.25 && c(1) < 0.5)
853  {
854  const auto sides = bi.sides_with_boundary_id(elem, 5);
855  CPPUNIT_ASSERT_EQUAL(static_cast<unsigned long>(2), sides.size());
856  CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(1), sides[0]);
857  CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(2), sides[1]);
858  }
859  }
860  }
861 #endif //LIBMESH_ENABLE_AMR
862 };
863 
void remove_id(boundary_id_type id, bool global=false)
Removes all entities (nodes, sides, edges, shellfaces) with boundary id id from their respective cont...
A class to stub for features that should be in libMesh, but haven&#39;t been written yet, to be thrown by "libmesh_not_implemented();".
void allow_children_on_boundary_side(const bool children_on_boundary)
Whether or not to allow directly setting boundary sides on child elements.
This is the EquationSystems class.
virtual void read(const std::string &name, void *mesh_data=nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false)=0
Interfaces for reading/writing a mesh to/from a file.
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:2558
std::string & nodeset_name(boundary_id_type id)
bool has_boundary_id(const Node *const node, const boundary_id_type id) const
ConstFunction that simply returns 0.
Definition: zero_function.h:38
std::vector< unsigned int > sides_with_boundary_id(const Elem *const elem, const boundary_id_type boundary_id) const
void testShellFaceConstraints()
void allow_renumbering(bool allow)
If false is passed in then this mesh will no longer be renumbered when being prepared for use...
Definition: mesh_base.h:1196
std::size_t n_edge_conds() const
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:2164
dof_id_type n_elem(const MeshBase::const_element_iterator &begin, const MeshBase::const_element_iterator &end)
Count up the number of elements of a specific type (as defined by an iterator range).
Definition: mesh_tools.C:969
libMesh::Parallel::Communicator * TestCommWorld
Definition: driver.C:171
bool refine_elements()
Only refines the user-requested elements.
std::size_t n_shellface_conds() const
void prepare_for_use(const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
Prepare a newly ecreated (or read) mesh for use.
Definition: mesh_base.C:759
CPPUNIT_TEST_SUITE_REGISTRATION(BoundaryInfoTest)
void barrier() const
unsigned int side_with_boundary_id(const Elem *const elem, const boundary_id_type boundary_id) const
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
const Parallel::Communicator & comm() const
void build_square(UnstructuredMesh &mesh, const unsigned int nx, const unsigned int ny, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 2D meshes.
void boundary_ids(const Node *node, std::vector< boundary_id_type > &vec_to_fill) const
Fills a user-provided std::vector with the boundary ids associated with Node node.
The libMesh namespace provides an interface to certain functionality in the library.
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
void renumber_id(boundary_id_type old_id, boundary_id_type new_id)
Changes all entities (nodes, sides, edges, shellfaces) with boundary id old_id to instead be labeled ...
virtual std::unique_ptr< MeshBase > clone() const =0
Virtual "copy constructor".
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id)=0
Add a new Node at Point p to the end of the vertex array, with processor_id procid.
std::size_t n_boundary_ids() const
void build_side_list(std::vector< dof_id_type > &element_id_list, std::vector< unsigned short int > &side_list, std::vector< boundary_id_type > &bc_id_list) const
Creates a list of element numbers, sides, and ids for those sides.
bool coarsen_elements()
Only coarsens the user-requested elements.
void testBoundaryOnChildrenBoundaryIDs()
dof_id_type n_active_local_elem() const
Definition: mesh_base.h:565
virtual bool is_serial() const
Definition: mesh_base.h:211
void testEdgeBoundaryConditions()
Implements (adaptive) mesh refinement algorithms for a MeshBase.
int8_t boundary_id_type
Definition: id_types.h:51
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:96
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition: elem.C:444
virtual dof_id_type max_elem_id() const =0
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:57
bool is_constrained_dof(const dof_id_type dof) const
Definition: dof_map.h:2258
void testBoundaryOnChildrenErrors()
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
std::string & sideset_name(boundary_id_type id)
virtual void write(const std::string &name) const =0
const std::set< boundary_id_type > & get_boundary_ids() const
void remove_side(const Elem *elem, const unsigned short int side)
Removes all boundary conditions associated with side side of element elem, if any exist...
virtual const Elem * query_elem_ptr(const dof_id_type i) const =0
void max(const T &r, T &o, Request &req) const
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
Add side side of element number elem with boundary id id to the boundary information data structure...
void add_shellface(const dof_id_type elem, const unsigned short int shellface, const boundary_id_type id)
Add shell face shellface of element number elem with boundary id id to the boundary information data ...
void build_line(UnstructuredMesh &mesh, const unsigned int nx, const Real xmin=0., const Real xmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 1D meshes.
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
bool is_children_on_boundary_side() const
void testBoundaryOnChildrenBoundarySides()
virtual const Node * node_ptr(const dof_id_type i) const =0
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
processor_id_type processor_id() const
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 Point & point(const unsigned int i) const
Definition: elem.h:2453
void build_cube(UnstructuredMesh &mesh, const unsigned int nx=0, const unsigned int ny=0, const unsigned int nz=0, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const Real zmin=0., const Real zmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
Builds a (elements) cube.
void testBoundaryOnChildrenElementsRefineCoarsen()
dof_id_type n_constrained_dofs() const
Definition: system.C:128
uint8_t dof_id_type
Definition: id_types.h:67
void add_edge(const dof_id_type elem, const unsigned short int edge, const boundary_id_type id)
Add edge edge of element number elem with boundary id id to the boundary information data structure...