libMesh
Public Member Functions | Protected Attributes | Private Member Functions | List of all members
ExodusTest< elem_type > Class Template Reference
Inheritance diagram for ExodusTest< elem_type >:
[legend]

Public Member Functions

void test_read_gold ()
 
void test_write ()
 
void setUp ()
 

Protected Attributes

std::unique_ptr< Mesh_mesh
 
std::string libmesh_suite_name
 

Private Member Functions

bool meshes_equal_enough (Mesh &other_mesh)
 

Detailed Description

template<ElemType elem_type>
class ExodusTest< elem_type >

Definition at line 13 of file exodus_test.C.

Member Function Documentation

◆ meshes_equal_enough()

template<ElemType elem_type>
bool ExodusTest< elem_type >::meshes_equal_enough ( Mesh other_mesh)
inlineprivate

Definition at line 17 of file exodus_test.C.

References libMesh::MeshBase::elem_ref(), libMesh::DofObject::id(), libMesh::make_range(), libMesh::DistributedMesh::max_elem_id(), libMesh::DistributedMesh::max_node_id(), libMesh::DistributedMesh::node_ptr(), libMesh::Elem::node_ref_range(), libMesh::TypeVector< T >::norm(), libMesh::DofObject::processor_id(), libMesh::DistributedMesh::renumber_elem(), libMesh::DistributedMesh::renumber_node(), libMesh::DistributedMesh::set_next_unique_id(), libMesh::DofObject::set_unique_id(), and libMesh::MeshBase::sub_point_locator().

18  {
19  // We'll need to fix up processor_id() and unique_id() values
20  // before we can operator== these meshes. But worse: our gold
21  // meshes might have been numbered differently to our generated
22  // meshes. Some of our generated mesh options practically
23  // *require* renumbering (e.g. after interior HEX20 nodes are
24  // deleted, ExodusII still wants to see a contiguous numbering),
25  // but ReplicatedMesh and DistributedMesh renumber differently.
26  //
27  // So, let's renumber too.
28 
29  MeshSerializer serialthis(*this->_mesh);
30  MeshSerializer serialother(other_mesh);
31 
32  const dof_id_type max_elem_id = this->_mesh->max_elem_id();
33  const dof_id_type max_node_id = this->_mesh->max_node_id();
34 
35  CPPUNIT_ASSERT_EQUAL(max_elem_id, other_mesh.max_elem_id());
36  CPPUNIT_ASSERT_EQUAL(max_node_id, other_mesh.max_node_id());
37 
38  auto locator = other_mesh.sub_point_locator();
39 
40  for (Elem * e1 : this->_mesh->element_ptr_range())
41  {
42  const Elem * e2c = (*locator)(e1->vertex_average());
43  CPPUNIT_ASSERT(e2c);
44  Elem & e2 = other_mesh.elem_ref(e2c->id());
45  e1->processor_id() = 0;
46  e2.processor_id() = 0;
47 
48  const dof_id_type e1_id = e1->id();
49  const dof_id_type e2_id = e2.id();
50  // Do a swap if necessary, using a free temporary id
51  if (e1_id != e2_id)
52  {
53  other_mesh.renumber_elem(e1_id, max_elem_id);
54  other_mesh.renumber_elem(e2_id, e1_id);
55  other_mesh.renumber_elem(max_elem_id, e2_id);
56  }
57 
58 #ifdef LIBMESH_ENABLE_UNIQUE_ID
59  e2.set_unique_id(e1->unique_id());
60 #endif
61  }
62 
63  for (Node * n1 : this->_mesh->node_ptr_range())
64  {
65  const Elem * e1c = (*locator)(*n1);
66  Node * n2 = nullptr;
67  for (const Node & n : e1c->node_ref_range())
68  {
69 #if LIBMESH_DEFAULT_QUADRUPLE_PRECISION
70  const Point diff = Point(*n1)-Point(n);
71 
72  // We're testing against ExodusII input, and if we're in
73  // triple or quadruple precision that means our lovely
74  // higher-precision node coordinates got truncated to double
75  // to be written. We need to adjust ours or they won't
76  // satisfy operator== later.
77 
78  // We're *also* testing against gold files that were
79  // calculated at double precision, so just casting a higher
80  // precision calculation to double won't give the exact same
81  // result, we have to account for error.
82  if (diff.norm() < 1e-15)
83  for (auto d : make_range(LIBMESH_DIM))
84  (*n1)(d) = double(n(d));
85 #endif
86  if (Point(*n1) == Point(n))
87  n2 = other_mesh.node_ptr(n.id());
88  }
89  CPPUNIT_ASSERT(n2);
90  n1->processor_id() = 0;
91  n2->processor_id() = 0;
92 
93  const dof_id_type n1_id = n1->id();
94  const dof_id_type n2_id = n2->id();
95  // Do a swap if necessary, using a free temporary id
96  if (n1_id != n2_id)
97  {
98  other_mesh.renumber_node(n1_id,max_node_id);
99  other_mesh.renumber_node(n2_id,n1_id);
100  other_mesh.renumber_node(max_node_id, n2_id);
101  }
102 
103 #ifdef LIBMESH_ENABLE_UNIQUE_ID
104  n2->set_unique_id(n1->unique_id());
105 #endif
106  }
107 
108 #ifdef LIBMESH_ENABLE_UNIQUE_ID
109  other_mesh.set_next_unique_id(this->_mesh->parallel_max_unique_id());
110  this->_mesh->set_next_unique_id(this->_mesh->parallel_max_unique_id());
111 #endif
112 
113  return *this->_mesh == other_mesh;
114  }
unique_id_type & set_unique_id()
Definition: dof_object.h:858
A Node is like a Point, but with more information.
Definition: node.h:52
auto norm() const -> decltype(std::norm(T()))
Definition: type_vector.h:907
std::unique_ptr< PointLocatorBase > sub_point_locator() const
Definition: mesh_base.C:1638
virtual void renumber_elem(dof_id_type old_id, dof_id_type new_id) override final
Changes the id of element old_id, both by changing elem(old_id)->id() and by moving elem(old_id) in t...
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:22
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
dof_id_type id() const
Definition: dof_object.h:828
virtual void set_next_unique_id(unique_id_type id) override
Sets the next available unique id to be used.
SimpleRange< NodeRefIter > node_ref_range()
Returns a range with all nodes of an element, usable in range-based for loops.
Definition: elem.h:2665
virtual void renumber_node(dof_id_type old_id, dof_id_type new_id) override final
Changes the id of node old_id, both by changing node(old_id)->id() and by moving node(old_id) in the ...
Temporarily serialize a DistributedMesh for non-distributed-mesh capable code paths.
virtual dof_id_type max_node_id() const override final
virtual const Elem & elem_ref(const dof_id_type i) const
Definition: mesh_base.h:639
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
virtual dof_id_type max_elem_id() const override final
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
uint8_t dof_id_type
Definition: id_types.h:67
virtual const Node * node_ptr(const dof_id_type i) const override final

