1 #include <libmesh/equation_systems.h>
2 #include <libmesh/linear_implicit_system.h>
3 #include <libmesh/mesh.h>
4 #include <libmesh/mesh_communication.h>
5 #include <libmesh/mesh_generation.h>
6 #include <libmesh/numeric_vector.h>
7 #include <libmesh/replicated_mesh.h>
8 #include <libmesh/dyna_io.h>
9 #include <libmesh/exodusII_io.h>
10 #include <libmesh/dof_map.h>
11
12 #include "test_comm.h"
13 #include "libmesh_cppunit.h"
14
15
16 using namespace libMesh;
17
18
x_plus_y(const Point & p,const Parameters &,const std::string &,const std::string &)19 Number x_plus_y (const Point& p,
20 const Parameters&,
21 const std::string&,
22 const std::string&)
23 {
24 const Real & x = p(0);
25 const Real & y = p(1);
26
27 return x + y;
28 }
29
30
31 class MeshInputTest : public CppUnit::TestCase {
32 public:
33 CPPUNIT_TEST_SUITE( MeshInputTest );
34
35 #if LIBMESH_DIM > 1
36 #ifdef LIBMESH_HAVE_EXODUS_API
37 CPPUNIT_TEST( testExodusCopyElementSolution );
38 CPPUNIT_TEST( testExodusReadHeader );
39 #ifndef LIBMESH_USE_COMPLEX_NUMBERS
40 // Eventually this will support complex numbers.
41 CPPUNIT_TEST( testExodusWriteElementDataFromDiscontinuousNodalData );
42 #endif // LIBMESH_USE_COMPLEX_NUMBERS
43 #endif // LIBMESH_HAVE_EXODUS_API
44 CPPUNIT_TEST( testDynaReadElem );
45 CPPUNIT_TEST( testDynaReadPatch );
46 #endif // LIBMESH_DIM > 1
47
48 CPPUNIT_TEST_SUITE_END();
49
50 private:
51
52 public:
setUp()53 void setUp()
54 {}
55
tearDown()56 void tearDown()
57 {}
58
59 #ifdef LIBMESH_HAVE_EXODUS_API
testExodusReadHeader()60 void testExodusReadHeader ()
61 {
62 // first scope: write file
63 {
64 ReplicatedMesh mesh(*TestCommWorld);
65 MeshTools::Generation::build_square (mesh, 3, 3, 0., 1., 0., 1.);
66 ExodusII_IO exii(mesh);
67 mesh.write("read_header_test.e");
68 }
69
70 // Make sure that the writing is done before the reading starts.
71 TestCommWorld->barrier();
72
73 // second scope: read header
74 // Note: The header information is read from file on processor 0
75 // and then broadcast to the other procs, so with this test we are
76 // checking both that the header information is read correctly and
77 // that it is correctly communicated to other procs.
78 {
79 ReplicatedMesh mesh(*TestCommWorld);
80 ExodusII_IO exii(mesh);
81 ExodusHeaderInfo header_info = exii.read_header("read_header_test.e");
82
83 // Make sure the header information is as expected.
84 CPPUNIT_ASSERT_EQUAL(std::string(header_info.title.data()), std::string("read_header_test.e"));
85 CPPUNIT_ASSERT_EQUAL(header_info.num_dim, 2);
86 CPPUNIT_ASSERT_EQUAL(header_info.num_elem, 9);
87 CPPUNIT_ASSERT_EQUAL(header_info.num_elem_blk, 1);
88 CPPUNIT_ASSERT_EQUAL(header_info.num_node_sets, 0);
89 CPPUNIT_ASSERT_EQUAL(header_info.num_side_sets, 4);
90 CPPUNIT_ASSERT_EQUAL(header_info.num_edge_blk, 0);
91 CPPUNIT_ASSERT_EQUAL(header_info.num_edge, 0);
92 }
93 }
94
testExodusCopyElementSolution()95 void testExodusCopyElementSolution ()
96 {
97 {
98 Mesh mesh(*TestCommWorld);
99
100 EquationSystems es(mesh);
101 System &sys = es.add_system<System> ("SimpleSystem");
102 sys.add_variable("e", CONSTANT, MONOMIAL);
103
104 MeshTools::Generation::build_square (mesh,
105 3, 3,
106 0., 1., 0., 1.);
107
108 es.init();
109 sys.project_solution(x_plus_y, nullptr, es.parameters);
110
111 ExodusII_IO exii(mesh);
112
113 // Don't try to write element data as nodal data
114 std::set<std::string> sys_list;
115 exii.write_equation_systems("mesh_with_soln.e", es, &sys_list);
116
117 // Just write it as element data
118 exii.write_element_data(es);
119 }
120
121 // copy_elemental_solution currently requires ReplicatedMesh
122 {
123 ReplicatedMesh mesh(*TestCommWorld);
124
125 EquationSystems es(mesh);
126 System &sys = es.add_system<System> ("SimpleSystem");
127 sys.add_variable("teste", CONSTANT, MONOMIAL);
128
129 ExodusII_IO exii(mesh);
130
131 if (mesh.processor_id() == 0)
132 exii.read("mesh_with_soln.e");
133 MeshCommunication().broadcast(mesh);
134 mesh.prepare_for_use();
135
136 es.init();
137
138 // Read the solution e into variable teste.
139 //
140 // With complex numbers, we'll only bother reading the real
141 // part.
142 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
143 exii.copy_elemental_solution(sys, "teste", "r_e");
144 #else
145 exii.copy_elemental_solution(sys, "teste", "e");
146 #endif
147
148 // Exodus only handles double precision
149 Real exotol = std::max(TOLERANCE*TOLERANCE, Real(1e-12));
150
151 for (Real x = Real(1.L/6.L); x < 1; x += Real(1.L/3.L))
152 for (Real y = Real(1.L/6.L); y < 1; y += Real(1.L/3.L))
153 {
154 Point p(x,y);
155 LIBMESH_ASSERT_FP_EQUAL(libmesh_real(sys.point_value(0,p)),
156 libmesh_real(x+y),
157 exotol);
158 }
159 }
160 }
161
162
163 #ifndef LIBMESH_USE_COMPLEX_NUMBERS
testExodusWriteElementDataFromDiscontinuousNodalData()164 void testExodusWriteElementDataFromDiscontinuousNodalData()
165 {
166 // first scope: write file
167 {
168 Mesh mesh(*TestCommWorld);
169
170 EquationSystems es(mesh);
171 System & sys = es.add_system<System> ("SimpleSystem");
172 sys.add_variable("u", FIRST, L2_LAGRANGE);
173
174 MeshTools::Generation::build_cube
175 (mesh, 2, 2, 2, 0., 1., 0., 1., 0., 1., HEX8);
176
177 es.init();
178
179 // Set solution u^e_i = i, for the ith vertex of a given element e.
180 const DofMap & dof_map = sys.get_dof_map();
181 std::vector<dof_id_type> dof_indices;
182 for (const auto & elem : mesh.element_ptr_range())
183 {
184 dof_map.dof_indices(elem, dof_indices, /*var_id=*/0);
185 for (unsigned int i=0; i<dof_indices.size(); ++i)
186 sys.solution->set(dof_indices[i], i);
187 }
188 sys.solution->close();
189
190 // Now write to file.
191 ExodusII_IO exii(mesh);
192
193 // Don't try to write element data as averaged nodal data.
194 std::set<std::string> sys_list;
195 exii.write_equation_systems("elemental_from_nodal.e", es, &sys_list);
196
197 // Write one elemental data field per vertex value.
198 sys_list = {"SimpleSystem"};
199
200 exii.write_element_data_from_discontinuous_nodal_data
201 (es, &sys_list, /*var_suffix=*/"_elem_corner_");
202 } // end first scope
203
204 // second scope: read values back in, verify they are correct.
205 {
206 std::vector<std::string> file_var_names =
207 {"u_elem_corner_0",
208 "u_elem_corner_1",
209 "u_elem_corner_2",
210 "u_elem_corner_3"};
211 std::vector<Real> expected_values = {0., 1., 2., 3.};
212
213 // copy_elemental_solution currently requires ReplicatedMesh
214 ReplicatedMesh mesh(*TestCommWorld);
215
216 EquationSystems es(mesh);
217 System & sys = es.add_system<System> ("SimpleSystem");
218 for (auto i : index_range(file_var_names))
219 sys.add_variable(file_var_names[i], CONSTANT, MONOMIAL);
220
221 ExodusII_IO exii(mesh);
222
223 if (mesh.processor_id() == 0)
224 exii.read("elemental_from_nodal.e");
225 MeshCommunication().broadcast(mesh);
226 mesh.prepare_for_use();
227
228 es.init();
229
230 for (auto i : index_range(file_var_names))
231 exii.copy_elemental_solution
232 (sys, sys.variable_name(i), file_var_names[i]);
233
234 // Check that the values we read back in are as expected.
235 for (const auto & elem : mesh.active_element_ptr_range())
236 for (auto i : index_range(file_var_names))
237 {
238 Real read_val = sys.point_value(i, elem->centroid());
239 LIBMESH_ASSERT_FP_EQUAL
240 (expected_values[i], read_val, TOLERANCE*TOLERANCE);
241 }
242 } // end second scope
243 } // end testExodusWriteElementDataFromDiscontinuousNodalData
244
245 #endif // !LIBMESH_USE_COMPLEX_NUMBERS
246 #endif // LIBMESH_HAVE_EXODUS_API
247
testDynaReadElem()248 void testDynaReadElem ()
249 {
250 Mesh mesh(*TestCommWorld);
251
252 DynaIO dyna(mesh);
253
254 // Make DynaIO::add_spline_constraints work on DistributedMesh
255 mesh.allow_renumbering(false);
256 mesh.allow_remote_element_removal(false);
257
258 if (mesh.processor_id() == 0)
259 dyna.read("1_quad.dyn");
260 MeshCommunication().broadcast (mesh);
261
262 mesh.prepare_for_use();
263
264 // We have 1 QUAD9 finite element, attached via a trivial map to 9
265 // spline Node+NodeElem objects
266 CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), dof_id_type(10));
267 CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), dof_id_type(18));
268
269 CPPUNIT_ASSERT_EQUAL(mesh.default_mapping_type(),
270 RATIONAL_BERNSTEIN_MAP);
271
272 unsigned char weight_index = mesh.default_mapping_data();
273
274 for (auto & elem : mesh.element_ptr_range())
275 {
276 if (elem->type() == NODEELEM)
277 continue;
278
279 CPPUNIT_ASSERT_EQUAL(elem->type(), QUAD9);
280 for (unsigned int n=0; n != 9; ++n)
281 CPPUNIT_ASSERT_EQUAL
282 (elem->node_ref(n).get_extra_datum<Real>(weight_index),
283 Real(0.75));
284
285 CPPUNIT_ASSERT_EQUAL(elem->point(0)(0), Real(0.5));
286 CPPUNIT_ASSERT_EQUAL(elem->point(0)(1), Real(0.5));
287 CPPUNIT_ASSERT_EQUAL(elem->point(1)(0), Real(1.5));
288 CPPUNIT_ASSERT_EQUAL(elem->point(1)(1), Real(0.5));
289 CPPUNIT_ASSERT_EQUAL(elem->point(2)(0), Real(1.5));
290 CPPUNIT_ASSERT_EQUAL(elem->point(2)(1), Real(1.5));
291 CPPUNIT_ASSERT_EQUAL(elem->point(3)(0), Real(0.5));
292 CPPUNIT_ASSERT_EQUAL(elem->point(3)(1), Real(1.5));
293 CPPUNIT_ASSERT(elem->has_affine_map());
294 #if LIBMESH_DIM > 2
295 for (unsigned int v=0; v != 4; ++v)
296 CPPUNIT_ASSERT_EQUAL(elem->point(v)(2), Real(0));
297 #endif
298 }
299 }
300
301
testDynaReadPatch()302 void testDynaReadPatch ()
303 {
304 Mesh mesh(*TestCommWorld);
305
306 // Make DynaIO::add_spline_constraints work on DistributedMesh
307 mesh.allow_renumbering(false);
308 mesh.allow_remote_element_removal(false);
309
310 DynaIO dyna(mesh);
311 if (mesh.processor_id() == 0)
312 dyna.read("25_quad.bxt");
313 MeshCommunication().broadcast (mesh);
314
315 mesh.prepare_for_use();
316
317 // We have 5^2 QUAD9 elements, with 11^2 nodes,
318 // tied to 49 Node/NodeElem spline nodes
319 CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), dof_id_type(25+49));
320 CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), dof_id_type(121+49));
321
322 CPPUNIT_ASSERT_EQUAL(mesh.default_mapping_type(),
323 RATIONAL_BERNSTEIN_MAP);
324
325 unsigned char weight_index = mesh.default_mapping_data();
326
327 for (const auto & elem : mesh.active_element_ptr_range())
328 {
329 if (elem->type() == NODEELEM)
330 continue;
331 LIBMESH_ASSERT_FP_EQUAL(libmesh_real(0.04), elem->volume(), TOLERANCE);
332
333 for (unsigned int n=0; n != 9; ++n)
334 CPPUNIT_ASSERT_EQUAL
335 (elem->node_ref(n).get_extra_datum<Real>(weight_index),
336 Real(1.0));
337
338 unsigned int n_neighbors = 0, n_neighbors_expected = 2;
339 for (unsigned int side=0; side != 4; ++side)
340 if (elem->neighbor_ptr(side))
341 n_neighbors++;
342 Point c = elem->centroid();
343
344 if (c(0) > 0.2 && c(0) < 0.8)
345 n_neighbors_expected++;
346 if (c(1) > 0.2 && c(1) < 0.8)
347 n_neighbors_expected++;
348
349 CPPUNIT_ASSERT_EQUAL(n_neighbors, n_neighbors_expected);
350 }
351
352 #ifdef LIBMESH_HAVE_SOLVER
353 #ifdef LIBMESH_ENABLE_CONSTRAINTS
354 // Now test whether we can assign the desired constraint equations
355 EquationSystems es(mesh);
356 System & sys = es.add_system<LinearImplicitSystem>("test");
357 sys.add_variable("u", SECOND); // to match QUAD9
358 es.init();
359 dyna.add_spline_constraints(sys.get_dof_map(), 0, 0);
360
361 // We should have a constraint on every FE dof
362 CPPUNIT_ASSERT_EQUAL(sys.get_dof_map().n_constrained_dofs(), dof_id_type(121));
363 #endif // LIBMESH_ENABLE_CONSTRAINTS
364 #endif // LIBMESH_HAVE_SOLVER
365 }
366 };
367
368 CPPUNIT_TEST_SUITE_REGISTRATION( MeshInputTest );
369