libMesh
Loading...
Searching...
No Matches
fe_test.h
Go to the documentation of this file.
1#ifndef FE_TEST_H
2#define FE_TEST_H
3
4#include "test_comm.h"
5
6#include <libmesh/cell_c0polyhedron.h>
7#include <libmesh/dof_map.h>
8#include <libmesh/elem.h>
9#include <libmesh/equation_systems.h>
10#include <libmesh/face_c0polygon.h>
11#include <libmesh/fe.h>
12#include <libmesh/fe_base.h>
13#include <libmesh/fe_interface.h>
14#include <libmesh/function_base.h>
15#include <libmesh/mesh.h>
16#include <libmesh/mesh_generation.h>
17#include <libmesh/mesh_modification.h>
18#include <libmesh/numeric_vector.h>
19#include <libmesh/system.h>
20#include <libmesh/quadrature_gauss.h>
21
22#include <vector>
23
24#include "libmesh_cppunit.h"
25
26#define FETEST \
27 CPPUNIT_TEST( testFEInterface ); \
28 CPPUNIT_TEST( testU ); \
29 CPPUNIT_TEST( testPartitionOfUnity ); \
30 CPPUNIT_TEST( testGradU ); \
31 CPPUNIT_TEST( testGradUComp ); \
32 CPPUNIT_TEST( testHessU ); \
33 CPPUNIT_TEST( testHessUComp ); \
34 CPPUNIT_TEST( testDualDoesntScreamAndDie ); \
35 CPPUNIT_TEST( testCustomReinit );
36
37using namespace libMesh;
38
39
40class SkewFunc : public FunctionBase<Real>
41{
42 std::unique_ptr<FunctionBase<Real>> clone () const override
43 { return std::make_unique<SkewFunc>(); }
44
45 Real operator() (const Point &,
46 const Real = 0.) override
47 { libmesh_not_implemented(); } // scalar-only API
48
49 // Skew in x based on y, y based on z
50 void operator() (const Point & p,
51 const Real,
52 DenseVector<Real> & output)
53 {
54 output.resize(3);
55 output(0) = p(0);
56#if LIBMESH_DIM > 0
57 output(0) += .1*p(1);
58 output(1) = p(1);
59#endif
60#if LIBMESH_DIM > 1
61 output(1) += .2*p(2);
62 output(2) = p(2);
63#endif
64 }
65};
66
67
68inline
69Number linear_test (const Point& p,
70 const Parameters&,
71 const std::string&,
72 const std::string&)
73{
74 const Real & x = p(0);
75 const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
76 const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
77
78 return x + 0.25*y + 0.0625*z;
79}
80
81inline
82Gradient linear_test_grad (const Point&,
83 const Parameters&,
84 const std::string&,
85 const std::string&)
86{
87 Gradient grad = 1;
88 if (LIBMESH_DIM > 1)
89 grad(1) = 0.25;
90 if (LIBMESH_DIM > 2)
91 grad(2) = 0.0625;
92
93 return grad;
94}
95
96
97inline
98Number quadratic_test (const Point& p,
99 const Parameters&,
100 const std::string&,
101 const std::string&)
102{
103 const Real & x = p(0);
104 const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
105 const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
106
107 return x*x + 0.5*y*y + 0.25*z*z + 0.125*x*y + 0.0625*x*z + 0.03125*y*z;
108}
109
110inline
111Gradient quadratic_test_grad (const Point & p,
112 const Parameters&,
113 const std::string&,
114 const std::string&)
115{
116 const Real & x = p(0);
117 const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
118 const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
119
120 Gradient grad = 2*x + 0.125*y + 0.0625*z;
121 if (LIBMESH_DIM > 1)
122 grad(1) = y + 0.125*x + 0.03125*z;
123 if (LIBMESH_DIM > 2)
124 grad(2) = 0.5*z + 0.0625*x + 0.03125*y;
125
126 return grad;
127}
128
129
130inline
131Number fe_cubic_test (const Point& p,
132 const Parameters&,
133 const std::string&,
134 const std::string&)
135{
136 const Real & x = p(0);
137 const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
138 const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
139
140 return x*(1-x)*(1-x) + x*x*(1-y) + x*(1-y)*(1-z) + y*(1-y)*z + z*(1-z)*(1-z);
141}
142
143inline
144Gradient fe_cubic_test_grad (const Point & p,
145 const Parameters&,
146 const std::string&,
147 const std::string&)
148{
149 const Real & x = p(0);
150 const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
151 const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
152
153 Gradient grad = 3*x*x-4*x+1 + 2*x*(1-y) + (1-y)*(1-z);
154 if (LIBMESH_DIM > 1)
155 grad(1) = -x*x - x*(1-z) + (1-2*y)*z;
156 if (LIBMESH_DIM > 2)
157 grad(2) = -x*(1-y) + y*(1-y) + 3*z*z-4*z+1;
158
159 return grad;
160}
161
162
163inline
164Number fe_quartic_test (const Point& p,
165 const Parameters&,
166 const std::string&,
167 const std::string&)
168{
169 const Real & x = p(0);
170 const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
171 const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
172
173 return x*x*(1-x)*(1-x) + x*x*z*(1-y) + x*(1-x)*(1-y)*(1-z) + (1-x)*y*(1-y)*z + z*z*(1-z)*(1-z);
174}
175
176inline
177Gradient fe_quartic_test_grad (const Point & p,
178 const Parameters&,
179 const std::string&,
180 const std::string&)
181{
182 const Real & x = p(0);
183 const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
184 const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
185
186 Gradient grad = 4*x*x*x-6*x*x+2*x + 2*x*z*(1-y) + (1-2*x)*(1-y)*(1-z) - y*(1-y)*z;
187 if (LIBMESH_DIM > 1)
188 grad(1) = -x*x*z - x*(1-x)*(1-z) + (1-x)*(1-2*y)*z;
189 if (LIBMESH_DIM > 2)
190 grad(2) = x*x*(1-y) - x*(1-x)*(1-y) + (1-x)*y*(1-y) + 4*z*z*z-6*z*z+2*z;
191
192 return grad;
193}
194
195
196
197// Higher order rational bases need uniform weights to exactly
198// represent linears; we can easily try out other functions on
199// tensor product elements.
200static const Real rational_w = 0.75;
201
202inline
203Number rational_test (const Point& p,
204 const Parameters&,
205 const std::string&,
206 const std::string&)
207{
208 const Real & x = p(0);
209 const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
210 const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
211
212 const Real denom = ((1-x)*(1-x)+x*x+2*rational_w*x*(1-x))*
213 ((1-y)*(1-y)+y*y+2*rational_w*y*(1-y))*
214 ((1-z)*(1-z)+z*z+2*rational_w*z*(1-z));
215
216 return (x + 0.25*y + 0.0625*z)/denom;
217}
218
219inline
220Gradient rational_test_grad (const Point& p,
221 const Parameters&,
222 const std::string&,
223 const std::string&)
224{
225 const Real & x = p(0);
226 const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
227 const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
228
229 const Real xpoly = (1-x)*(1-x)+x*x+2*rational_w*x*(1-x);
230 const Real xderiv = -2*(1-x)+2*x+2*rational_w*(1-2*x);
231 const Real ypoly = (1-y)*(1-y)+y*y+2*rational_w*y*(1-y);
232 const Real yderiv = -2*(1-y)+2*y+2*rational_w*(1-2*y);
233 const Real zpoly = (1-z)*(1-z)+z*z+2*rational_w*z*(1-z);
234 const Real zderiv = -2*(1-z)+2*z+2*rational_w*(1-2*z);
235
236 const Real denom = xpoly * ypoly * zpoly;
237
238 const Real numer = (x + 0.25*y + 0.0625*z);
239
240 Gradient grad_n = 1, grad_d = xderiv * ypoly * zpoly;
241 if (LIBMESH_DIM > 1)
242 {
243 grad_n(1) = 0.25;
244 grad_d(1) = xpoly * yderiv * zpoly;
245 }
246 if (LIBMESH_DIM > 2)
247 {
248 grad_n(2) = 0.0625;
249 grad_d(2) = xpoly * ypoly * zderiv;
250 }
251
252 Gradient grad = (grad_n - numer * grad_d / denom) / denom;
253
254 return grad;
255}
256
257
258#define FE_CAN_TEST_CUBIC \
259 (((family != LAGRANGE && family != L2_LAGRANGE) || \
260 (elem_type != TRI7 && elem_type != TET14 && \
261 elem_type != PRISM20 && elem_type != PRISM21 && \
262 elem_type != PYRAMID18)) && order > 2)
263
264
265// Note: type name is used to name the test
266struct Default { };
267
268template <Order order, FEFamily family, ElemType elem_type, unsigned int build_nx, typename CaseName = Default>
269class FETestBase : public CppUnit::TestCase {
270
271protected:
272 std::string libmesh_suite_name;
273
275 std::string _case_name;
276
277 unsigned int _dim, _nx, _ny, _nz;
278 Elem *_elem;
279 std::vector<dof_id_type> _dof_indices;
280 System * _sys;
281 std::unique_ptr<Mesh> _mesh;
282 std::unique_ptr<EquationSystems> _es;
283 std::unique_ptr<FEBase> _fe;
284 std::unique_ptr<QGauss> _qrule;
285
286 Real _value_tol, _grad_tol, _hess_tol;
287
288
289 static RealGradient true_gradient(Point p)
290 {
291 Parameters dummy;
292
293 Gradient true_grad;
294 RealGradient returnval;
295
296 if (family == RATIONAL_BERNSTEIN && order > 1)
297 true_grad = rational_test_grad(p, dummy, "", "");
298 else if (order > 3)
299 true_grad = fe_quartic_test_grad(p, dummy, "", "");
300 else if (FE_CAN_TEST_CUBIC)
301 true_grad = fe_cubic_test_grad(p, dummy, "", "");
302 else if (order > 1)
303 {
304 const Real & x = p(0);
305 const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
306 const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
307
308 true_grad = Gradient(2*x+0.125*y+0.0625*z,
309 y+0.125*x+0.03125*z,
310 0.5*z+0.0625*x+0.03125*y);
311 }
312 else
313 true_grad = Gradient(1.0, 0.25, 0.0625);
314
315 for (unsigned int d=0; d != LIBMESH_DIM; ++d)
316 {
317 CPPUNIT_ASSERT(true_grad(d) ==
318 Number(libmesh_real(true_grad(d))));
319
320 returnval(d) = libmesh_real(true_grad(d));
321 }
322
323 return returnval;
324 }
325
326
327 static RealTensor true_hessian(Point p)
328 {
329 const Real & x = p(0);
330 const Real & y = LIBMESH_DIM > 1 ? p(1) : 0;
331 const Real & z = LIBMESH_DIM > 2 ? p(2) : 0;
332
333 if (order > 3)
334 return RealTensor
335 { 12*x*x-12*x+2+2*z*(1-y)-2*(1-y)*(1-z), -2*x*z-(1-2*x)*(1-z)-(1-2*y)*z, 2*x*(1-y)-(1-2*x)*(1-y)-y*(1-y),
336 -2*x*z-(1-2*x)*(1-z)-(1-2*y)*z, -2*(1-x)*z, -x*x+x*(1-x)+(1-x)*(1-2*y),
337 2*x*(1-y)-(1-2*x)*(1-y)-y*(1-y), -x*x+x*(1-x)+(1-x)*(1-2*y), 12*z*z-12*z+2 };
338 else if (FE_CAN_TEST_CUBIC)
339 return RealTensor
340 { 6*x-4+2*(1-y), -2*x+z-1, y-1,
341 -2*x+z-1, -2*z, x+1-2*y,
342 y-1, x+1-2*y, 6*z-4 };
343 else if (order > 1)
344 return RealTensor
345 { 2, 0.125, 0.0625,
346 0.125, 1, 0.03125,
347 0.0625, 0.03125, 0.5 };
348
349 return RealTensor
350 { 0, 0, 0,
351 0, 0, 0,
352 0, 0, 0 };
353 }
354
355public:
356 void setUp()
357 {
358 _mesh = std::make_unique<Mesh>(*TestCommWorld);
359 _dim = Elem::type_to_dim_map[elem_type];
360 const unsigned int build_ny = (_dim > 1) * build_nx;
361 const unsigned int build_nz = (_dim > 2) * build_nx;
362
363 unsigned char weight_index = 0;
364
365 if (family == RATIONAL_BERNSTEIN)
366 {
367 // Make sure we can handle non-zero weight indices
368 _mesh->add_node_integer("buffer integer");
369
370 // Default weight to 1.0 so we don't get NaN from contains_point
371 // checks with a default GhostPointNeighbors ghosting
372 const Real default_weight = 1.0;
373 weight_index = cast_int<unsigned char>
374 (_mesh->add_node_datum<Real>("rational_weight", true,
375 &default_weight));
376 libmesh_assert_not_equal_to(weight_index, 0);
377
378 // Set mapping data but not type, since here we're testing
379 // rational basis functions in the FE space but testing then
380 // with Lagrange bases for the mapping space.
381 _mesh->set_default_mapping_data(weight_index);
382 }
383
384 if (elem_type == C0POLYGON)
385 {
386 // Build a pentagon by hand for testing purposes. Make this
387 // one non-skewed so a MONOMIAL basis will be able to exactly
388 // represent the polynomials we're projecting.
389
390 _mesh->add_point(Point(0, 0), 0);
391 _mesh->add_point(Point(1, 0), 1);
392 _mesh->add_point(Point(1+cos(2*libMesh::pi/5), sin(2*libMesh::pi/5)), 2);
393 _mesh->add_point(Point(0.5, sin(2*libMesh::pi/5)+sin(libMesh::pi/5)), 3);
394 _mesh->add_point(Point(-cos(2*libMesh::pi/5), sin(2*libMesh::pi/5)), 4);
395
396 std::unique_ptr<Elem> polygon = std::make_unique<C0Polygon>(5);
397 for (auto i : make_range(5))
398 polygon->set_node(i, _mesh->node_ptr(i));
399 polygon->set_id() = 0;
400
401 _mesh->add_elem(std::move(polygon));
402 _mesh->prepare_for_use();
403 }
404 else if (elem_type == C0POLYHEDRON)
405 {
406 // Build a plain cube by hand for testing purposes. We could
407 // upgrade this, but it should be to something non-skewed so a
408 // MONOMIAL basis will be able to exactly represent the
409 // polynomials we're projecting.
410 //
411 // See elem_test.h for the limitations here.
412
413 _mesh->add_point(Point(0, 0, 0), 0);
414 _mesh->add_point(Point(1, 0, 0), 1);
415 _mesh->add_point(Point(1, 1, 0), 2);
416 _mesh->add_point(Point(0, 1, 0), 3);
417 _mesh->add_point(Point(0, 0, 1), 4);
418 _mesh->add_point(Point(1, 0, 1), 5);
419 _mesh->add_point(Point(1, 1, 1), 6);
420 _mesh->add_point(Point(0, 1, 1), 7);
421
422 // Pick case
423 std::vector<std::vector<unsigned int>> nodes_on_side;
424 if (_case_name == "Default")
425 nodes_on_side = { {0, 1, 2, 3}, // min z
426 {0, 1, 5, 4}, // min y
427 {2, 6, 5, 1}, // max x
428 {2, 3, 7, 6}, // max y
429 {0, 4, 7, 3}, // min x
430 {5, 6, 7, 4} }; // max z
431 else if (_case_name == "MidNode")
432 nodes_on_side = { {0, 1, 2, 3}, // min z
433 {0, 1, 5, 4}, // min y
434 {1, 2, 6, 5}, // max x
435 {2, 3, 7, 6}, // max y
436 {3, 0, 4, 7}, // min x
437 {4, 5, 6, 7} }; // max z
438 else
439 libmesh_error_msg("Unknown case name: " + _case_name);
440
441 // Build all the sides.
442 std::vector<std::shared_ptr<Polygon>> sides(nodes_on_side.size());
443
444 for (auto s : index_range(nodes_on_side))
445 {
446 const auto & nodes_on_s = nodes_on_side[s];
447 sides[s] = std::make_shared<C0Polygon>(nodes_on_s.size());
448 for (auto i : index_range(nodes_on_s))
449 sides[s]->set_node(i, _mesh->node_ptr(nodes_on_s[i]));
450 }
451
452 std::unique_ptr<libMesh::Node> mid_elem_node;
453 std::unique_ptr<Elem> polyhedron = std::make_unique<C0Polyhedron>(sides, mid_elem_node);
454 _mesh->add_elem(std::move(polyhedron));
455 if (mid_elem_node)
456 _mesh->add_node(std::move(mid_elem_node));
457 _mesh->prepare_for_use();
458 }
459 else
460 {
461 MeshTools::Generation::build_cube (*_mesh,
462 build_nx, build_ny, build_nz,
463 0., 1., 0., 1., 0., 1.,
464 elem_type);
465 }
466
467 // For debugging purposes it can be helpful to only consider one
468 // element even when we're using an element type that requires
469 // more than one element to fill out a square or cube.
470#if 0
471 for (dof_id_type i = 0; i != _mesh->max_elem_id(); ++i)
472 {
473 Elem * elem = _mesh->query_elem_ptr(i);
474 if (elem && elem->id())
475 _mesh->delete_elem(elem);
476 }
477 _mesh->prepare_for_use();
478 CPPUNIT_ASSERT_EQUAL(_mesh->n_elem(), dof_id_type(1));
479#endif
480
481 // Permute our elements randomly and rotate and skew our mesh so
482 // we test all sorts of orientations ... except with Hermite
483 // elements, which are only designed to support meshes with a
484 // single orientation shared by all elements. We're also not
485 // rotating and/or skewing the rational elements, since our test
486 // solution was designed for a specific weighted mesh.
487 if (family != HERMITE &&
488 family != RATIONAL_BERNSTEIN)
489 {
490 MeshTools::Modification::permute_elements(*_mesh);
491
492 // Not yet testing manifolds embedded in higher-D space
493 if (_dim > 1)
494 MeshTools::Modification::rotate(*_mesh, 4,
495 8*(_dim>2), 16*(_dim>2));
496
497 SkewFunc skew_func;
498 MeshTools::Modification::redistribute(*_mesh, skew_func);
499 }
500
501 // Set rational weights so we can exactly match our test solution
502 if (family == RATIONAL_BERNSTEIN)
503 {
504 for (auto elem : _mesh->active_element_ptr_range())
505 {
506 // Workaround for nvc++ bug
507 Elem & e = *elem;
508
509 const unsigned int nv = elem->n_vertices();
510 const unsigned int nn = elem->n_nodes();
511 // We want interiors in lower dimensional elements treated
512 // like edges or faces as appropriate.
513 const unsigned int n_edges =
514 (elem->type() == EDGE3) ? 1 : elem->n_edges();
515 const unsigned int n_faces =
516 (elem->type() == QUAD9) ? 1 : elem->n_faces();
517 const unsigned int nve = std::min(nv + n_edges, nn);
518 const unsigned int nvef = std::min(nve + n_faces, nn);
519
520 for (unsigned int i = 0; i != nv; ++i)
521 e.node_ref(i).set_extra_datum<Real>(weight_index, 1.);
522 for (unsigned int i = nv; i != nve; ++i)
523 e.node_ref(i).set_extra_datum<Real>(weight_index, rational_w);
524 const Real w2 = rational_w * rational_w;
525 for (unsigned int i = nve; i != nvef; ++i)
526 e.node_ref(i).set_extra_datum<Real>(weight_index, w2);
527 const Real w3 = rational_w * w2;
528 for (unsigned int i = nvef; i != nn; ++i)
529 e.node_ref(i).set_extra_datum<Real>(weight_index, w3);
530 }
531 }
532
533 _mesh->complete_preparation();
534 _es = std::make_unique<EquationSystems>(*_mesh);
535 _sys = &(_es->add_system<System> ("SimpleSystem"));
536 _sys->add_variable("u", order, family);
537 _es->init();
538
539 if (family == RATIONAL_BERNSTEIN && order > 1)
540 {
541 _sys->project_solution(rational_test, rational_test_grad, _es->parameters);
542 }
543 else if (order > 3)
544 {
545 _sys->project_solution(fe_quartic_test, fe_quartic_test_grad, _es->parameters);
546 }
547 // Lagrange "cubic" on Tet7 only supports a bubble function, not
548 // all of P^3
549 else if (FE_CAN_TEST_CUBIC)
550 {
551 _sys->project_solution(fe_cubic_test, fe_cubic_test_grad, _es->parameters);
552 }
553 else if (order > 1)
554 {
555 _sys->project_solution(quadratic_test, quadratic_test_grad, _es->parameters);
556 }
557 else
558 {
559 _sys->project_solution(linear_test, linear_test_grad, _es->parameters);
560 }
561
562 FEType fe_type = _sys->variable_type(0);
563 _fe = FEBase::build(_dim, fe_type);
564
565 // Create quadrature rule for use in computing dual shape coefficients
566 _qrule = std::make_unique<QGauss>(_dim, fe_type.default_quadrature_order());
567 _fe->attach_quadrature_rule(_qrule.get());
568
569 auto rng = _mesh->active_local_element_ptr_range();
570 this->_elem = rng.begin() == rng.end() ? nullptr : *(rng.begin());
571
572 _sys->get_dof_map().dof_indices(this->_elem, _dof_indices);
573
574 _nx = 10;
575 _ny = (_dim > 1) ? _nx : 1;
576 _nz = (_dim > 2) ? _nx : 1;
577
578 this->_value_tol = TOLERANCE * sqrt(TOLERANCE);
579
580 // We see 6.5*tol*sqrt(tol) errors on cubic Hermites with the fe_cubic
581 // hermite test function
582 // On Tri7 we see 10*tol*sqrt(tol) errors, even!
583 // On Tet14 Monomial gives us at least 13*tol*sqrt(tol) in some
584 // cases
585 this->_grad_tol = 15 * TOLERANCE * sqrt(TOLERANCE);
586
587 this->_hess_tol = sqrt(TOLERANCE); // FIXME: we see some ~1e-5 errors?!?
588
589 // Prerequest everything we'll want to calculate later.
590 _fe->get_phi();
591 _fe->get_dphi();
592 _fe->get_dphidx();
593#if LIBMESH_DIM > 1
594 _fe->get_dphidy();
595#endif
596#if LIBMESH_DIM > 2
597 _fe->get_dphidz();
598#endif
599
600#if LIBMESH_ENABLE_SECOND_DERIVATIVES
601
602 // Szabab elements don't have second derivatives yet
603 if (family == SZABAB)
604 return;
605
606 _fe->get_d2phi();
607 _fe->get_d2phidx2();
608#if LIBMESH_DIM > 1
609 _fe->get_d2phidxdy();
610 _fe->get_d2phidy2();
611#endif
612#if LIBMESH_DIM > 2
613 _fe->get_d2phidxdz();
614 _fe->get_d2phidydz();
615 _fe->get_d2phidz2();
616#endif
617
618#endif
619 }
620
621 void tearDown() {}
622};
623
624
625
626template <Order order, FEFamily family, ElemType elem_type, typename CaseName>
627class FETest : public FETestBase<order, family, elem_type, 1, CaseName> {
628
629public:
630
631 template <typename Functor>
632 void testLoop(Functor f)
633 {
634 // Handle the "more processors than elements" case
635 if (!this->_elem)
636 return;
637
638 // These tests require exceptions to be enabled because a
639 // TypeTensor::solve() call down in Elem::contains_point()
640 // actually throws a non-fatal exception for a certain Point which
641 // is not in the Elem. When exceptions are not enabled, this test
642 // simply aborts.
643#ifdef LIBMESH_ENABLE_EXCEPTIONS
644 for (unsigned int i=0; i != this->_nx; ++i)
645 for (unsigned int j=0; j != this->_ny; ++j)
646 for (unsigned int k=0; k != this->_nz; ++k)
647 {
648 Point p = Real(i)/this->_nx;
649 if (j > 0)
650 p(1) = Real(j)/this->_ny;
651 if (k > 0)
652 p(2) = Real(k)/this->_nz;
653 if (!this->_elem->contains_point(p))
654 continue;
655
656 // If at a singular node, cannot use FEMap::map
657 if (this->_elem->local_singular_node(p) != invalid_uint)
658 continue;
659
660 std::vector<Point> master_points
661 (1, FEMap::inverse_map(this->_dim, this->_elem, p));
662
663 // Reinit at point to test against analytic solution
664 this->_fe->reinit(this->_elem, &master_points);
665
666 f(p);
667 }
668#endif // LIBMESH_ENABLE_EXCEPTIONS
669 }
670
671 void testPartitionOfUnity()
672 {
673 if (!this->_elem)
674 return;
675
676 this->_fe->reinit(this->_elem);
677
678 bool satisfies_partition_of_unity = true;
679 for (const auto qp : make_range(this->_qrule->n_points()))
680 {
681 Real phi_sum = 0;
682 for (std::size_t d = 0; d != this->_dof_indices.size(); ++d)
683 phi_sum += this->_fe->get_phi()[d][qp];
684 if (phi_sum < (1 - TOLERANCE) || phi_sum > (1 + TOLERANCE))
685 {
686 satisfies_partition_of_unity = false;
687 break;
688 }
689 }
690
691 switch (this->_fe->get_family())
692 {
693 case MONOMIAL:
694 {
695 switch (this->_fe->get_order())
696 {
697 case CONSTANT:
698 CPPUNIT_ASSERT(satisfies_partition_of_unity);
699 break;
700
701 default:
702 CPPUNIT_ASSERT(!satisfies_partition_of_unity);
703 break;
704 }
705 break;
706 }
707 case SZABAB:
708 case HIERARCHIC:
709 case L2_HIERARCHIC:
710 {
711 switch (this->_fe->get_order())
712 {
713 case FIRST:
714 CPPUNIT_ASSERT(satisfies_partition_of_unity);
715 break;
716
717 default:
718 CPPUNIT_ASSERT(!satisfies_partition_of_unity);
719 break;
720 }
721 break;
722 }
723
724 case XYZ:
725 case CLOUGH:
726 case HERMITE:
727 {
728 CPPUNIT_ASSERT(!satisfies_partition_of_unity);
729 break;
730 }
731
732 case LAGRANGE:
733 case L2_LAGRANGE:
734 case BERNSTEIN:
735 case RATIONAL_BERNSTEIN:
736 {
737 CPPUNIT_ASSERT(satisfies_partition_of_unity);
738 break;
739 }
740
741 default:
742 CPPUNIT_FAIL("Uncovered FEFamily");
743 }
744 }
745
746
747 void testFEInterface()
748 {
749 LOG_UNIT_TEST;
750
751 // Handle the "more processors than elements" case
752 if (!this->_elem)
753 return;
754
755 this->_fe->reinit(this->_elem);
756
757 const FEType fe_type = this->_sys->variable_type(0);
758
759 unsigned int my_n_dofs = 0;
760
761 switch (this->_elem->dim())
762 {
763 case 0:
764 my_n_dofs = FE<0,family>::n_dofs(this->_elem, order);
765 break;
766 case 1:
767 my_n_dofs = FE<1,family>::n_dofs(this->_elem, order);
768 break;
769 case 2:
770 my_n_dofs = FE<2,family>::n_dofs(this->_elem, order);
771 break;
772 case 3:
773 my_n_dofs = FE<3,family>::n_dofs(this->_elem, order);
774 break;
775 default:
776 libmesh_error();
777 }
778
779 CPPUNIT_ASSERT_EQUAL(
780 FEInterface::n_shape_functions(fe_type, this->_elem),
781 my_n_dofs);
782
783 CPPUNIT_ASSERT_EQUAL(
784 FEInterface::get_continuity(fe_type),
785 this->_fe->get_continuity());
786
787 CPPUNIT_ASSERT_EQUAL(
788 FEInterface::is_hierarchic(fe_type),
789 this->_fe->is_hierarchic());
790 }
791
792 void testU()
793 {
794 LOG_UNIT_TEST;
795
796 auto f = [this](Point p)
797 {
798 Parameters dummy;
799
800 Number u = 0;
801 for (std::size_t d = 0; d != this->_dof_indices.size(); ++d)
802 u += this->_fe->get_phi()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
803
804 Number true_u;
805
806 if (family == RATIONAL_BERNSTEIN && order > 1)
807 true_u = rational_test(p, dummy, "", "");
808 else if (order > 3)
809 true_u = fe_quartic_test(p, dummy, "", "");
810 else if (FE_CAN_TEST_CUBIC)
811 true_u = fe_cubic_test(p, dummy, "", "");
812 else if (order > 1)
813 true_u = p(0)*p(0) + 0.5*p(1)*p(1) + 0.25*p(2)*p(2) +
814 0.125*p(0)*p(1) + 0.0625*p(0)*p(2) + 0.03125*p(1)*p(2);
815 else
816 true_u = p(0) + 0.25*p(1) + 0.0625*p(2);
817
818 LIBMESH_ASSERT_NUMBERS_EQUAL (true_u, u, this->_value_tol);
819 };
820
821 testLoop(f);
822 }
823
824 void testDualDoesntScreamAndDie()
825 {
826 LOG_UNIT_TEST;
827
828 // Handle the "more processors than elements" case
829 if (!this->_elem)
830 return;
831
832 // Request dual calculations
833 this->_fe->get_dual_phi();
834
835 // reinit using the default quadrature rule in order to calculate the dual coefficients
836 this->_fe->reinit(this->_elem);
837 }
838
839
840 void testGradU()
841 {
842 LOG_UNIT_TEST;
843
844 auto f = [this](Point p)
845 {
846 Parameters dummy;
847
848 Gradient grad_u = 0;
849 for (std::size_t d = 0; d != this->_dof_indices.size(); ++d)
850 grad_u += this->_fe->get_dphi()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
851
852 RealGradient true_grad = this->true_gradient(p);
853
854 LIBMESH_ASSERT_NUMBERS_EQUAL
855 (grad_u(0), true_grad(0), this->_grad_tol);
856 if (this->_dim > 1)
857 LIBMESH_ASSERT_NUMBERS_EQUAL
858 (grad_u(1), true_grad(1), this->_grad_tol);
859 if (this->_dim > 2)
860 LIBMESH_ASSERT_NUMBERS_EQUAL
861 (grad_u(2), true_grad(2), this->_grad_tol);
862 };
863
864 testLoop(f);
865 }
866
867 void testGradUComp()
868 {
869 LOG_UNIT_TEST;
870
871 auto f = [this](Point p)
872 {
873 Parameters dummy;
874
875 Number grad_u_x = 0, grad_u_y = 0, grad_u_z = 0;
876 for (std::size_t d = 0; d != this->_dof_indices.size(); ++d)
877 {
878 grad_u_x += this->_fe->get_dphidx()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
879#if LIBMESH_DIM > 1
880 grad_u_y += this->_fe->get_dphidy()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
881#endif
882#if LIBMESH_DIM > 2
883 grad_u_z += this->_fe->get_dphidz()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
884#endif
885 }
886
887 RealGradient true_grad = this->true_gradient(p);
888
889 LIBMESH_ASSERT_NUMBERS_EQUAL(grad_u_x,
890 true_grad(0), this->_grad_tol);
891 if (this->_dim > 1)
892 LIBMESH_ASSERT_NUMBERS_EQUAL
893 (grad_u_y, true_grad(1), this->_grad_tol);
894 if (this->_dim > 2)
895 LIBMESH_ASSERT_NUMBERS_EQUAL
896 (grad_u_z, true_grad(2), this->_grad_tol);
897 };
898
899 testLoop(f);
900 }
901
902
903 void testHessU()
904 {
905 LOG_UNIT_TEST;
906
907 // Szabab elements don't have second derivatives yet
908 if (family == SZABAB)
909 return;
910
911#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
912 auto f = [this](Point p)
913 {
914 Tensor hess_u;
915 for (std::size_t d = 0; d != this->_dof_indices.size(); ++d)
916 hess_u += this->_fe->get_d2phi()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
917
918 // TODO: Yeah we'll test the ugly expressions later.
919 if (family == RATIONAL_BERNSTEIN && order > 1)
920 return;
921
922 RealTensor true_hess = this->true_hessian(p);
923
924 LIBMESH_ASSERT_NUMBERS_EQUAL
925 (true_hess(0,0), hess_u(0,0), this->_hess_tol);
926 if (this->_dim > 1)
927 {
928 LIBMESH_ASSERT_NUMBERS_EQUAL
929 (hess_u(0,1), hess_u(1,0), this->_hess_tol);
930 LIBMESH_ASSERT_NUMBERS_EQUAL
931 (true_hess(0,1), hess_u(0,1), this->_hess_tol);
932 LIBMESH_ASSERT_NUMBERS_EQUAL
933 (true_hess(1,1), hess_u(1,1), this->_hess_tol);
934 }
935 if (this->_dim > 2)
936 {
937 LIBMESH_ASSERT_NUMBERS_EQUAL
938 (hess_u(0,2), hess_u(2,0), this->_hess_tol);
939 LIBMESH_ASSERT_NUMBERS_EQUAL
940 (hess_u(1,2), hess_u(2,1), this->_hess_tol);
941 LIBMESH_ASSERT_NUMBERS_EQUAL
942 (true_hess(0,2), hess_u(0,2), this->_hess_tol);
943 LIBMESH_ASSERT_NUMBERS_EQUAL
944 (true_hess(1,2), hess_u(1,2), this->_hess_tol);
945 LIBMESH_ASSERT_NUMBERS_EQUAL
946 (true_hess(2,2), hess_u(2,2), this->_hess_tol);
947 }
948 };
949
950 testLoop(f);
951#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
952 }
953
954 void testHessUComp()
955 {
956#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
957 LOG_UNIT_TEST;
958
959 // Szabab elements don't have second derivatives yet
960 if (family == SZABAB)
961 return;
962
963 auto f = [this](Point p)
964 {
965 Number hess_u_xx = 0, hess_u_xy = 0, hess_u_yy = 0,
966 hess_u_xz = 0, hess_u_yz = 0, hess_u_zz = 0;
967 for (std::size_t d = 0; d != this->_dof_indices.size(); ++d)
968 {
969 hess_u_xx += this->_fe->get_d2phidx2()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
970#if LIBMESH_DIM > 1
971 hess_u_xy += this->_fe->get_d2phidxdy()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
972 hess_u_yy += this->_fe->get_d2phidy2()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
973#endif
974#if LIBMESH_DIM > 2
975 hess_u_xz += this->_fe->get_d2phidxdz()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
976 hess_u_yz += this->_fe->get_d2phidydz()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
977 hess_u_zz += this->_fe->get_d2phidz2()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
978#endif
979 }
980
981 // TODO: Yeah we'll test the ugly expressions later.
982 if (family == RATIONAL_BERNSTEIN && order > 1)
983 return;
984
985 RealTensor true_hess = this->true_hessian(p);
986
987 LIBMESH_ASSERT_NUMBERS_EQUAL
988 (true_hess(0,0), hess_u_xx, this->_hess_tol);
989 if (this->_dim > 1)
990 {
991 LIBMESH_ASSERT_NUMBERS_EQUAL
992 (true_hess(0,1), hess_u_xy, this->_hess_tol);
993 LIBMESH_ASSERT_NUMBERS_EQUAL
994 (true_hess(1,1), hess_u_yy, this->_hess_tol);
995 }
996 if (this->_dim > 2)
997 {
998 LIBMESH_ASSERT_NUMBERS_EQUAL
999 (true_hess(0,2), hess_u_xz, this->_hess_tol);
1000 LIBMESH_ASSERT_NUMBERS_EQUAL
1001 (true_hess(1,2), hess_u_yz, this->_hess_tol);
1002 LIBMESH_ASSERT_NUMBERS_EQUAL
1003 (true_hess(2,2), hess_u_zz, this->_hess_tol);
1004 }
1005 };
1006
1007 testLoop(f);
1008#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1009 }
1010
1011 void testCustomReinit()
1012 {
1013 LOG_UNIT_TEST;
1014
1015 std::vector<Point> q_points;
1016 std::vector<Real> weights;
1017 q_points.resize(3); weights.resize(3);
1018 q_points[0](0) = 0.0; q_points[0](1) = 0.0; weights[0] = Real(1)/6;
1019 q_points[1](0) = 1.0; q_points[1](1) = 0.0; weights[1] = weights[0];
1020 q_points[2](0) = 0.0; q_points[2](1) = 1.0; weights[2] = weights[0];
1021
1022 FEType fe_type = this->_sys->variable_type(0);
1023 std::unique_ptr<FEBase> fe (FEBase::build(this->_dim, fe_type));
1024 const int extraorder = 3;
1025 std::unique_ptr<QBase> qrule (fe_type.default_quadrature_rule (this->_dim, extraorder));
1026 fe->attach_quadrature_rule (qrule.get());
1027
1028 const std::vector<Point> & q_pos = fe->get_xyz();
1029
1030 for (const auto & elem : this->_mesh->active_local_element_ptr_range()) {
1031 fe->reinit (elem, &q_points, &weights);
1032 CPPUNIT_ASSERT_EQUAL(q_points.size(), std::size_t(3));
1033 CPPUNIT_ASSERT_EQUAL(q_pos.size(), std::size_t(3)); // 6? bug?
1034 }
1035 }
1036
1037};
1038
1039
1040#define INSTANTIATE_FETEST_CASE(order, family, elemtype, case_name) \
1041 class FETest_##order##_##family##_##elemtype##_##case_name : public FETest<order, family, elemtype, case_name> { \
1042 public: \
1043 FETest_##order##_##family##_##elemtype##_##case_name() : \
1044 FETest<order,family,elemtype,case_name>() { \
1045 if (unitlog->summarized_logs_enabled()) \
1046 this->libmesh_suite_name = "FETest"; \
1047 else \
1048 this->libmesh_suite_name = "FETest_" #order "_" #family "_" #elemtype "_" #case_name; \
1049 this->_case_name = #case_name; \
1050 } \
1051 CPPUNIT_TEST_SUITE( FETest_##order##_##family##_##elemtype##_##case_name ); \
1052 FETEST \
1053 CPPUNIT_TEST_SUITE_END(); \
1054 }; \
1055 \
1056 CPPUNIT_TEST_SUITE_REGISTRATION( FETest_##order##_##family##_##elemtype##_##case_name );
1057
1058#define INSTANTIATE_FETEST(order, family, elemtype) \
1059 INSTANTIATE_FETEST_CASE(order, family, elemtype, Default)
1060
1061#endif
FETestBase::_ny
unsigned int _ny
Definition fe_test.h:277
FETestBase::_dim
unsigned int _dim
Definition fe_test.h:277
FETestBase::_hess_tol
Real _hess_tol
Definition fe_test.h:286
FETestBase::_case_name
std::string _case_name
Name of the case, can be used in the test to switch cases.
Definition fe_test.h:275
FETestBase::libmesh_suite_name
std::string libmesh_suite_name
Definition fe_test.h:272
FETestBase::_es
std::unique_ptr< EquationSystems > _es
Definition fe_test.h:282
FETestBase::_fe
std::unique_ptr< FEBase > _fe
Definition fe_test.h:283
FETestBase::_nz
unsigned int _nz
Definition fe_test.h:277
FETestBase::_mesh
std::unique_ptr< Mesh > _mesh
Definition fe_test.h:281
FETestBase::tearDown
void tearDown()
Definition fe_test.h:621
FETestBase::_dof_indices
std::vector< dof_id_type > _dof_indices
Definition fe_test.h:279
FETestBase::_grad_tol
Real _grad_tol
Definition fe_test.h:286
FETestBase::_qrule
std::unique_ptr< QGauss > _qrule
Definition fe_test.h:284
FETestBase::setUp
void setUp()
Definition fe_test.h:356
FETestBase::true_hessian
static RealTensor true_hessian(Point p)
Definition fe_test.h:327
FETestBase::_sys
System * _sys
Definition fe_test.h:280
FETestBase::_elem
Elem * _elem
Definition fe_test.h:278
FETestBase::_value_tol
Real _value_tol
Definition fe_test.h:286
FETestBase::_nx
unsigned int _nx
Definition fe_test.h:277
FETestBase::true_gradient
static RealGradient true_gradient(Point p)
Definition fe_test.h:289
FETest::testPartitionOfUnity
void testPartitionOfUnity()
Definition fe_test.h:671
FETest::testGradU
void testGradU()
Definition fe_test.h:840
FETest::testLoop
void testLoop(Functor f)
Definition fe_test.h:632
FETest::testHessUComp
void testHessUComp()
Definition fe_test.h:954
FETest::testCustomReinit
void testCustomReinit()
Definition fe_test.h:1011
FETest::testFEInterface
void testFEInterface()
Definition fe_test.h:747
FETest::testGradUComp
void testGradUComp()
Definition fe_test.h:867
FETest::testDualDoesntScreamAndDie
void testDualDoesntScreamAndDie()
Definition fe_test.h:824
FETest::testHessU
void testHessU()
Definition fe_test.h:903
FETest::testU
void testU()
Definition fe_test.h:792
SkewFunc::operator()
Real operator()(const Point &, const Real=0.) override
Definition fe_test.h:45
SkewFunc::clone
std::unique_ptr< FunctionBase< Real > > clone() const override
Definition fe_test.h:42
libMesh::DenseVector
Defines a dense vector for use in Finite Element-type computations.
Definition dense_vector.h:59
libMesh::DenseVector::resize
void resize(const unsigned int n)
Resize the vector.
Definition dense_vector.h:396
libMesh::DofMap::dof_indices
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition dof_map.C:2201
libMesh::DofObject::set_extra_datum
void set_extra_datum(const unsigned int index, const T value)
Sets the value on this object of the extra datum associated with index, which should have been obtain...
Definition dof_object.h:1102
libMesh::DofObject::id
dof_id_type id() const
Definition dof_object.h:819
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition elem.h:96
libMesh::Elem::n_vertices
virtual unsigned int n_vertices() const =0
libMesh::Elem::node_ref
const Node & node_ref(const unsigned int i) const
Definition elem.h:2538
libMesh::Elem::contains_point
virtual bool contains_point(const Point &p, Real tol=TOLERANCE) const
Definition elem.C:2784
libMesh::Elem::dim
virtual unsigned short dim() const =0
libMesh::Elem::type_to_dim_map
static const unsigned int type_to_dim_map[INVALID_ELEM]
This array maps the integer representation of the ElemType enum to the geometric dimension of the ele...
Definition elem.h:628
libMesh::Elem::local_singular_node
virtual unsigned int local_singular_node(const Point &, const Real=TOLERANCE *TOLERANCE) const
Definition elem.h:1840
libMesh::FEGenericBase::build
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
libMesh::FEInterface::is_hierarchic
static bool is_hierarchic(const FEType &fe_type)
Returns whether or not the input FEType's higher-order shape functions are always hierarchic.
Definition fe_interface.C:2689
libMesh::FEInterface::get_continuity
static FEContinuity get_continuity(const FEType &fe_type)
Returns the input FEType's FEContinuity based on the underlying FEFamily and potentially the Order,...
Definition fe_interface.C:2726
libMesh::FEInterface::n_shape_functions
static unsigned int n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition fe_interface.C:272
libMesh::FEMap::inverse_map
static Point inverse_map(const unsigned int dim, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true, const bool extra_checks=true)
Definition fe_map.C:1512
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition fe_type.h:197
libMesh::FEType::default_quadrature_rule
std::unique_ptr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition fe_type.C:34
libMesh::FEType::default_quadrature_order
Order default_quadrature_order() const
Definition fe_type.h:415
libMesh::FE::n_dofs
static unsigned int n_dofs(const ElemType t, const Order o)
libMesh::FunctionBase
Base class for functors that can be evaluated at a point and (optionally) time.
Definition function_base.h:64
libMesh::Parameters
This class provides the ability to map between arbitrary, user-defined strings and several data types...
Definition parameters.h:75
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:40
libMesh::System
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition system.h:100
libMesh::System::project_solution
void project_solution(FunctionBase< Number > *f, FunctionBase< Gradient > *g=nullptr, std::optional< ConstElemRange > active_local_range=std::nullopt, std::optional< std::vector< unsigned int > > variable_numbers=std::nullopt) const
Projects arbitrary functions onto the current solution.
Definition system_projection.C:1065
libMesh::System::current_local_solution
std::unique_ptr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition system.h:1667
libMesh::System::variable_type
const FEType & variable_type(const unsigned int i) const
Definition system.C:2721
libMesh::System::add_variable
unsigned int add_variable(std::string_view 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:1344
libMesh::System::get_dof_map
const DofMap & get_dof_map() const
Definition system.h:2417
libMesh::TensorValue
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
Definition tensor_value.h:45
libMesh::VectorValue
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
Definition vector_value.h:48
TestCommWorld
Communicator * TestCommWorld
rational_w
static const Real rational_w
Definition fe_test.h:200
fe_cubic_test
Number fe_cubic_test(const Point &p, const Parameters &, const std::string &, const std::string &)
Definition fe_test.h:131
rational_test
Number rational_test(const Point &p, const Parameters &, const std::string &, const std::string &)
Definition fe_test.h:203
quadratic_test_grad
Gradient quadratic_test_grad(const Point &p, const Parameters &, const std::string &, const std::string &)
Definition fe_test.h:111
quadratic_test
Number quadratic_test(const Point &p, const Parameters &, const std::string &, const std::string &)
Definition fe_test.h:98
fe_quartic_test_grad
Gradient fe_quartic_test_grad(const Point &p, const Parameters &, const std::string &, const std::string &)
Definition fe_test.h:177
linear_test_grad
Gradient linear_test_grad(const Point &, const Parameters &, const std::string &, const std::string &)
Definition fe_test.h:82
fe_quartic_test
Number fe_quartic_test(const Point &p, const Parameters &, const std::string &, const std::string &)
Definition fe_test.h:164
rational_test_grad
Gradient rational_test_grad(const Point &p, const Parameters &, const std::string &, const std::string &)
Definition fe_test.h:220
fe_cubic_test_grad
Gradient fe_cubic_test_grad(const Point &p, const Parameters &, const std::string &, const std::string &)
Definition fe_test.h:144
linear_test
Number linear_test(const Point &p, const Parameters &, const std::string &, const std::string &)
Definition fe_test.h:69
libmesh_cppunit.h
libMesh::MeshTools::Generation::build_cube
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.
Definition mesh_generation.C:317
libMesh::MeshTools::Modification::redistribute
void redistribute(MeshBase &mesh, const FunctionBase< Real > &mapfunc)
Deterministically perturb the nodal locations.
Definition mesh_modification.C:291
libMesh::MeshTools::Modification::rotate
RealTensorValue rotate(MeshBase &mesh, const Real phi, const Real theta=0., const Real psi=0.)
Rotates the mesh in 3D space.
Definition mesh_modification.C:371
libMesh::MeshTools::Modification::permute_elements
void permute_elements(MeshBase &mesh)
Randomly permute the nodal ordering of each element (without twisting the element mapping).
Definition mesh_modification.C:233
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::index_range
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:153
libMesh::Gradient
NumberVectorValue Gradient
Definition exact_solution.h:51
libMesh::libmesh_real
T libmesh_real(T a)
Definition libmesh_common.h:189
libMesh::invalid_uint
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition libmesh.h:303
libMesh::pi
const Real pi
.
Definition libmesh.h:292
libMesh::TOLERANCE
static constexpr Real TOLERANCE
Definition libmesh_common.h:151
libMesh::dof_id_type
uint8_t dof_id_type
Definition id_types.h:67
libMesh::RATIONAL_BERNSTEIN
@ RATIONAL_BERNSTEIN
Definition enum_fe_family.h:69
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:144
libMesh::make_range
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:176
test_comm.h
Generated on Fri Aug 21 2026 13:32:31 for libMesh by doxygen 1.9.8