libMesh
boundary_mesh.C
Go to the documentation of this file.
1 // Ignore unused parameter warnings coming from cppunit headers
2 #include <libmesh/ignore_warnings.h>
3 #include <cppunit/extensions/HelperMacros.h>
4 #include <cppunit/TestCase.h>
5 #include <libmesh/restore_warnings.h>
6 
7 #include <libmesh/mesh.h>
8 #include <libmesh/mesh_generation.h>
9 #include <libmesh/mesh_refinement.h>
10 #include <libmesh/remote_elem.h>
11 #include <libmesh/replicated_mesh.h>
12 
13 #include "test_comm.h"
14 
15 // THE CPPUNIT_TEST_SUITE_END macro expands to code that involves
16 // std::auto_ptr, which in turn produces -Wdeprecated-declarations
17 // warnings. These can be ignored in GCC as long as we wrap the
18 // offending code in appropriate pragmas. We can't get away with a
19 // single ignore_warnings.h inclusion at the beginning of this file,
20 // since the libmesh headers pull in a restore_warnings.h at some
21 // point. We also don't bother restoring warnings at the end of this
22 // file since it's not a header.
23 #include <libmesh/ignore_warnings.h>
24 
25 using namespace libMesh;
26 
27 class BoundaryMeshTest : public CppUnit::TestCase {
32 public:
33  CPPUNIT_TEST_SUITE( BoundaryMeshTest );
34 
35  CPPUNIT_TEST( testMesh );
36 
37  CPPUNIT_TEST_SUITE_END();
38 
39 protected:
40 
45 
46  void build_mesh()
47  {
48  _mesh.reset(new Mesh(*TestCommWorld));
49  _all_boundary_mesh.reset(new Mesh(*TestCommWorld));
50  _left_boundary_mesh.reset(new Mesh(*TestCommWorld));
51  _internal_boundary_mesh.reset(new Mesh(*TestCommWorld));
52 
54  0.2, 0.8, 0.2, 0.7, QUAD9);
55 
56  // We'll need to skip repartitioning with DistributedMesh for now;
57  // otherwise the boundary meshes' interior parents might get
58  // shuffled off to different processors.
59  if (!_mesh->is_serial())
60  {
61  _mesh->skip_partitioning(true);
62  _left_boundary_mesh->skip_partitioning(true);
63  _all_boundary_mesh->skip_partitioning(true);
64  _internal_boundary_mesh->skip_partitioning(true);
65  }
66 
67  // Set subdomain ids for specific elements. This allows us to later
68  // build an internal sideset with respect to a given
69  // subdomain. The element subdomains look like:
70  // ___________________
71  // | 2 | 2 | 2 |
72  // |_____|_____|_____|
73  // | 2 | 2 | 2 |
74  // |_____|_____|_____|
75  // | 2 | 2 | 2 |
76  // |_____|_____|_____|
77  // | 1 | 1 | 2 |
78  // |_____|_____|_____|
79  // | 1 | 1 | 2 |
80  // |_____|_____|_____|
81  //
82  // and we will create an internal sideset along the border between
83  // subdomains 1 and 2.
84 
85  for (MeshBase::element_iterator elem_it =
86  _mesh->active_elements_begin(),
87  elem_end =
88  _mesh->active_elements_end();
89  elem_it != elem_end; ++elem_it)
90  {
91  Elem * elem = *elem_it;
92  if (elem)
93  {
94  const Point c = elem->centroid();
95  if (c(0) < 0.6 &&
96  c(1) < 0.4)
97  elem->subdomain_id() = 1;
98  else
99  elem->subdomain_id() = 2;
100  }
101  }
102 
103  // Get the border of the square
104  _mesh->get_boundary_info().sync(*_all_boundary_mesh);
105 
106  std::set<boundary_id_type> left_id, right_id;
107  left_id.insert(3);
108  right_id.insert(1);
109 
110  // Add the right side of the square to the square; this should
111  // make it a mixed dimension mesh
112  _mesh->get_boundary_info().add_elements(right_id, *_mesh);
113  _mesh->prepare_for_use();
114 
115  // Add the left side of the square to its own boundary mesh.
116  _mesh->get_boundary_info().sync(left_id, *_left_boundary_mesh);
117 
118  // We create an internal sideset ID that does not conflict with
119  // sidesets 0-3 that get created by build_square().
120  boundary_id_type bid = 5;
121 
122  // To test the "relative to" feature, we add the same sides to the
123  // same sideset twice, from elements in subdomain 2 the second
124  // time. These should not show up in the BoundaryMesh, i.e. there
125  // should not be overlapped elems in the BoundaryMesh.
126  BoundaryInfo & bi = _mesh->get_boundary_info();
127 
128  for (MeshBase::element_iterator elem_it =
129  _mesh->active_elements_begin(),
130  elem_end =
131  _mesh->active_elements_end();
132  elem_it != elem_end; ++elem_it)
133  {
134  Elem * elem = *elem_it;
135  if (elem)
136  {
137  const Point c = elem->centroid();
138  if (c(0) < 0.6 &&
139  c(1) < 0.4)
140  {
141  if (c(0) > 0.4)
142  bi.add_side(elem, 1, bid);
143  if (c(1) > 0.3)
144  bi.add_side(elem, 2, bid);
145  }
146  else
147  {
148  if (c(0) < 0.75 &&
149  c(1) < 0.4)
150  bi.add_side(elem, 3, bid);
151  if (c(0) < 0.6 &&
152  c(1) < 0.5)
153  bi.add_side(elem, 0, bid);
154  }
155  }
156  }
157 
158 
159  // Create a BoundaryMesh from the internal sideset relative to subdomain 1.
160  {
161  std::set<boundary_id_type> requested_boundary_ids;
162  requested_boundary_ids.insert(bid);
163  std::set<subdomain_id_type> subdomains_relative_to;
164  subdomains_relative_to.insert(1);
165  _mesh->get_boundary_info().sync(requested_boundary_ids,
166  *_internal_boundary_mesh,
167  subdomains_relative_to);
168  }
169  }
170 
171 public:
172  void setUp()
173  {
174  this->build_mesh();
175  }
176 
177  void testMesh()
178  {
179  // There'd better be 3*5 + 5 elements in the interior plus right
180  // boundary
181  CPPUNIT_ASSERT_EQUAL(static_cast<dof_id_type>(20),
182  _mesh->n_elem());
183 
184  // There'd better be only 7*11 nodes in the interior
185  CPPUNIT_ASSERT_EQUAL(static_cast<dof_id_type>(77),
186  _mesh->n_nodes());
187 
188  // There'd better be only 2*(3+5) elements on the full boundary
189  CPPUNIT_ASSERT_EQUAL(static_cast<dof_id_type>(16),
190  _all_boundary_mesh->n_elem());
191 
192  // There'd better be only 2*2*(3+5) nodes on the full boundary
193  CPPUNIT_ASSERT_EQUAL(static_cast<dof_id_type>(32),
194  _all_boundary_mesh->n_nodes());
195 
196  // There'd better be only 5 elements on the left boundary
197  CPPUNIT_ASSERT_EQUAL(static_cast<dof_id_type>(5),
198  _left_boundary_mesh->n_elem());
199 
200  // There'd better be only 2*5+1 nodes on the left boundary
201  CPPUNIT_ASSERT_EQUAL(static_cast<dof_id_type>(11),
202  _left_boundary_mesh->n_nodes());
203 
204  // There are only four elements in the internal sideset mesh.
205  CPPUNIT_ASSERT_EQUAL(static_cast<dof_id_type>(4),
206  _internal_boundary_mesh->n_elem());
207 
208  // There are 2*n_elem + 1 nodes in the internal sideset mesh.
209  CPPUNIT_ASSERT_EQUAL(static_cast<dof_id_type>(9),
210  _internal_boundary_mesh->n_nodes());
211 
212  this->sanityCheck();
213  }
214 
215  void sanityCheck()
216  {
217  // Sanity check all the elements
219  _mesh->active_elements_begin();
220  const MeshBase::const_element_iterator elem_end =
221  _mesh->active_elements_end();
222  for (; elem_it != elem_end; ++elem_it)
223  {
224  const Elem * elem = *elem_it;
225 
226  const Elem * pip = elem->interior_parent();
227 
228  // On a DistributedMesh we might not be able to see the
229  // interior_parent of a non-local element
230  if (pip == remote_elem)
231  {
232  CPPUNIT_ASSERT(elem->processor_id() != TestCommWorld->rank());
233  continue;
234  }
235 
236  // All the edges should have interior parents; none of the
237  // quads should.
238  if (pip)
239  {
240  CPPUNIT_ASSERT_EQUAL(elem->type(), EDGE3);
241  CPPUNIT_ASSERT_EQUAL(pip->type(), QUAD9);
242  CPPUNIT_ASSERT_EQUAL(pip->level(), elem->level());
243 
244  // We only added right edges
245  CPPUNIT_ASSERT_DOUBLES_EQUAL(elem->centroid()(0), 0.8,
247  }
248  else
249  {
250  CPPUNIT_ASSERT_EQUAL(elem->type(), QUAD9);
251  }
252  }
253 
254  MeshBase::const_element_iterator left_bdy_elem_it =
255  _left_boundary_mesh->active_elements_begin();
256  const MeshBase::const_element_iterator left_bdy_elem_end =
257  _left_boundary_mesh->active_elements_end();
258  for (; left_bdy_elem_it != left_bdy_elem_end; ++left_bdy_elem_it)
259  {
260  const Elem * elem = *left_bdy_elem_it;
261 
262  CPPUNIT_ASSERT_EQUAL(elem->type(), EDGE3);
263 
264  const Elem * pip = elem->interior_parent();
265 
266  // On a DistributedMesh we might not be able to see the
267  // interior_parent of a non-local element
268  if (pip == remote_elem)
269  {
270  CPPUNIT_ASSERT(elem->processor_id() != TestCommWorld->rank());
271  continue;
272  }
273 
274  // All the edges should have interior parents
275  CPPUNIT_ASSERT(pip);
276  CPPUNIT_ASSERT_EQUAL(pip->type(), QUAD9);
277  CPPUNIT_ASSERT_EQUAL(pip->level(), elem->level());
278 
279  // We only added left edges
280  CPPUNIT_ASSERT_DOUBLES_EQUAL(elem->centroid()(0), 0.2,
282  }
283 
284 
285  MeshBase::const_element_iterator all_bdy_elem_it =
286  _left_boundary_mesh->active_elements_begin();
287  const MeshBase::const_element_iterator all_bdy_elem_end =
288  _left_boundary_mesh->active_elements_end();
289  for (; all_bdy_elem_it != all_bdy_elem_end; ++all_bdy_elem_it)
290  {
291  const Elem * elem = *all_bdy_elem_it;
292 
293  CPPUNIT_ASSERT_EQUAL(elem->type(), EDGE3);
294 
295  const Elem * pip = elem->interior_parent();
296 
297  // On a DistributedMesh we might not be able to see the
298  // interior_parent of a non-local element
299  if (pip == remote_elem)
300  {
301  CPPUNIT_ASSERT(elem->processor_id() != TestCommWorld->rank());
302  continue;
303  }
304 
305  // All the edges should have interior parents
306  CPPUNIT_ASSERT(pip);
307  CPPUNIT_ASSERT_EQUAL(pip->type(), QUAD9);
308  CPPUNIT_ASSERT_EQUAL(pip->level(), elem->level());
309  }
310 
311 
312  // Sanity check for the internal sideset mesh.
314  internal_elem_it = _internal_boundary_mesh->active_elements_begin(),
315  internal_elem_end = _internal_boundary_mesh->active_elements_end();
316 
317  for (; internal_elem_it != internal_elem_end; ++internal_elem_it)
318  {
319  const Elem * elem = *internal_elem_it;
320 
321  CPPUNIT_ASSERT_EQUAL(elem->type(), EDGE3);
322 
323  // All of the elements in the internal sideset mesh should
324  // have the same subdomain id as the parent Elems (i.e. 1)
325  // they came from.
326  CPPUNIT_ASSERT_EQUAL(static_cast<subdomain_id_type>(1),
327  elem->subdomain_id());
328  }
329  }
330 
331 };
332 
334 
335 #ifdef LIBMESH_ENABLE_AMR
336 
344 public:
345  CPPUNIT_TEST_SUITE( BoundaryRefinedMeshTest );
346 
347  CPPUNIT_TEST( testMesh );
348 
349  CPPUNIT_TEST_SUITE_END();
350 
351  // Yes, this is necessary. Somewhere in those macros is a protected/private
352 public:
353 
354  void setUp()
355  {
356  this->build_mesh();
357 
358  // Need to refine interior mesh before separate boundary meshes,
359  // if we want to get interior_parent links right.
360  MeshRefinement(*_mesh).uniformly_refine(1);
361  MeshRefinement(*_left_boundary_mesh).uniformly_refine(1);
362  MeshRefinement(*_all_boundary_mesh).uniformly_refine(1);
363  }
364 
365  void testMesh()
366  {
367  // There'd better be 3*5*4 + 5*2 active elements in the interior
368  // plus right boundary
369  CPPUNIT_ASSERT_EQUAL(static_cast<dof_id_type>(70),
370  _mesh->n_active_elem());
371 
372  // Plus the original 20 now-inactive elements
373  CPPUNIT_ASSERT_EQUAL(static_cast<dof_id_type>(90),
374  _mesh->n_elem());
375 
376  // There'd better be only 13*21 nodes in the interior
377  CPPUNIT_ASSERT_EQUAL(static_cast<dof_id_type>(273),
378  _mesh->n_nodes());
379 
380  // There'd better be only 2*2*(3+5) active elements on the full boundary
381  CPPUNIT_ASSERT_EQUAL(static_cast<dof_id_type>(32),
382  _all_boundary_mesh->n_active_elem());
383 
384  // Plus the original 16 now-inactive elements
385  CPPUNIT_ASSERT_EQUAL(static_cast<dof_id_type>(48),
386  _all_boundary_mesh->n_elem());
387 
388  // There'd better be only 2*2*2*(3+5) nodes on the full boundary
389  CPPUNIT_ASSERT_EQUAL(static_cast<dof_id_type>(64),
390  _all_boundary_mesh->n_nodes());
391 
392  // There'd better be only 2*5 active elements on the left boundary
393  CPPUNIT_ASSERT_EQUAL(static_cast<dof_id_type>(10),
394  _left_boundary_mesh->n_active_elem());
395 
396  // Plus the original 5 now-inactive elements
397  CPPUNIT_ASSERT_EQUAL(static_cast<dof_id_type>(15),
398  _left_boundary_mesh->n_elem());
399 
400  // There'd better be only 2*2*5+1 nodes on the left boundary
401  CPPUNIT_ASSERT_EQUAL(static_cast<dof_id_type>(21),
402  _left_boundary_mesh->n_nodes());
403 
404  this->sanityCheck();
405  }
406 
407 };
408 
410 
411 #endif
The definition of the element_iterator struct.
Definition: mesh_base.h:1476
UniquePtr< Mesh > _left_boundary_mesh
Definition: boundary_mesh.C:43
virtual ElemType type() const =0
libMesh::Parallel::Communicator * TestCommWorld
Definition: driver.C:28
The definition of the const_element_iterator struct.
Definition: mesh_base.h:1494
const Elem * interior_parent() const
Definition: elem.C:951
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
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.
static const Real TOLERANCE
The libMesh namespace provides an interface to certain functionality in the library.
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
This is the MeshRefinement class.
int8_t boundary_id_type
Definition: id_types.h:51
UniquePtr< Mesh > _mesh
Definition: boundary_mesh.C:41
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:56
CPPUNIT_TEST_SUITE_REGISTRATION(BoundaryMeshTest)
subdomain_id_type subdomain_id() const
Definition: elem.h:1951
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...
unsigned int level() const
Definition: elem.h:2388
virtual Point centroid() const
Definition: elem.C:446
unsigned int rank() const
Definition: parallel.h:724
UniquePtr< Mesh > _internal_boundary_mesh
Definition: boundary_mesh.C:44
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
UniquePtr< Mesh > _all_boundary_mesh
Definition: boundary_mesh.C:42
processor_id_type processor_id() const
Definition: dof_object.h:694
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.
const RemoteElem * remote_elem
Definition: remote_elem.C:57