libMesh
parsed_fem_function_test.C
Go to the documentation of this file.
1 // libmesh includes
2 #include "libmesh/elem.h"
3 #include "libmesh/equation_systems.h"
4 #include "libmesh/mesh.h"
5 #include "libmesh/mesh_generation.h"
6 #include "libmesh/numeric_vector.h"
7 #include "libmesh/parsed_fem_function.h"
8 #include "libmesh/system.h"
9 
10 #ifdef LIBMESH_HAVE_FPARSER
11 
12 // test includes
13 #include "test_comm.h"
14 #include "libmesh_cppunit.h"
15 
16 // C++ includes
17 #include <memory>
18 
19 
20 using namespace libMesh;
21 
22 class ParsedFEMFunctionTest : public CppUnit::TestCase
23 {
24 public:
25  void setUp() {
26 #if LIBMESH_DIM > 2
27  mesh = std::make_unique<Mesh>(*TestCommWorld);
29  es = std::make_unique<EquationSystems>(*mesh);
30  sys = &(es->add_system<System> ("SimpleSystem"));
31  sys->add_variable("x2");
32  sys->add_variable("x3");
33  sys->add_variable("c05");
34  sys->add_variable("y4");
35  sys->add_variable("xy");
36  sys->add_variable("yz");
37  sys->add_variable("xyz");
38 
39  es->init();
40 
41  NumericVector<Number> & sol = *sys->solution;
42  Elem *elem = mesh->query_elem_ptr(0);
43 
44  if (elem && elem->processor_id() == TestCommWorld->rank())
45  {
46  // Set x2 = 2*x
47  sol.set(elem->node_ref(1).dof_number(0,0,0), 2);
48  sol.set(elem->node_ref(2).dof_number(0,0,0), 2);
49  sol.set(elem->node_ref(5).dof_number(0,0,0), 2);
50  sol.set(elem->node_ref(6).dof_number(0,0,0), 2);
51 
52  // Set x3 = 3*x
53  sol.set(elem->node_ref(1).dof_number(0,1,0), 3);
54  sol.set(elem->node_ref(2).dof_number(0,1,0), 3);
55  sol.set(elem->node_ref(5).dof_number(0,1,0), 3);
56  sol.set(elem->node_ref(6).dof_number(0,1,0), 3);
57 
58  // Set c05 = 0.5
59  sol.set(elem->node_ref(0).dof_number(0,2,0), 0.5);
60  sol.set(elem->node_ref(1).dof_number(0,2,0), 0.5);
61  sol.set(elem->node_ref(2).dof_number(0,2,0), 0.5);
62  sol.set(elem->node_ref(3).dof_number(0,2,0), 0.5);
63  sol.set(elem->node_ref(4).dof_number(0,2,0), 0.5);
64  sol.set(elem->node_ref(5).dof_number(0,2,0), 0.5);
65  sol.set(elem->node_ref(6).dof_number(0,2,0), 0.5);
66  sol.set(elem->node_ref(7).dof_number(0,2,0), 0.5);
67 
68  // Set y4 = 4*y
69  sol.set(elem->node_ref(2).dof_number(0,3,0), 4);
70  sol.set(elem->node_ref(3).dof_number(0,3,0), 4);
71  sol.set(elem->node_ref(6).dof_number(0,3,0), 4);
72  sol.set(elem->node_ref(7).dof_number(0,3,0), 4);
73 
74  // Set xy = x*y
75  sol.set(elem->node_ref(2).dof_number(0,4,0), 1);
76  sol.set(elem->node_ref(6).dof_number(0,4,0), 1);
77 
78  // Set yz = y*z
79  sol.set(elem->node_ref(6).dof_number(0,5,0), 1);
80  sol.set(elem->node_ref(7).dof_number(0,5,0), 1);
81 
82  // Set xyz = x*y*z
83  sol.set(elem->node_ref(6).dof_number(0,6,0), 1);
84  }
85 
86  sol.close();
87  sys->update();
88 
89  c = std::make_unique<FEMContext>(*sys);
90  s = std::make_unique<FEMContext>(*sys);
91  if (elem && elem->processor_id() == TestCommWorld->rank())
92  {
93  c->get_element_fe(0)->get_phi();
94  c->get_element_fe(0)->get_dphi();
95 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
96  c->get_element_fe(0)->get_d2phi();
97 #endif
98  c->pre_fe_reinit(*sys, elem);
99  c->elem_fe_reinit();
100 
101  s->get_side_fe(0)->get_normals(); // Prerequest
102  s->pre_fe_reinit(*sys, elem);
103  s->side = 3;
104  s->side_fe_reinit();
105  }
106 #endif
107  }
108 
109  void tearDown() {
110  c.reset();
111  s.reset();
112  es.reset();
113  mesh.reset();
114  }
115 
116  LIBMESH_CPPUNIT_TEST_SUITE(ParsedFEMFunctionTest);
117 
118 #if LIBMESH_DIM > 2
119  CPPUNIT_TEST(testValues);
120  CPPUNIT_TEST(testGradients);
121 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
122  CPPUNIT_TEST(testHessians);
123 #endif
124  CPPUNIT_TEST(testInlineGetter);
125  CPPUNIT_TEST(testInlineSetter);
126  CPPUNIT_TEST(testNormals);
127 #endif
128 
129  CPPUNIT_TEST_SUITE_END();
130 
131 
132 private:
133  std::unique_ptr<UnstructuredMesh> mesh;
134  std::unique_ptr<EquationSystems> es;
136  std::unique_ptr<FEMContext> c, s;
137 
138  void testValues()
139  {
140  LOG_UNIT_TEST;
141 
142  if (c->has_elem() &&
143  c->get_elem().processor_id() == TestCommWorld->rank())
144  {
145  // Test that we can copy these into vectors
146  std::vector<ParsedFEMFunction<Number>> pfvec;
147 
148  {
149  ParsedFEMFunction<Number> x2(*sys, "x2");
150  ParsedFEMFunction<Number> xy8(*sys, "x2*y4");
151 
152  // Test that move constructor works
153  ParsedFEMFunction<Number> xy8_stolen(std::move(xy8));
154 
155  pfvec.push_back(xy8_stolen);
156 
157  LIBMESH_ASSERT_FP_EQUAL
158  (2.0, libmesh_real(xy8_stolen(*c,Point(0.5,0.5,0.5))), TOLERANCE*TOLERANCE);
159  }
160 
161  LIBMESH_ASSERT_FP_EQUAL
162  (2.0, libmesh_real(pfvec[0](*c,Point(0.5,0.5,0.5))), TOLERANCE*TOLERANCE);
163  }
164  }
165 
166 
168  {
169  LOG_UNIT_TEST;
170 
171  if (c->has_elem() &&
172  c->get_elem().processor_id() == TestCommWorld->rank())
173  {
174  ParsedFEMFunction<Number> c2(*sys, "grad_x_x2");
175 
176  // Test that copy/move assignment fails to compile. Note:
177  // ParsedFEMFunction is neither move-assignable nor
178  // copy-assignable because it contains a const reference.
179  // ParsedFEMFunction<Number> c2_assigned(*sys, "grad_y_xyz");
180  // c2_assigned = c2;
181 
182  LIBMESH_ASSERT_FP_EQUAL
183  (2.0, libmesh_real(c2(*c,Point(0.35,0.45,0.55))), TOLERANCE*TOLERANCE);
184 
185  ParsedFEMFunction<Number> xz(*sys, "grad_y_xyz");
186 
187  LIBMESH_ASSERT_FP_EQUAL
188  (0.1875, libmesh_real(xz(*c,Point(0.25,0.35,0.75))), TOLERANCE*TOLERANCE);
189 
190  ParsedFEMFunction<Number> xyz(*sys, "grad_y_xyz*grad_x_xy");
191 
192  LIBMESH_ASSERT_FP_EQUAL
193  (0.09375, libmesh_real(xyz(*c,Point(0.25,0.5,0.75))), TOLERANCE*TOLERANCE);
194  }
195  }
196 
197 
199  {
200  LOG_UNIT_TEST;
201 
202  if (c->has_elem() &&
203  c->get_elem().processor_id() == TestCommWorld->rank())
204  {
205  ParsedFEMFunction<Number> c1(*sys, "hess_xy_xy");
206 
207  LIBMESH_ASSERT_FP_EQUAL
208  (1.0, libmesh_real(c1(*c,Point(0.35,0.45,0.55))), TOLERANCE*TOLERANCE);
209 
210  ParsedFEMFunction<Number> x(*sys, "hess_yz_xyz");
211 
212  LIBMESH_ASSERT_FP_EQUAL
213  (0.25, libmesh_real(x(*c,Point(0.25,0.35,0.55))), TOLERANCE*TOLERANCE);
214 
215  ParsedFEMFunction<Number> xz(*sys, "hess_yz_xyz*grad_y_yz");
216 
217  LIBMESH_ASSERT_FP_EQUAL
218  (0.1875, libmesh_real(xz(*c,Point(0.25,0.4,0.75))), TOLERANCE*TOLERANCE);
219  }
220  }
221 
223  {
224  LOG_UNIT_TEST;
225 
226  if (c->has_elem() &&
227  c->get_elem().processor_id() == TestCommWorld->rank())
228  {
229  ParsedFEMFunction<Number> ax2(*sys, "a:=4.5;a*x2");
230 
231  LIBMESH_ASSERT_FP_EQUAL
232  (2.25, libmesh_real(ax2(*c,Point(0.25,0.25,0.25))), TOLERANCE*TOLERANCE);
233 
234  LIBMESH_ASSERT_FP_EQUAL
236 
238  (*sys, "a := 4 ; b := a/2+1; c:=b-a+3.5; c*x2*y4");
239 
240  LIBMESH_ASSERT_FP_EQUAL
241  (5.0, libmesh_real(cxy8(*c,Point(0.5,0.5,0.5))), TOLERANCE*TOLERANCE);
242 
243  LIBMESH_ASSERT_FP_EQUAL
245 
246  LIBMESH_ASSERT_FP_EQUAL
248  }
249  }
250 
252  {
253  LOG_UNIT_TEST;
254 
255  if (c->has_elem() &&
256  c->get_elem().processor_id() == TestCommWorld->rank())
257  {
258  ParsedFEMFunction<Number> ax2(*sys, "a:=4.5;a*x2");
259  ax2.set_inline_value("a", 2.5);
260 
261  LIBMESH_ASSERT_FP_EQUAL
262  (1.25, libmesh_real(ax2(*c,Point(0.25,0.25,0.25))), TOLERANCE*TOLERANCE);
263 
264  LIBMESH_ASSERT_FP_EQUAL
266 
268  (*sys, "a := 4 ; b := a/2+1; c:=b-a+3.5; c*x2*y4");
269 
270  cxy8.set_inline_value("a", 2);
271 
272  LIBMESH_ASSERT_FP_EQUAL
273  (7.0, libmesh_real(cxy8(*c,Point(0.5,0.5,0.5))), TOLERANCE*TOLERANCE);
274 
275  LIBMESH_ASSERT_FP_EQUAL
277 
278  LIBMESH_ASSERT_FP_EQUAL
280 
281  }
282  }
283 
284  void testNormals()
285  {
286  LOG_UNIT_TEST;
287 
288  if (s->has_elem() &&
289  s->get_elem().processor_id() == TestCommWorld->rank())
290  {
291  ParsedFEMFunction<Number> nx(*sys, "n_x");
292 
293  ParsedFEMFunction<Number> ny(*sys, "n_y");
294 
295  ParsedFEMFunction<Number> nz(*sys, "n_z");
296 
297  const std::vector<Point> & xyz = s->get_side_fe(0)->get_xyz();
298 
299  // On side 3 of a hex the normal direction is +y
300  for (std::size_t qp=0; qp != xyz.size(); ++qp)
301  {
302  LIBMESH_ASSERT_FP_EQUAL
303  (0.0, libmesh_real(nx(*s,xyz[qp])), TOLERANCE*TOLERANCE);
304  LIBMESH_ASSERT_FP_EQUAL
305  (1.0, libmesh_real(ny(*s,xyz[qp])), TOLERANCE*TOLERANCE);
306  LIBMESH_ASSERT_FP_EQUAL
307  (0.0, libmesh_real(nz(*s,xyz[qp])), TOLERANCE*TOLERANCE);
308  }
309  }
310  }
311 
312 };
313 
315 
316 #endif // #ifdef LIBMESH_HAVE_FPARSER
T libmesh_real(T a)
CPPUNIT_TEST_SUITE_REGISTRATION(ParsedFEMFunctionTest)
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:1025
libMesh::Parallel::Communicator * TestCommWorld
Definition: driver.C:159
static constexpr Real TOLERANCE
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
processor_id_type rank() const
The libMesh namespace provides an interface to certain functionality in the library.
Output get_inline_value(std::string_view inline_var_name) const
const Node & node_ref(const unsigned int i) const
Definition: elem.h:2353
ParsedFEMFunction provides support for FParser-based parsed functions in FEMSystem.
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:96
std::unique_ptr< EquationSystems > es
void set_inline_value(std::string_view inline_var_name, Output newval)
Changes the value of an inline variable.
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
std::unique_ptr< UnstructuredMesh > mesh
virtual const Elem * query_elem_ptr(const dof_id_type i) const =0
std::unique_ptr< FEMContext > s
virtual void set(const numeric_index_type i, const T value)=0
Sets v(i) = value.
processor_id_type processor_id() const
Definition: dof_object.h:898
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
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.