libMesh
slit_mesh_test.C
Go to the documentation of this file.
1 #include <libmesh/equation_systems.h>
2 #include <libmesh/mesh.h>
3 #include <libmesh/mesh_generation.h>
4 #include <libmesh/edge_edge2.h>
5 #include <libmesh/face_quad4.h>
6 #include <libmesh/cell_hex8.h>
7 #include <libmesh/dof_map.h>
8 #include <libmesh/linear_implicit_system.h>
9 #include <libmesh/mesh_refinement.h>
10 
11 #include "test_comm.h"
12 #include "libmesh_cppunit.h"
13 
14 
15 using namespace libMesh;
16 
17 class SlitFunc : public FEMFunctionBase<Number>
18 {
19 public:
20 
21  SlitFunc() {}
22 
23  ~SlitFunc () {}
24 
25  virtual void init_context (const FEMContext &) override {}
26 
27  virtual std::unique_ptr<FEMFunctionBase<Number>>
28  clone () const override
29  {
30  return libmesh_make_unique<SlitFunc>();
31  }
32 
33  virtual Number operator() (const FEMContext & c,
34  const Point & p,
35  const Real /*time*/ = 0.) override
36  {
37  const Real & x = p(0);
38  const Real & y = p(1);
39  const Point centroid = c.get_elem().centroid();
40  const Real sign = centroid(1)/std::abs(centroid(1));
41 
42  return (1 - std::abs(1-x)) * (1-std::abs(y)) * sign;
43  }
44 
45  virtual void operator() (const FEMContext & c,
46  const Point & p,
47  const Real time,
48  DenseVector<Number> & output) override
49  {
50  for (unsigned int i=0; i != output.size(); ++i)
51  output(i) = (*this)(c, p, time);
52  }
53 };
54 
55 
56 
57 
58 
59 
60 
61 class SlitMeshTest : public CppUnit::TestCase {
66 public:
67  CPPUNIT_TEST_SUITE( SlitMeshTest );
68 
69 #if LIBMESH_DIM > 1
70  CPPUNIT_TEST( testMesh );
71 #endif
72 
73  CPPUNIT_TEST_SUITE_END();
74 
75 protected:
76 
77  Mesh* _mesh;
78 
79  void build_mesh()
80  {
81  _mesh = new Mesh(*TestCommWorld);
82 
83  // (0,1) (1,1) (2,1)
84  // x---------------x---------------x
85  // | | |
86  // | | |
87  // | | |
88  // | | |
89  // | | |
90  // x---------------x---------------x
91  // (0,0) (1,0) (2,0)
92  // | | |
93  // | | |
94  // | | |
95  // | | |
96  // x---------------x---------------x
97  // (0,-1) (1,-1) (2,-1)
98 
99  _mesh->set_mesh_dimension(2);
100 
101  _mesh->add_point( Point(0.0, 0.0), 0 );
102  _mesh->add_point( Point(1.0, 0.0), 1 );
103  _mesh->add_point( Point(1.0, 1.0), 2 );
104  _mesh->add_point( Point(0.0, 1.0), 3 );
105  _mesh->add_point( Point(0.0,-1.0), 4 );
106  _mesh->add_point( Point(1.0,-1.0), 5 );
107  _mesh->add_point( Point(1.0, 0.0), 6 );
108  _mesh->add_point( Point(2.0, 0.0), 7 );
109  _mesh->add_point( Point(2.0, 1.0), 8 );
110  _mesh->add_point( Point(2.0,-1.0), 9 );
111 
112  {
113  Elem* elem_top_left = new Quad4;
114  elem_top_left->set_node(0) = _mesh->node_ptr(0);
115  elem_top_left->set_node(1) = _mesh->node_ptr(1);
116  elem_top_left->set_node(2) = _mesh->node_ptr(2);
117  elem_top_left->set_node(3) = _mesh->node_ptr(3);
118  elem_top_left->set_id() = 0;
119  _mesh->add_elem(elem_top_left);
120 
121  Elem* elem_bottom_left = new Quad4;
122  elem_bottom_left->set_node(0) = _mesh->node_ptr(4);
123  elem_bottom_left->set_node(1) = _mesh->node_ptr(5);
124  elem_bottom_left->set_node(2) = _mesh->node_ptr(6);
125  elem_bottom_left->set_node(3) = _mesh->node_ptr(0);
126  elem_bottom_left->set_id() = 1;
127  _mesh->add_elem(elem_bottom_left);
128 
129  Elem* elem_top_right = new Quad4;
130  elem_top_right->set_node(0) = _mesh->node_ptr(1);
131  elem_top_right->set_node(1) = _mesh->node_ptr(7);
132  elem_top_right->set_node(2) = _mesh->node_ptr(8);
133  elem_top_right->set_node(3) = _mesh->node_ptr(2);
134  elem_top_right->set_id() = 2;
135  _mesh->add_elem(elem_top_right);
136 
137  Elem* elem_bottom_right = new Quad4;
138  elem_bottom_right->set_node(0) = _mesh->node_ptr(5);
139  elem_bottom_right->set_node(1) = _mesh->node_ptr(9);
140  elem_bottom_right->set_node(2) = _mesh->node_ptr(7);
141  elem_bottom_right->set_node(3) = _mesh->node_ptr(6);
142  elem_bottom_right->set_id() = 3;
143  _mesh->add_elem(elem_bottom_right);
144  }
145 
146  // libMesh shouldn't renumber, or our based-on-initial-id
147  // assertions later may fail.
148  _mesh->allow_renumbering(false);
149 
150  _mesh->prepare_for_use();
151  }
152 
153 public:
154  void setUp()
155  {
156 #if LIBMESH_DIM > 1
157  this->build_mesh();
158 #endif
159  }
160 
161  void tearDown()
162  {
163  delete _mesh;
164  }
165 
166  void testMesh()
167  {
168  // There'd better be 4 elements
169  CPPUNIT_ASSERT_EQUAL( (dof_id_type)4, _mesh->n_elem() );
170 
171  // There'd better still be a full 10 nodes
172  CPPUNIT_ASSERT_EQUAL( (dof_id_type)10, _mesh->n_nodes() );
173 
174  /* The middle nodes should still be distinct between the top and
175  * bottom elements */
176  if (_mesh->query_elem_ptr(0) && _mesh->query_elem_ptr(1))
177  CPPUNIT_ASSERT( _mesh->elem_ref(0).node_id(1) != _mesh->elem_ref(1).node_id(2) );
178  if (_mesh->query_elem_ptr(2) && _mesh->query_elem_ptr(3))
179  CPPUNIT_ASSERT( _mesh->elem_ref(2).node_id(0) != _mesh->elem_ref(3).node_id(3) );
180 
181  /* The middle nodes should still be shared between left and right
182  * elements on top and bottom */
183  if (_mesh->query_elem_ptr(0) && _mesh->query_elem_ptr(2))
184  CPPUNIT_ASSERT_EQUAL( _mesh->elem_ref(0).node_id(1),
185  _mesh->elem_ref(2).node_id(0) );
186  if (_mesh->query_elem_ptr(1) && _mesh->query_elem_ptr(3))
187  CPPUNIT_ASSERT_EQUAL( _mesh->elem_ref(1).node_id(2),
188  _mesh->elem_ref(3).node_id(3) );
189  }
190 
191 };
192 
193 class SlitMeshRefinedMeshTest : public SlitMeshTest {
200 public:
201  CPPUNIT_TEST_SUITE( SlitMeshRefinedMeshTest );
202 
203 #if LIBMESH_DIM > 1
204  CPPUNIT_TEST( testMesh );
205 #endif
206 
207  CPPUNIT_TEST_SUITE_END();
208 
209  // Yes, this is necessary. Somewhere in those macros is a protected/private
210 public:
211 
212  void setUp()
213  {
214 #if LIBMESH_DIM > 1
215  this->build_mesh();
216 
217 #ifdef LIBMESH_ENABLE_AMR
218  MeshRefinement(*_mesh).uniformly_refine(1);
219 #endif
220 #endif
221  }
222 
223  void testMesh()
224  {
225 #ifdef LIBMESH_ENABLE_AMR
226  // We should have 20 total and 16 active elements.
227  CPPUNIT_ASSERT_EQUAL( (dof_id_type)20, _mesh->n_elem() );
228  CPPUNIT_ASSERT_EQUAL( (dof_id_type)16, _mesh->n_active_elem() );
229 
230  // We should have 28 nodes, not 25 or 26
231  CPPUNIT_ASSERT_EQUAL( (dof_id_type)28, _mesh->n_nodes() );
232 #endif
233  }
234 };
235 
236 class SlitMeshRefinedSystemTest : public SlitMeshTest {
242 public:
243  CPPUNIT_TEST_SUITE( SlitMeshRefinedSystemTest );
244 
245 #if LIBMESH_DIM > 1
246  CPPUNIT_TEST( testMesh );
247 
248  CPPUNIT_TEST( testSystem );
249 
250 #ifdef LIBMESH_ENABLE_UNIQUE_ID
251  CPPUNIT_TEST( testRestart );
252 #endif
253 #endif
254 
255  CPPUNIT_TEST_SUITE_END();
256 
257 protected:
258 
259  System* _sys;
260  EquationSystems* _es;
261 
262 public:
263 
264  void setUp()
265  {
266 #if LIBMESH_DIM > 1
267  this->build_mesh();
268 
269  // libMesh *should* renumber now, or a DistributedMesh might not
270  // have contiguous ids, which is a requirement to write xda files.
271  _mesh->allow_renumbering(true);
272 
273  _es = new EquationSystems(*_mesh);
274  _sys = &_es->add_system<System> ("SimpleSystem");
275  _sys->add_variable("u", FIRST);
276 
277  _es->init();
278  SlitFunc slitfunc;
279  _sys->project_solution(&slitfunc);
280 
281 #ifdef LIBMESH_ENABLE_AMR
282  MeshRefinement(*_mesh).uniformly_refine(1);
283  _es->reinit();
284  MeshRefinement(*_mesh).uniformly_refine(1);
285  _es->reinit();
286 #endif
287 #endif
288  }
289 
290  void tearDown()
291  {
292  delete _es;
293  // _sys is owned by _es
294  delete _mesh;
295  }
296 
297  void testMesh()
298  {
299 #ifdef LIBMESH_ENABLE_AMR
300  // We should have 84 total and 64 active elements.
301  CPPUNIT_ASSERT_EQUAL( (dof_id_type)(4+16+64), _mesh->n_elem() );
302  CPPUNIT_ASSERT_EQUAL( (dof_id_type)64, _mesh->n_active_elem() );
303 
304  // We should have 88 nodes
305  CPPUNIT_ASSERT_EQUAL( (dof_id_type)88, _mesh->n_nodes() );
306 #endif
307  }
308 
309  void testSystem()
310  {
311  SlitFunc slitfunc;
312 
313  unsigned int dim = 2;
314 
315  CPPUNIT_ASSERT_EQUAL( _sys->n_vars(), 1u );
316 
317  FEMContext context(*_sys);
318  FEBase * fe = NULL;
319  context.get_element_fe( 0, fe, dim );
320  const std::vector<Point> & xyz = fe->get_xyz();
321  fe->get_phi();
322 
323  for (const auto & elem : _mesh->active_local_element_ptr_range())
324  {
325  context.pre_fe_reinit(*_sys, elem);
326  context.elem_fe_reinit();
327 
328  const unsigned int n_qp = xyz.size();
329 
330  for (unsigned int qp=0; qp != n_qp; ++qp)
331  {
332  const Number exact_val = slitfunc(context, xyz[qp]);
333 
334  const Number discrete_val = context.interior_value(0, qp);
335 
336  LIBMESH_ASSERT_FP_EQUAL(libmesh_real(exact_val),
337  libmesh_real(discrete_val),
338  TOLERANCE*TOLERANCE);
339  }
340  }
341  }
342 
343  void testRestart()
344  {
345  SlitFunc slitfunc;
346 
347  _mesh->write("slit_mesh.xda");
348  _es->write("slit_solution.xda",
349  EquationSystems::WRITE_DATA |
350  EquationSystems::WRITE_SERIAL_FILES);
351 
352  Mesh mesh2(*TestCommWorld);
353  mesh2.read("slit_mesh.xda");
354  EquationSystems es2(mesh2);
355  es2.read("slit_solution.xda");
356 
357  System & sys2 = es2.get_system<System> ("SimpleSystem");
358 
359  unsigned int dim = 2;
360 
361  CPPUNIT_ASSERT_EQUAL( sys2.n_vars(), 1u );
362 
363  FEMContext context(sys2);
364  FEBase * fe = NULL;
365  context.get_element_fe( 0, fe, dim );
366  const std::vector<Point> & xyz = fe->get_xyz();
367  fe->get_phi();
368 
369  // While we're in the middle of a unique id based test case, let's
370  // make sure our unique ids were all read in correctly too.
371  std::unique_ptr<PointLocatorBase> locator = _mesh->sub_point_locator();
372 
373  if (!_mesh->is_serial())
374  locator->enable_out_of_mesh_mode();
375 
376  for (const auto & elem : mesh2.active_local_element_ptr_range())
377  {
378  const Elem * mesh1_elem = (*locator)(elem->centroid());
379  if (mesh1_elem)
380  {
381  CPPUNIT_ASSERT_EQUAL( elem->unique_id(),
382  mesh1_elem->unique_id() );
383 
384  for (unsigned int n=0; n != elem->n_nodes(); ++n)
385  {
386  const Node & node = elem->node_ref(n);
387  const Node & mesh1_node = mesh1_elem->node_ref(n);
388  CPPUNIT_ASSERT_EQUAL( node.unique_id(),
389  mesh1_node.unique_id() );
390  }
391  }
392 
393  context.pre_fe_reinit(sys2, elem);
394  context.elem_fe_reinit();
395 
396  const unsigned int n_qp = xyz.size();
397 
398  for (unsigned int qp=0; qp != n_qp; ++qp)
399  {
400  const Number exact_val = slitfunc(context, xyz[qp]);
401 
402  const Number discrete_val = context.interior_value(0, qp);
403 
404  LIBMESH_ASSERT_FP_EQUAL(libmesh_real(exact_val),
405  libmesh_real(discrete_val),
406  TOLERANCE*TOLERANCE);
407  }
408  }
409  }
410 };
411 
412 CPPUNIT_TEST_SUITE_REGISTRATION( SlitMeshTest );
413 CPPUNIT_TEST_SUITE_REGISTRATION( SlitMeshRefinedMeshTest );
414 CPPUNIT_TEST_SUITE_REGISTRATION( SlitMeshRefinedSystemTest );
libMesh::System
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:100
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:67
libMesh::Number
Real Number
Definition: libmesh_common.h:195
libMesh::System::n_vars
unsigned int n_vars() const
Definition: system.h:2155
SlitFunc::clone
virtual std::unique_ptr< FEMFunctionBase< Number > > clone() const override
Definition: slit_mesh_test.C:28
SlitMeshRefinedSystemTest
Definition: slit_mesh_test.C:236
libMesh::Mesh
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
libMesh::FEAbstract::get_xyz
const std::vector< Point > & get_xyz() const
Definition: fe_abstract.h:237
libMesh::EquationSystems::add_system
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.
Definition: equation_systems.C:345
libMesh::DistributedMesh::active_local_element_ptr_range
virtual SimpleRange< element_iterator > active_local_element_ptr_range() override
Definition: distributed_mesh.h:372
libMesh::EquationSystems::WRITE_DATA
Definition: equation_systems.h:93
libMesh::FEMFunctionBase
FEMFunctionBase is a base class from which users can derive in order to define "function-like" object...
Definition: dirichlet_boundaries.h:43
libMesh::libmesh_real
T libmesh_real(T a)
Definition: libmesh_common.h:166
libMesh::DofObject::set_id
dof_id_type & set_id()
Definition: dof_object.h:776
libMesh::DistributedMesh::add_point
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) override
functions for adding /deleting nodes elements.
Definition: distributed_mesh.C:625
libMesh::MeshBase::elem_ref
virtual const Elem & elem_ref(const dof_id_type i) const
Definition: mesh_base.h:521
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
Definition: factoryfunction.C:55
SlitMeshRefinedSystemTest::tearDown
void tearDown()
Definition: slit_mesh_test.C:290
libMesh::FEGenericBase
This class forms the foundation from which generic finite elements may be derived.
Definition: exact_error_estimator.h:39
libMesh::EquationSystems::get_system
const T_sys & get_system(const std::string &name) const
Definition: equation_systems.h:757
libMesh::FEMContext::get_elem
const Elem & get_elem() const
Accessor for current Elem object.
Definition: fem_context.h:896
libMesh::TOLERANCE
static const Real TOLERANCE
Definition: libmesh_common.h:128
SlitMeshTest::tearDown
void tearDown()
Definition: slit_mesh_test.C:161
SlitMeshRefinedMeshTest::setUp
void setUp()
Definition: slit_mesh_test.C:212
SlitMeshRefinedSystemTest::setUp
void setUp()
Definition: slit_mesh_test.C:264
SlitMeshTest
Definition: slit_mesh_test.C:61
dim
unsigned int dim
Definition: adaptivity_ex3.C:113
libMesh::Elem::centroid
virtual Point centroid() const
Definition: elem.C:345
libMesh::TypeVector::size
auto size() const -> decltype(std::norm(T()))
Definition: type_vector.h:944
libMesh::MeshRefinement
Implements (adaptive) mesh refinement algorithms for a MeshBase.
Definition: mesh_refinement.h:61
libMesh::DistributedMesh::n_nodes
virtual dof_id_type n_nodes() const override
Definition: distributed_mesh.h:229
std::abs
MetaPhysicL::DualNumber< T, D > abs(const MetaPhysicL::DualNumber< T, D > &in)
SlitMeshRefinedSystemTest::_sys
System * _sys
Definition: slit_mesh_test.C:259
SlitMeshRefinedSystemTest::testSystem
void testSystem()
Definition: slit_mesh_test.C:309
SlitMeshRefinedSystemTest::testRestart
void testRestart()
Definition: slit_mesh_test.C:343
libMesh::DistributedMesh::n_elem
virtual dof_id_type n_elem() const override
Definition: distributed_mesh.h:232
SlitMeshTest::setUp
void setUp()
Definition: slit_mesh_test.C:154
SlitMeshRefinedSystemTest::_es
EquationSystems * _es
Definition: slit_mesh_test.C:260
SlitMeshTest::testMesh
void testMesh()
Definition: slit_mesh_test.C:166
libMesh::System::add_variable
unsigned int add_variable(const std::string &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:1069
TestCommWorld
libMesh::Parallel::Communicator * TestCommWorld
Definition: driver.C:111
libMesh::EquationSystems::init
virtual void init()
Initialize all the systems.
Definition: equation_systems.C:96
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
libMesh::UnstructuredMesh::read
virtual void read(const std::string &name, void *mesh_data=nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false) override
Reads the file specified by name.
Definition: unstructured_mesh.C:620
SlitFunc::SlitFunc
SlitFunc()
Definition: slit_mesh_test.C:21
libMesh::DistributedMesh::add_elem
virtual Elem * add_elem(Elem *e) override
Add elem e to the end of the element array.
Definition: distributed_mesh.C:436
libMesh::Node
A Node is like a Point, but with more information.
Definition: node.h:52
SlitMeshTest::build_mesh
void build_mesh()
Definition: slit_mesh_test.C:79
libMesh::Quad4
The QUAD4 is an element in 2D composed of 4 nodes.
Definition: face_quad4.h:51
libMesh::DofObject::unique_id
unique_id_type unique_id() const
Definition: dof_object.h:784
libMesh::DenseVector::size
virtual unsigned int size() const override
Definition: dense_vector.h:92
SlitMeshTest::_mesh
Mesh * _mesh
Definition: slit_mesh_test.C:77
CPPUNIT_TEST_SUITE_REGISTRATION
CPPUNIT_TEST_SUITE_REGISTRATION(SlitMeshTest)
libMesh::EquationSystems::read
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.
Definition: equation_systems_io.C:90
libMesh::EquationSystems
This is the EquationSystems class.
Definition: equation_systems.h:74
libMesh::Elem::set_node
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:2059
SlitFunc::~SlitFunc
~SlitFunc()
Definition: slit_mesh_test.C:23
libMesh::EquationSystems::WRITE_SERIAL_FILES
Definition: equation_systems.h:96
libMesh::EquationSystems::reinit
virtual void reinit()
Handle any mesh changes and reinitialize all the systems on the updated mesh.
Definition: equation_systems.C:121
libMesh::DistributedMesh::query_elem_ptr
virtual const Elem * query_elem_ptr(const dof_id_type i) const override
Definition: distributed_mesh.C:404
SlitMeshRefinedMeshTest::testMesh
void testMesh()
Definition: slit_mesh_test.C:223
libmesh_cppunit.h
SlitMeshRefinedSystemTest::testMesh
void testMesh()
Definition: slit_mesh_test.C:297
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:100
libMesh::EquationSystems::write
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.
Definition: equation_systems_io.C:378
libMesh::FEGenericBase::get_phi
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:206
libMesh::MeshBase::allow_renumbering
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:1025
libMesh::Elem::node_id
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:1977
libMesh::Elem::node_ref
const Node & node_ref(const unsigned int i) const
Definition: elem.h:2031
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:121
test_comm.h
libMesh::MeshBase::prepare_for_use
void prepare_for_use(const bool skip_renumber_nodes_and_elements=false, const bool skip_find_neighbors=false)
Prepare a newly ecreated (or read) mesh for use.
Definition: mesh_base.C:318
libMesh::System::project_solution
void project_solution(FunctionBase< Number > *f, FunctionBase< Gradient > *g=nullptr) const
Projects arbitrary functions onto the current solution.
Definition: system_projection.C:950
libMesh::MeshBase::set_mesh_dimension
void set_mesh_dimension(unsigned char d)
Resets the logical dimension of the mesh.
Definition: mesh_base.h:218
SlitFunc::init_context
virtual void init_context(const FEMContext &) override
Prepares a context object for use.
Definition: slit_mesh_test.C:25
libMesh::MeshRefinement::uniformly_refine
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.
Definition: mesh_refinement.C:1678
libMesh::FIRST
Definition: enum_order.h:42
SlitFunc
Definition: slit_mesh_test.C:17
SlitMeshRefinedMeshTest
Definition: slit_mesh_test.C:193
libMesh::DenseVector< Number >
libMesh::DistributedMesh::node_ptr
virtual const Node * node_ptr(const dof_id_type i) const override
Definition: distributed_mesh.C:328
libMesh::FEMContext
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:62
Generated on Sat Jan 25 2020 12:07:06 for libMesh by   doxygen 1.8.16