◆ setUp()

template<ElemType elem_type>
void PerElemTest< elem_type >::setUp ( )
inlineinherited

Definition at line 26 of file elem_test.h.

References libMesh::Elem::build(), libMesh::MeshTools::Generation::build_cube(), libMesh::C0POLYGON, libMesh::C0POLYHEDRON, dim, libMesh::index_range(), libMesh::INFHEX16, libMesh::INFHEX18, libMesh::INFHEX8, libMesh::INFPRISM12, libMesh::INFPRISM6, libMesh::INFQUAD4, libMesh::INFQUAD6, libMesh::make_range(), libMesh::Real, and libMesh::Elem::set_node().

27  {
28  const Real minpos = 1.5, maxpos = 5.5;
29  const unsigned int N = 2;
30 
31  _mesh = std::make_unique<Mesh>(*TestCommWorld);
32 
33  std::unique_ptr<Elem> test_elem;
34 
35  if (elem_type != C0POLYHEDRON)
36  test_elem = Elem::build(elem_type);
37 
38 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
39 #if LIBMESH_DIM > 1
40  if (test_elem.get() && test_elem->infinite())
41  {
42  Elem * elem = _mesh->add_elem(std::move(test_elem));
43 
44  const auto add_point =
45  [this, elem](const unsigned int i,
46  const Real x,
47  const Real y,
48  const Real
49 #if LIBMESH_DIM == 3
50  z
51 #endif
52  )
53  {
54 #if LIBMESH_DIM == 2
55  auto node = _mesh->add_point(Point(x, y), i);
56 #else
57  auto node = _mesh->add_point(Point(x, y, z), i);
58 #endif
59  elem->set_node(i, node);
60  };
61 
62  const Real halfpos = (minpos + maxpos) / 2.;
63 
64  if (elem_type == INFQUAD4 || elem_type == INFQUAD6 ||
65  elem_type == INFHEX8 || elem_type == INFHEX16 || elem_type == INFHEX18)
66  {
67  const bool is_quad = (elem_type == INFQUAD4 || elem_type == INFQUAD6);
68 
69  add_point(0, minpos, minpos, minpos);
70  add_point(1, maxpos, minpos, minpos);
71  add_point(2+is_quad, maxpos, maxpos, minpos);
72  add_point(3-is_quad, minpos, maxpos, minpos);
73 
74  if (elem_type == INFQUAD6)
75  {
76  add_point(4, halfpos, minpos, minpos);
77  add_point(5, halfpos, maxpos, minpos);
78  }
79  }
80  if (elem_type == INFHEX8 || elem_type == INFHEX16 || elem_type == INFHEX18)
81  {
82  add_point(4, minpos, minpos, maxpos);
83  add_point(5, maxpos, minpos, maxpos);
84  add_point(6, maxpos, maxpos, maxpos);
85  add_point(7, minpos, maxpos, maxpos);
86 
87  if (elem_type == INFHEX16 || elem_type == INFHEX18)
88  {
89  add_point(8, halfpos, minpos, minpos);
90  add_point(9, maxpos, halfpos, minpos);
91  add_point(10, halfpos, maxpos, minpos);
92  add_point(11, minpos, halfpos, minpos);
93  add_point(12, halfpos, minpos, maxpos);
94  add_point(13, maxpos, halfpos, maxpos);
95  add_point(14, halfpos, maxpos, maxpos);
96  add_point(15, minpos, halfpos, maxpos);
97  }
98  if (elem_type == INFHEX18)
99  {
100  add_point(16, halfpos, halfpos, minpos);
101  add_point(17, halfpos, halfpos, maxpos);
102  }
103  }
104  if (elem_type == INFPRISM6 || elem_type == INFPRISM12)
105  {
106  add_point(0, minpos, minpos, minpos);
107  add_point(1, maxpos, minpos, minpos);
108  add_point(2, halfpos, maxpos, minpos);
109  add_point(3, minpos, minpos, maxpos);
110  add_point(4, maxpos, minpos, maxpos);
111  add_point(5, halfpos, maxpos, maxpos);
112 
113  if (elem_type == INFPRISM12)
114  {
115  add_point(6, halfpos, minpos, minpos);
116  add_point(7, (halfpos + maxpos) / 2., halfpos, minpos);
117  add_point(8, (halfpos + minpos) / 2., halfpos, minpos);
118  add_point(9, halfpos, minpos, maxpos);
119  add_point(10, (halfpos + maxpos) / 2., halfpos, maxpos);
120  add_point(11, (halfpos + minpos) / 2., halfpos, maxpos);
121  }
122  }
123 
124  _mesh->prepare_for_use();
125  }
126  else
127 #endif // LIBMESH_DIM > 1
128 #endif // LIBMESH_ENABLE_INFINITE_ELEMENTS
129  if (elem_type == C0POLYGON)
130  {
131  // We're not going to implement build_square for e.g.
132  // pentagons any time soon.
133  //
134  // We should probably implement build_dual_mesh, though, and
135  // run it on a perturbed or triangle input so we don't just
136  // get quads in the output...
137 
138  _mesh->add_point(Point(0, 0), 0);
139  _mesh->add_point(Point(1, 0), 1);
140  _mesh->add_point(Point(1.5, 0.5), 2);
141  _mesh->add_point(Point(1, 1), 3);
142  _mesh->add_point(Point(0, 1), 4);
143 
144  std::unique_ptr<Elem> polygon = std::make_unique<C0Polygon>(5);
145  for (auto i : make_range(5))
146  polygon->set_node(i, _mesh->node_ptr(i));
147  polygon->set_id() = 0;
148 
149  _mesh->add_elem(std::move(polygon));
150  _mesh->prepare_for_use();
151  }
152  else if (elem_type == C0POLYHEDRON)
153  {
154  // There's even less point in having a build_cube for general
155  // polyhedra, so again we'll hand-make one for testing and
156  // we'll plan on making a build_dual_mesh for the future.
157 
158 #if 0
159  // Try a truncated cube for relatively simple verification.
160 
161  // We have some differing orientations on the top sides, to
162  // test handling of that.
163  // Lower octagon points, counterclockwise
164  _mesh->add_point(Point(1/Real(3), 0, 0), 0);
165  _mesh->add_point(Point(2/Real(3), 0, 0), 1);
166  _mesh->add_point(Point(1, 1/Real(3), 0), 2);
167  _mesh->add_point(Point(1, 2/Real(3), 0), 3);
168  _mesh->add_point(Point(2/Real(3), 1, 0), 4);
169  _mesh->add_point(Point(1/Real(3), 1, 0), 5);
170  _mesh->add_point(Point(0, 2/Real(3), 0), 6);
171  _mesh->add_point(Point(0, 1/Real(3), 0), 7);
172 
173  // Lower-middle points, counterclockwise
174  _mesh->add_point(Point(0, 0, 1/Real(3)), 8);
175  _mesh->add_point(Point(1, 0, 1/Real(3)), 9);
176  _mesh->add_point(Point(1, 1, 1/Real(3)), 10);
177  _mesh->add_point(Point(0, 1, 1/Real(3)), 11);
178 
179  // Upper-middle points, counterclockwise
180  _mesh->add_point(Point(0, 0, 2/Real(3)), 12);
181  _mesh->add_point(Point(1, 0, 2/Real(3)), 13);
182  _mesh->add_point(Point(1, 1, 2/Real(3)), 14);
183  _mesh->add_point(Point(0, 1, 2/Real(3)), 15);
184 
185  // Upper octagon points, counterclockwise
186  _mesh->add_point(Point(1/Real(3), 0, 1), 16);
187  _mesh->add_point(Point(2/Real(3), 0, 1), 17);
188  _mesh->add_point(Point(1, 1/Real(3), 1), 18);
189  _mesh->add_point(Point(1, 2/Real(3), 1), 19);
190  _mesh->add_point(Point(2/Real(3), 1, 1), 20);
191  _mesh->add_point(Point(1/Real(3), 1, 1), 21);
192  _mesh->add_point(Point(0, 2/Real(3), 1), 22);
193  _mesh->add_point(Point(0, 1/Real(3), 1), 23);
194 
195  const std::vector<std::vector<unsigned int>> nodes_on_side =
196  { {0, 1, 2, 3, 4, 5, 6, 7}, // min z
197  {0, 1, 9, 13, 17, 16, 12, 8}, // min y
198  {2, 3, 10, 14, 19, 18, 13, 9}, // max x
199  {4, 5, 11, 15, 21, 20, 14, 10}, // max y
200  {6, 7, 8, 12, 23, 22, 15, 11}, // min x
201  {16, 17, 18, 19, 20, 21, 22, 23}, // max z
202  {7, 0, 8}, // max nothing
203  {1, 2, 9}, // max x
204  {3, 4, 10}, // max xy
205  {5, 6, 11}, // max y
206  {23, 16, 12}, // max z
207  {17, 18, 13}, // max xz
208  {19, 20, 14}, // max xyz
209  {21, 22, 15} }; // max yz
210 #endif
211 
212 #if 1
213  // Or try a plain cube, for simpler debugging
214  _mesh->add_point(Point(0, 0, 0), 0);
215  _mesh->add_point(Point(1, 0, 0), 1);
216  _mesh->add_point(Point(1, 1, 0), 2);
217  _mesh->add_point(Point(0, 1, 0), 3);
218  _mesh->add_point(Point(0, 0, 1), 4);
219  _mesh->add_point(Point(1, 0, 1), 5);
220  _mesh->add_point(Point(1, 1, 1), 6);
221  _mesh->add_point(Point(0, 1, 1), 7);
222 
223  // With some combinations of face triangulations, even a
224  // simple cube has no tetrahedralization that doesn't have
225  // either interior discontinuities or a 0-volume pancake tet!
226  //
227  // The initial "natural" way to orient our sides is commented
228  // out here; permutations that fix diagonals for us are used
229  // instead.
230  const std::vector<std::vector<unsigned int>> nodes_on_side =
231  { {0, 1, 2, 3}, // min z
232  {0, 1, 5, 4}, // min y
233  // {1, 2, 6, 5}, // max x - bad
234  {2, 6, 5, 1}, // max x
235  {2, 3, 7, 6}, // max y
236  // {3, 0, 4, 7}, // min x - bad
237  {0, 4, 7, 3}, // min x
238  // {4, 5, 6, 7} }; // max z - bad
239  {5, 6, 7, 4} }; // max z
240 #endif
241 
242  // Build all the sides.
243  std::vector<std::shared_ptr<Polygon>> sides(nodes_on_side.size());
244 
245  for (auto s : index_range(nodes_on_side))
246  {
247  const auto & nodes_on_s = nodes_on_side[s];
248  sides[s] = std::make_shared<C0Polygon>(nodes_on_s.size());
249  for (auto i : index_range(nodes_on_s))
250  sides[s]->set_node(i, _mesh->node_ptr(nodes_on_s[i]));
251  }
252 
253  std::unique_ptr<Elem> polyhedron = std::make_unique<C0Polyhedron>(sides);
254  _mesh->add_elem(std::move(polyhedron));
255  _mesh->prepare_for_use();
256  }
257  else
258  {
259  const unsigned int dim = test_elem->dim();
260  const unsigned int use_x = dim > 0;
261  const unsigned int use_y = dim > 1;
262  const unsigned int use_z = dim > 2;
263 
265  N*use_x, N*use_y, N*use_z,
266  minpos, maxpos,
267  minpos, use_y*maxpos,
268  minpos, use_z*maxpos,
269  elem_type);
270  }
271 
272  // Use non-default subdomain ids so we can properly test their
273  // preservation
274  for (const auto & elem :
275  this->_mesh->element_ptr_range())
276  elem->subdomain_id() = 10 + (elem->id() % 10);
277 
278  // And make sure our mesh's cache knows about them all for later
279  // tests comparisons
280  this->_mesh->cache_elem_data();
281  }
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:2558
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:22
unsigned int dim
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
dof_id_type id() const
Definition: dof_object.h:828
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
subdomain_id_type subdomain_id() const
Definition: elem.h:2582
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117
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.

