libMesh
write_vec_and_scalar.C
Go to the documentation of this file.
1 #include "libmesh/exodusII_io.h"
2 
3 // Ignore unused parameter warnings coming from cppunit headers
4 #include <libmesh/ignore_warnings.h>
5 #include <cppunit/extensions/HelperMacros.h>
6 #include <cppunit/TestCase.h>
7 #include <libmesh/restore_warnings.h>
8 
9 // Basic include files
10 #include "libmesh/equation_systems.h"
11 #include "libmesh/mesh.h"
12 #include "libmesh/mesh_generation.h"
13 #include "libmesh/node.h"
14 #include "libmesh/string_to_enum.h"
15 #include "libmesh/exodusII_io.h"
16 #include "libmesh/explicit_system.h"
17 #include "libmesh/numeric_vector.h"
18 
19 #include "test_comm.h"
20 
21 // THE CPPUNIT_TEST_SUITE_END macro expands to code that involves
22 // std::auto_ptr, which in turn produces -Wdeprecated-declarations
23 // warnings. These can be ignored in GCC as long as we wrap the
24 // offending code in appropriate pragmas. We can't get away with a
25 // single ignore_warnings.h inclusion at the beginning of this file,
26 // since the libmesh headers pull in a restore_warnings.h at some
27 // point. We also don't bother restoring warnings at the end of this
28 // file since it's not a header.
29 #include <libmesh/ignore_warnings.h>
30 
31 // Bring in everything from the libMesh namespace
32 using namespace libMesh;
33 
34 class WriteVecAndScalar : public CppUnit::TestCase
35 {
39 public:
40  CPPUNIT_TEST_SUITE(WriteVecAndScalar);
41 
42  CPPUNIT_TEST(testWrite);
43 
44  CPPUNIT_TEST_SUITE_END();
45 
46  void testWrite()
47  {
49 
50  // We set our initial conditions based on build_square node ids
51  mesh.allow_renumbering(false);
52 
54  mesh, 1, 1, -1., 1., -1., 1., Utility::string_to_enum<ElemType>("TRI6"));
55 
56  // Create an equation systems object.
57  EquationSystems equation_systems(mesh);
58 
59  ExplicitSystem & system = equation_systems.add_system<ExplicitSystem>("Tester");
60  system.add_variable("u", FIRST, NEDELEC_ONE);
61  system.add_variable("v", FIRST, LAGRANGE);
62 
63  // Initialize the system
64  equation_systems.init();
65 
66  // Mimic the initial conditions of u(i)=2*i in the serial case,
67  // but indexed by node id rather than dof id.
68  std::vector<Real> initial_vector({10, 0, 12, 8, 4, 2, 16, 6, 14});
69 
70  NumericVector<Number> & sys_solution = *(system.solution);
71  for (const auto & node : mesh.node_ptr_range())
72  {
74  if (node->n_comp(0,0))
75  {
76  dof_id = node->dof_number(0,0,0);
77  }
78  else
79  {
80  CPPUNIT_ASSERT(node->n_comp(0,1));
81  dof_id = node->dof_number(0,1,0);
82  }
83  if (dof_id >= sys_solution.first_local_index() &&
84  dof_id < sys_solution.last_local_index())
85  sys_solution.set(dof_id, initial_vector[node->id()]);
86  if (mesh.n_processors() == 1)
87  CPPUNIT_ASSERT_EQUAL(initial_vector[node->id()], dof_id*Real(2));
88  }
89  sys_solution.close();
90 
91 #ifdef LIBMESH_HAVE_EXODUS_API
92 
93  // We write the file in the ExodusII format.
94  ExodusII_IO(mesh).write_equation_systems("out.e", equation_systems);
95 
96  // Make sure that the writing is done before the reading starts.
98 
99  // Now read it back in
100  Mesh read_mesh(*TestCommWorld);
101  ExodusII_IO exio(read_mesh);
102  exio.read("out.e");
103 
104  read_mesh.prepare_for_use();
105  EquationSystems es2(read_mesh);
106  ExplicitSystem & sys2 = es2.add_system<ExplicitSystem>("Test");
107 
108  const std::vector<std::string> & nodal_vars(exio.get_nodal_var_names());
109 
110  for (const auto & var_name : nodal_vars)
111  sys2.add_variable(var_name, SECOND, LAGRANGE);
112 
113  es2.init();
114 
115  for (const auto & var_name : nodal_vars)
116  exio.copy_nodal_solution(sys2, var_name, var_name, 1);
117 
118  // VariableGroup optimization means that DoFs iterate over each of
119  // the 3 variables on each node before proceeding to the next
120  // node.
121  // u_x, u_y, v
122  const std::vector<Real> gold_vector({-1, 3, 10,
123  0, 1, 12,
124  2, 0, 14,
125  0, 1.5, 11,
126  0.5, 1, 13,
127  0.5, 1.5, 12,
128  3, 4, 16,
129  3, 1.5, 15,
130  0.5, 4, 13});
131 
132  // Translation from node id to dof indexing order as gets done in
133  // serial
134  const std::vector<dof_id_type>
135  node_reordering({0, 3, 1, 8, 5, 4, 6, 7, 2});
136 
137  Real tol = 1e-12;
138  NumericVector<Number> & sys2_soln(*sys2.solution);
139  for (const auto & node : read_mesh.node_ptr_range())
140  {
141  if (node->processor_id() != mesh.processor_id())
142  continue;
143 
144  // We write out real, imaginary, and magnitude for complex
145  // numbers
146 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
147  const unsigned int v2 = 3;
148 #else
149  const unsigned int v2 = 1;
150 #endif
151  CPPUNIT_ASSERT_EQUAL(node->n_vars(0), 3*v2);
152 
153  const dof_id_type gold_i_ux = node_reordering[node->id()] * 3;
154  const dof_id_type gold_i_uy = gold_i_ux + 1;
155  const dof_id_type gold_i_v = gold_i_uy + 1;
156 
157  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(sys2_soln(node->dof_number(0,0,0))),
158  gold_vector[gold_i_ux], tol);
159  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(sys2_soln(node->dof_number(0,v2,0))),
160  gold_vector[gold_i_uy], tol);
161  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(sys2_soln(node->dof_number(0,2*v2,0))),
162  gold_vector[gold_i_v], tol);
163 
164  // Let's check imaginary parts and magnitude for good measure
165 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
166  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(sys2_soln(node->dof_number(0,1,0))),
167  0.0, tol);
168  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(sys2_soln(node->dof_number(0,2,0))),
169  std::abs(gold_vector[gold_i_ux]), tol);
170  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(sys2_soln(node->dof_number(0,4,0))),
171  0.0, tol);
172  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(sys2_soln(node->dof_number(0,5,0))),
173  std::abs(gold_vector[gold_i_uy]), tol);
174  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(sys2_soln(node->dof_number(0,7,0))),
175  0.0, tol);
176  CPPUNIT_ASSERT_DOUBLES_EQUAL(libmesh_real(sys2_soln(node->dof_number(0,8,0))),
177  std::abs(gold_vector[gold_i_v]), tol);
178 #endif
179  }
180 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
181  }
182 };
183 
T libmesh_real(T a)
double abs(double a)
This is the EquationSystems class.
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 ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
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
processor_id_type n_processors() const
MeshBase & mesh
void copy_nodal_solution(System &system, std::string var_name, unsigned int timestep=1)
Backward compatibility version of function that takes a single variable name.
Definition: exodusII_io.C:78
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.
The libMesh namespace provides an interface to certain functionality in the library.
dof_id_type dof_id
Definition: xdr_io.C:48
CPPUNIT_TEST_SUITE_REGISTRATION(WriteVecAndScalar)
void barrier() const
Pause execution until all processors reach a certain point.
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.
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
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=libmesh_nullptr)
This method implements writing a mesh with data to a specified file where the data is taken from the ...
Definition: mesh_output.C:31
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual SimpleRange< node_iterator > node_ptr_range() libmesh_override
virtual void read(const std::string &name) libmesh_override
This method implements reading a mesh from a specified file.
Definition: exodusII_io.C:116
virtual void init()
Initialize all the systems.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
const std::vector< std::string > & get_nodal_var_names()
Return list of the nodal variable names.
Definition: exodusII_io.C:943
The ExplicitSystem provides only "right hand side" storage, which should be sufficient for solving mo...
processor_id_type processor_id() const
uint8_t dof_id_type
Definition: id_types.h:64