libMesh
slit_mesh_test.C
Go to the documentation of this file.
1 // Ignore unused parameter warnings coming from cppuint 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/equation_systems.h>
8 #include <libmesh/mesh.h>
9 #include <libmesh/mesh_generation.h>
10 #include <libmesh/edge_edge2.h>
11 #include <libmesh/face_quad4.h>
12 #include <libmesh/cell_hex8.h>
13 #include <libmesh/dof_map.h>
14 #include <libmesh/linear_implicit_system.h>
15 #include <libmesh/mesh_refinement.h>
16 
17 #include "test_comm.h"
18 
19 // THE CPPUNIT_TEST_SUITE_END macro expands to code that involves
20 // std::auto_ptr, which in turn produces -Wdeprecated-declarations
21 // warnings. These can be ignored in GCC as long as we wrap the
22 // offending code in appropriate pragmas. We can't get away with a
23 // single ignore_warnings.h inclusion at the beginning of this file,
24 // since the libmesh headers pull in a restore_warnings.h at some
25 // point. We also don't bother restoring warnings at the end of this
26 // file since it's not a header.
27 #include <libmesh/ignore_warnings.h>
28 
29 using namespace libMesh;
30 
31 class SlitFunc : public FEMFunctionBase<Number>
32 {
33 public:
34 
35  SlitFunc() {}
36 
37  ~SlitFunc () {}
38 
39  virtual void init_context (const FEMContext &) libmesh_override {}
40 
42  clone () const libmesh_override
43  {
45  }
46 
47  virtual Number operator() (const FEMContext & c,
48  const Point & p,
49  const Real /*time*/ = 0.) libmesh_override
50  {
51  const Real & x = p(0);
52  const Real & y = p(1);
53  const Point centroid = c.get_elem().centroid();
54  const Real sign = centroid(1)/std::abs(centroid(1));
55 
56  return (1 - std::abs(1-x)) * (1-std::abs(y)) * sign;
57  }
58 
59  virtual void operator() (const FEMContext & c,
60  const Point & p,
61  const Real time,
62  DenseVector<Number> & output) libmesh_override
63  {
64  for (unsigned int i=0; i != output.size(); ++i)
65  output(i) = (*this)(c, p, time);
66  }
67 };
68 
69 
70 
71 
72 
73 
74 
75 class SlitMeshTest : public CppUnit::TestCase {
80 public:
81  CPPUNIT_TEST_SUITE( SlitMeshTest );
82 
83  CPPUNIT_TEST( testMesh );
84 
85  CPPUNIT_TEST_SUITE_END();
86 
87 protected:
88 
90 
91  void build_mesh()
92  {
93  _mesh = new Mesh(*TestCommWorld);
94 
95  // (0,1) (1,1) (2,1)
96  // x---------------x---------------x
97  // | | |
98  // | | |
99  // | | |
100  // | | |
101  // | | |
102  // x---------------x---------------x
103  // (0,0) (1,0) (2,0)
104  // | | |
105  // | | |
106  // | | |
107  // | | |
108  // x---------------x---------------x
109  // (0,-1) (1,-1) (2,-1)
110 
111  _mesh->set_mesh_dimension(2);
112 
113  _mesh->add_point( Point(0.0, 0.0), 0 );
114  _mesh->add_point( Point(1.0, 0.0), 1 );
115  _mesh->add_point( Point(1.0, 1.0), 2 );
116  _mesh->add_point( Point(0.0, 1.0), 3 );
117  _mesh->add_point( Point(0.0,-1.0), 4 );
118  _mesh->add_point( Point(1.0,-1.0), 5 );
119  _mesh->add_point( Point(1.0, 0.0), 6 );
120  _mesh->add_point( Point(2.0, 0.0), 7 );
121  _mesh->add_point( Point(2.0, 1.0), 8 );
122  _mesh->add_point( Point(2.0,-1.0), 9 );
123 
124  {
125  Elem* elem_top_left = new Quad4;
126  elem_top_left->set_node(0) = _mesh->node_ptr(0);
127  elem_top_left->set_node(1) = _mesh->node_ptr(1);
128  elem_top_left->set_node(2) = _mesh->node_ptr(2);
129  elem_top_left->set_node(3) = _mesh->node_ptr(3);
130  elem_top_left->set_id() = 0;
131  _mesh->add_elem(elem_top_left);
132 
133  Elem* elem_bottom_left = new Quad4;
134  elem_bottom_left->set_node(0) = _mesh->node_ptr(4);
135  elem_bottom_left->set_node(1) = _mesh->node_ptr(5);
136  elem_bottom_left->set_node(2) = _mesh->node_ptr(6);
137  elem_bottom_left->set_node(3) = _mesh->node_ptr(0);
138  elem_bottom_left->set_id() = 1;
139  _mesh->add_elem(elem_bottom_left);
140 
141  Elem* elem_top_right = new Quad4;
142  elem_top_right->set_node(0) = _mesh->node_ptr(1);
143  elem_top_right->set_node(1) = _mesh->node_ptr(7);
144  elem_top_right->set_node(2) = _mesh->node_ptr(8);
145  elem_top_right->set_node(3) = _mesh->node_ptr(2);
146  elem_top_right->set_id() = 2;
147  _mesh->add_elem(elem_top_right);
148 
149  Elem* elem_bottom_right = new Quad4;
150  elem_bottom_right->set_node(0) = _mesh->node_ptr(5);
151  elem_bottom_right->set_node(1) = _mesh->node_ptr(9);
152  elem_bottom_right->set_node(2) = _mesh->node_ptr(7);
153  elem_bottom_right->set_node(3) = _mesh->node_ptr(6);
154  elem_bottom_right->set_id() = 3;
155  _mesh->add_elem(elem_bottom_right);
156  }
157 
158  // libMesh shouldn't renumber, or our based-on-initial-id
159  // assertions later may fail.
160  _mesh->allow_renumbering(false);
161 
162  _mesh->prepare_for_use();
163  }
164 
165 public:
166  void setUp()
167  {
168  this->build_mesh();
169  }
170 
171  void tearDown()
172  {
173  delete _mesh;
174  }
175 
176  void testMesh()
177  {
178  // There'd better be 4 elements
179  CPPUNIT_ASSERT_EQUAL( (dof_id_type)4, _mesh->n_elem() );
180 
181  // There'd better still be a full 10 nodes
182  CPPUNIT_ASSERT_EQUAL( (dof_id_type)10, _mesh->n_nodes() );
183 
184  /* The middle nodes should still be distinct between the top and
185  * bottom elements */
186  if (_mesh->query_elem_ptr(0) && _mesh->query_elem_ptr(1))
187  CPPUNIT_ASSERT( _mesh->elem_ref(0).node_id(1) != _mesh->elem_ref(1).node_id(2) );
188  if (_mesh->query_elem_ptr(2) && _mesh->query_elem_ptr(3))
189  CPPUNIT_ASSERT( _mesh->elem_ref(2).node_id(0) != _mesh->elem_ref(3).node_id(3) );
190 
191  /* The middle nodes should still be shared between left and right
192  * elements on top and bottom */
193  if (_mesh->query_elem_ptr(0) && _mesh->query_elem_ptr(2))
194  CPPUNIT_ASSERT_EQUAL( _mesh->elem_ref(0).node_id(1),
195  _mesh->elem_ref(2).node_id(0) );
196  if (_mesh->query_elem_ptr(1) && _mesh->query_elem_ptr(3))
197  CPPUNIT_ASSERT_EQUAL( _mesh->elem_ref(1).node_id(2),
198  _mesh->elem_ref(3).node_id(3) );
199  }
200 
201 };
202 
210 public:
211  CPPUNIT_TEST_SUITE( SlitMeshRefinedMeshTest );
212 
213  CPPUNIT_TEST( testMesh );
214 
215  CPPUNIT_TEST_SUITE_END();
216 
217  // Yes, this is necessary. Somewhere in those macros is a protected/private
218 public:
219 
220  void setUp()
221  {
222  this->build_mesh();
223 
224 #ifdef LIBMESH_ENABLE_AMR
225  MeshRefinement(*_mesh).uniformly_refine(1);
226 #endif
227  }
228 
229  void testMesh()
230  {
231 #ifdef LIBMESH_ENABLE_AMR
232  // We should have 20 total and 16 active elements.
233  CPPUNIT_ASSERT_EQUAL( (dof_id_type)20, _mesh->n_elem() );
234  CPPUNIT_ASSERT_EQUAL( (dof_id_type)16, _mesh->n_active_elem() );
235 
236  // We should have 28 nodes, not 25 or 26
237  CPPUNIT_ASSERT_EQUAL( (dof_id_type)28, _mesh->n_nodes() );
238 #endif
239  }
240 };
241 
248 public:
249  CPPUNIT_TEST_SUITE( SlitMeshRefinedSystemTest );
250 
251  CPPUNIT_TEST( testMesh );
252 
253  CPPUNIT_TEST( testSystem );
254 
255 #ifdef LIBMESH_ENABLE_UNIQUE_ID
256  CPPUNIT_TEST( testRestart );
257 #endif
258 
259  CPPUNIT_TEST_SUITE_END();
260 
261 protected:
262 
265 
266 public:
267 
268  void setUp()
269  {
270  this->build_mesh();
271 
272  // libMesh *should* renumber now, or a DistributedMesh might not
273  // have contiguous ids, which is a requirement to write xda files.
274  _mesh->allow_renumbering(true);
275 
276  _es = new EquationSystems(*_mesh);
277  _sys = &_es->add_system<System> ("SimpleSystem");
278  _sys->add_variable("u", FIRST);
279 
280  _es->init();
281  SlitFunc slitfunc;
282  _sys->project_solution(&slitfunc);
283 
284 #ifdef LIBMESH_ENABLE_AMR
285  MeshRefinement(*_mesh).uniformly_refine(1);
286  _es->reinit();
287  MeshRefinement(*_mesh).uniformly_refine(1);
288  _es->reinit();
289 #endif
290  }
291 
292  void tearDown()
293  {
294  delete _es;
295  // _sys is owned by _es
296  delete _mesh;
297  }
298 
299  void testMesh()
300  {
301 #ifdef LIBMESH_ENABLE_AMR
302  // We should have 84 total and 64 active elements.
303  CPPUNIT_ASSERT_EQUAL( (dof_id_type)(4+16+64), _mesh->n_elem() );
304  CPPUNIT_ASSERT_EQUAL( (dof_id_type)64, _mesh->n_active_elem() );
305 
306  // We should have 88 nodes
307  CPPUNIT_ASSERT_EQUAL( (dof_id_type)88, _mesh->n_nodes() );
308 #endif
309  }
310 
311  void testSystem()
312  {
313  SlitFunc slitfunc;
314 
315  unsigned int dim = 2;
316 
317  CPPUNIT_ASSERT_EQUAL( _sys->n_vars(), 1u );
318 
319  FEMContext context(*_sys);
320  FEBase * fe = NULL;
321  context.get_element_fe( 0, fe, dim );
322  const std::vector<Point> & xyz = fe->get_xyz();
323  fe->get_phi();
324 
326  _mesh->active_local_elements_begin();
327  const MeshBase::const_element_iterator end_el =
328  _mesh->active_local_elements_end();
329 
330  for (; el != end_el; ++el)
331  {
332  const Elem * elem = *el;
333  context.pre_fe_reinit(*_sys, elem);
334  context.elem_fe_reinit();
335 
336  const unsigned int n_qp = xyz.size();
337 
338  for (unsigned int qp=0; qp != n_qp; ++qp)
339  {
340  const Number exact_val = slitfunc(context, xyz[qp]);
341 
342  const Number discrete_val = context.interior_value(0, qp);
343 
344  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(exact_val),
345  libmesh_real(discrete_val),
347  }
348  }
349  }
350 
351  void testRestart()
352  {
353  SlitFunc slitfunc;
354 
355  _mesh->write("slit_mesh.xda");
356  _es->write("slit_solution.xda",
359 
360  Mesh mesh2(*TestCommWorld);
361  mesh2.read("slit_mesh.xda");
362  EquationSystems es2(mesh2);
363  es2.read("slit_solution.xda");
364 
365  System & sys2 = es2.get_system<System> ("SimpleSystem");
366 
367  unsigned int dim = 2;
368 
369  CPPUNIT_ASSERT_EQUAL( sys2.n_vars(), 1u );
370 
371  FEMContext context(sys2);
372  FEBase * fe = NULL;
373  context.get_element_fe( 0, fe, dim );
374  const std::vector<Point> & xyz = fe->get_xyz();
375  fe->get_phi();
376 
377  // While we're in the middle of a unique id based test case, let's
378  // make sure our unique ids were all read in correctly too.
379  UniquePtr<PointLocatorBase> locator = _mesh->sub_point_locator();
380 
381  if (!_mesh->is_serial())
382  locator->enable_out_of_mesh_mode();
383 
384  for (const auto & elem : mesh2.active_local_element_ptr_range())
385  {
386  const Elem * mesh1_elem = (*locator)(elem->centroid());
387  if (mesh1_elem)
388  {
389  CPPUNIT_ASSERT_EQUAL( elem->unique_id(),
390  mesh1_elem->unique_id() );
391 
392  for (unsigned int n=0; n != elem->n_nodes(); ++n)
393  {
394  const Node & node = elem->node_ref(n);
395  const Node & mesh1_node = mesh1_elem->node_ref(n);
396  CPPUNIT_ASSERT_EQUAL( node.unique_id(),
397  mesh1_node.unique_id() );
398  }
399  }
400 
401  context.pre_fe_reinit(sys2, elem);
402  context.elem_fe_reinit();
403 
404  const unsigned int n_qp = xyz.size();
405 
406  for (unsigned int qp=0; qp != n_qp; ++qp)
407  {
408  const Number exact_val = slitfunc(context, xyz[qp]);
409 
410  const Number discrete_val = context.interior_value(0, qp);
411 
412  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(exact_val),
413  libmesh_real(discrete_val),
415  }
416  }
417  }
418 };
419 
T libmesh_real(T a)
double abs(double a)
This is the EquationSystems class.
virtual void read(const std::string &name, void *mesh_data=libmesh_nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false) libmesh_override
Reads the file specified by name.
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:1941
A Node is like a Point, but with more information.
Definition: node.h:52
void write(const std::string &name, const XdrMODE, const unsigned int write_flags=(WRITE_DATA), bool partition_agnostic=true) const
Write the systems to disk using the XDR data format.
virtual dof_id_type n_elem() const libmesh_override
The QUAD4 is an element in 2D composed of 4 nodes.
Definition: face_quad4.h:49
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:749
libMesh::Parallel::Communicator * TestCommWorld
Definition: driver.C:28
The definition of the const_element_iterator struct.
Definition: mesh_base.h:1494
unsigned int dim
virtual SimpleRange< element_iterator > active_local_element_ptr_range() libmesh_override
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1101
CPPUNIT_TEST_SUITE_REGISTRATION(SlitMeshTest)
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
static const Real TOLERANCE
void read(const std::string &name, const XdrMODE, const unsigned int read_flags=(READ_HEADER|READ_DATA), bool partition_agnostic=true)
Read & initialize the systems from disk using the XDR data format.
void project_solution(FunctionBase< Number > *f, FunctionBase< Gradient > *g=libmesh_nullptr) const
Projects arbitrary functions onto the current solution.
The libMesh namespace provides an interface to certain functionality in the library.
const Elem & get_elem() const
Accessor for current Elem object.
Definition: fem_context.h:871
dof_id_type & set_id()
Definition: dof_object.h:641
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
PetscErrorCode Vec x
This is the MeshRefinement class.
virtual const Elem * query_elem_ptr(const dof_id_type i) const libmesh_override
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
virtual dof_id_type n_nodes() const libmesh_override
virtual void reinit()
Reinitialize all the systems.
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:208
const Node & node_ref(const unsigned int i) const
Definition: elem.h:1896
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
virtual System & add_system(const std::string &system_type, const std::string &name)
Add the system of type system_type named name to the systems array.
virtual const Node * node_ptr(const dof_id_type i) const libmesh_override
void prepare_for_use(const bool skip_renumber_nodes_and_elements=false, const bool skip_find_neighbors=false)
Prepare a newly created (or read) mesh for use.
Definition: mesh_base.C:174
void set_mesh_dimension(unsigned char d)
Resets the logical dimension of the mesh.
Definition: mesh_base.h:199
virtual Elem * add_elem(Elem *e) libmesh_override
Add elem e to the end of the element array.
virtual const Elem & elem_ref(const dof_id_type i) const
Definition: mesh_base.h:490
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void build_mesh()
virtual Point centroid() const
Definition: elem.C:446
virtual UniquePtr< FEMFunctionBase< Number > > clone() const
const T_sys & get_system(const std::string &name) const
unsigned int n_vars() const
Definition: system.h:2086
Real size() const
Definition: type_vector.h:898
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:1831
virtual void init()
Initialize all the systems.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
virtual void init_context(const FEMContext &)
Prepares a context object for use.
unique_id_type unique_id() const
Definition: dof_object.h:649
FEMFunctionBase is a base class from which users can derive in order to define "function-like" object...
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
const std::vector< Point > & get_xyz() const
Definition: fe_abstract.h:230
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) libmesh_override
functions for adding /deleting nodes elements.
This class forms the foundation from which generic finite elements may be derived.
uint8_t dof_id_type
Definition: id_types.h:64
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.