◆ test_read_gold()

template<ElemType elem_type>
void ExodusTest< elem_type >::test_read_gold ( )
inline

Definition at line 118 of file exodus_test.C.

References libMesh::MeshCommunication::broadcast(), libMesh::Utility::enum_to_string(), libMesh::MeshBase::prepare_for_use(), libMesh::ParallelObject::processor_id(), libMesh::ExodusII_IO::read(), and TestCommWorld.

119  {
120  LOG_UNIT_TEST;
121 
122  Mesh input_mesh(*TestCommWorld);
123 
124  ExodusII_IO exii(input_mesh);
125  if (input_mesh.processor_id() == 0)
126  exii.read("meshes/exodus_elements/read_exodus_" +
127  Utility::enum_to_string(elem_type) + ".e");
128 
129  MeshCommunication().broadcast(input_mesh);
130  input_mesh.prepare_for_use();
131 
132  CPPUNIT_ASSERT(this->meshes_equal_enough(input_mesh));
133  }
libMesh::Parallel::Communicator * TestCommWorld
Definition: driver.C:171
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
bool meshes_equal_enough(Mesh &other_mesh)
Definition: exodus_test.C:17
This is the MeshCommunication class.
std::string enum_to_string(const T e)
void broadcast(MeshBase &) const
Finds all the processors that may contain elements that neighbor my elements.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50

◆ test_write()

template<ElemType elem_type>
void ExodusTest< elem_type >::test_write ( )
inline

Definition at line 135 of file exodus_test.C.

References libMesh::MeshCommunication::broadcast(), libMesh::Utility::enum_to_string(), libMesh::MeshBase::prepare_for_use(), libMesh::ParallelObject::processor_id(), libMesh::ExodusII_IO::read(), TestCommWorld, and libMesh::ExodusII_IO::write().

136  {
137  LOG_UNIT_TEST;
138 
139  // This is a *buffered* write; we use scope to make sure the
140  // ExodusII_IO object gets destructed (and thus is guaranteed to
141  // finish writing and close the file) before we try to read what
142  // was written.
143  {
144  ExodusII_IO exii(*this->_mesh);
145  exii.write("write_exodus_" +
146  Utility::enum_to_string(elem_type) + ".e");
147  }
148 
149  Mesh input_mesh(*TestCommWorld);
150  ExodusII_IO exii_input(input_mesh);
151  if (input_mesh.processor_id() == 0)
152  exii_input.read("write_exodus_" +
153  Utility::enum_to_string(elem_type) + ".e");
154 
155  MeshCommunication().broadcast(input_mesh);
156  input_mesh.prepare_for_use();
157 
158  CPPUNIT_ASSERT(this->meshes_equal_enough(input_mesh));
159  }
libMesh::Parallel::Communicator * TestCommWorld
Definition: driver.C:171
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:22
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
bool meshes_equal_enough(Mesh &other_mesh)
Definition: exodus_test.C:17
This is the MeshCommunication class.
std::string enum_to_string(const T e)
void broadcast(MeshBase &) const
Finds all the processors that may contain elements that neighbor my elements.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50

Member Data Documentation

◆ _mesh

template<ElemType elem_type>
std::unique_ptr<Mesh> PerElemTest< elem_type >::_mesh
protectedinherited

Definition at line 22 of file elem_test.h.

◆ libmesh_suite_name

template<ElemType elem_type>
std::string PerElemTest< elem_type >::libmesh_suite_name
protectedinherited

Definition at line 23 of file elem_test.h.


The documentation for this class was generated from the following file: