6 #include <libmesh/dof_map.h> 7 #include <libmesh/elem.h> 8 #include <libmesh/equation_systems.h> 9 #include <libmesh/fe_base.h> 10 #include <libmesh/fe_interface.h> 11 #include <libmesh/function_base.h> 12 #include <libmesh/mesh.h> 13 #include <libmesh/mesh_generation.h> 14 #include <libmesh/mesh_modification.h> 15 #include <libmesh/numeric_vector.h> 16 #include <libmesh/system.h> 17 #include <libmesh/quadrature_gauss.h> 24 CPPUNIT_TEST( testFEInterface ); \ 25 CPPUNIT_TEST( testU ); \ 26 CPPUNIT_TEST( testPartitionOfUnity ); \ 27 CPPUNIT_TEST( testGradU ); \ 28 CPPUNIT_TEST( testGradUComp ); \ 29 CPPUNIT_TEST( testHessU ); \ 30 CPPUNIT_TEST( testHessUComp ); \ 31 CPPUNIT_TEST( testDualDoesntScreamAndDie ); \ 32 CPPUNIT_TEST( testCustomReinit ); 39 std::unique_ptr<FunctionBase<Real>>
clone ()
const override 40 {
return std::make_unique<SkewFunc>(); }
43 const Real = 0.)
override 44 { libmesh_not_implemented(); }
47 void operator() (
const Point & p,
71 const Real & x = p(0);
72 const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
73 const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
75 return x + 0.25*y + 0.0625*z;
100 const Real & x = p(0);
101 const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
102 const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
104 return x*x + 0.5*y*y + 0.25*z*z + 0.125*x*y + 0.0625*x*z + 0.03125*y*z;
113 const Real & x = p(0);
114 const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
115 const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
117 Gradient grad = 2*x + 0.125*y + 0.0625*z;
119 grad(1) = y + 0.125*x + 0.03125*z;
121 grad(2) = 0.5*z + 0.0625*x + 0.03125*y;
133 const Real & x = p(0);
134 const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
135 const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
137 return x*(1-x)*(1-x) + x*x*(1-y) + x*(1-y)*(1-z) + y*(1-y)*z + z*(1-z)*(1-z);
146 const Real & x = p(0);
147 const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
148 const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
150 Gradient grad = 3*x*x-4*x+1 + 2*x*(1-y) + (1-y)*(1-z);
152 grad(1) = -x*x - x*(1-z) + (1-2*y)*z;
154 grad(2) = -x*(1-y) + y*(1-y) + 3*z*z-4*z+1;
166 const Real & x = p(0);
167 const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
168 const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
170 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);
179 const Real & x = p(0);
180 const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
181 const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
183 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;
185 grad(1) = -x*x*z - x*(1-x)*(1-z) + (1-x)*(1-2*y)*z;
187 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;
205 const Real & x = p(0);
206 const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
207 const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
213 return (x + 0.25*y + 0.0625*z)/denom;
222 const Real & x = p(0);
223 const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
224 const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
233 const Real denom = xpoly * ypoly * zpoly;
235 const Real numer = (x + 0.25*y + 0.0625*z);
237 Gradient grad_n = 1, grad_d = xderiv * ypoly * zpoly;
241 grad_d(1) = xpoly * yderiv * zpoly;
246 grad_d(2) = xpoly * ypoly * zderiv;
249 Gradient grad = (grad_n - numer * grad_d / denom) / denom;
255 #define FE_CAN_TEST_CUBIC \ 256 (((family != LAGRANGE && family != L2_LAGRANGE) || \ 257 (elem_type != TRI7 && elem_type != TET14 && \ 258 elem_type != PRISM20 && elem_type != PRISM21 && \ 259 elem_type != PYRAMID18)) && order > 2) 262 template <Order order, FEFamily family, ElemType elem_type,
unsigned int build_nx>
268 unsigned int _dim, _nx, _ny,
_nz;
273 std::unique_ptr<EquationSystems>
_es;
274 std::unique_ptr<FEBase>
_fe;
291 else if (FE_CAN_TEST_CUBIC)
295 const Real & x = p(0);
296 const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
297 const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
299 true_grad =
Gradient(2*x+0.125*y+0.0625*z,
301 0.5*z+0.0625*x+0.03125*y);
304 true_grad =
Gradient(1.0, 0.25, 0.0625);
306 for (
unsigned int d=0; d != LIBMESH_DIM; ++d)
308 CPPUNIT_ASSERT(true_grad(d) ==
320 const Real & x = p(0);
321 const Real & y = LIBMESH_DIM > 1 ? p(1) : 0;
322 const Real & z = LIBMESH_DIM > 2 ? p(2) : 0;
326 { 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),
327 -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),
328 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 };
329 else if (FE_CAN_TEST_CUBIC)
331 { 6*x-4+2*(1-y), -2*x+z-1, y-1,
332 -2*x+z-1, -2*z, x+1-2*y,
333 y-1, x+1-2*y, 6*z-4 };
338 0.0625, 0.03125, 0.5 };
349 _mesh = std::make_unique<Mesh>(*TestCommWorld);
350 const std::unique_ptr<Elem> test_elem =
Elem::build(elem_type);
351 _dim = test_elem->dim();
352 const unsigned int build_ny = (_dim > 1) * build_nx;
353 const unsigned int build_nz = (_dim > 2) * build_nx;
355 unsigned char weight_index = 0;
360 _mesh->add_node_integer(
"buffer integer");
364 const Real default_weight = 1.0;
365 weight_index = cast_int<unsigned char>
366 (_mesh->add_node_datum<
Real>(
"rational_weight",
true,
368 libmesh_assert_not_equal_to(weight_index, 0);
373 _mesh->set_default_mapping_data(weight_index);
377 build_nx, build_ny, build_nz,
378 0., 1., 0., 1., 0., 1.,
385 for (
dof_id_type i = 0; i != _mesh->max_elem_id(); ++i)
387 Elem * elem = _mesh->query_elem_ptr(i);
388 if (elem && elem->
id())
389 _mesh->delete_elem(elem);
391 _mesh->prepare_for_use();
392 CPPUNIT_ASSERT_EQUAL(_mesh->n_elem(),
dof_id_type(1));
409 8*(_dim>2), 16*(_dim>2));
418 for (
auto elem : _mesh->active_element_ptr_range())
420 const unsigned int nv = elem->n_vertices();
421 const unsigned int nn = elem->n_nodes();
424 const unsigned int n_edges =
425 (elem->type() ==
EDGE3) ? 1 : elem->n_edges();
426 const unsigned int n_faces =
427 (elem->type() ==
QUAD9) ? 1 : elem->n_faces();
428 const unsigned int nve = std::min(nv + n_edges, nn);
429 const unsigned int nvef = std::min(nve + n_faces, nn);
431 for (
unsigned int i = 0; i != nv; ++i)
432 elem->node_ref(i).set_extra_datum<
Real>(weight_index, 1.);
433 for (
unsigned int i = nv; i != nve; ++i)
434 elem->node_ref(i).set_extra_datum<
Real>(weight_index,
rational_w);
436 for (
unsigned int i = nve; i != nvef; ++i)
437 elem->node_ref(i).set_extra_datum<
Real>(weight_index, w2);
439 for (
unsigned int i = nvef; i != nn; ++i)
440 elem->node_ref(i).set_extra_datum<
Real>(weight_index, w3);
444 _es = std::make_unique<EquationSystems>(*_mesh);
445 _sys = &(_es->add_system<
System> (
"SimpleSystem"));
459 else if (FE_CAN_TEST_CUBIC)
477 _fe->attach_quadrature_rule(_qrule.get());
479 auto rng = _mesh->active_local_element_ptr_range();
480 this->_elem = rng.begin() == rng.end() ? nullptr : *(rng.begin());
485 _ny = (_dim > 1) ? _nx : 1;
486 _nz = (_dim > 2) ? _nx : 1;
508 #if LIBMESH_ENABLE_SECOND_DERIVATIVES 517 _fe->get_d2phidxdy();
521 _fe->get_d2phidxdz();
522 _fe->get_d2phidydz();
534 template <Order order, FEFamily family, ElemType elem_type>
539 template <
typename Functor>
551 #ifdef LIBMESH_ENABLE_EXCEPTIONS 552 for (
unsigned int i=0; i != this->_nx; ++i)
553 for (
unsigned int j=0; j != this->_ny; ++j)
554 for (
unsigned int k=0; k != this->_nz; ++k)
558 p(1) =
Real(j)/this->_ny;
560 p(2) =
Real(k)/this->_nz;
561 if (!this->_elem->contains_point(p))
565 if (this->_elem->local_singular_node(p) !=
invalid_uint)
568 std::vector<Point> master_points
572 this->_fe->reinit(this->_elem, &master_points);
576 #endif // LIBMESH_ENABLE_EXCEPTIONS 584 this->_fe->reinit(this->_elem);
586 bool satisfies_partition_of_unity =
true;
587 for (
const auto qp :
make_range(this->_qrule->n_points()))
590 for (std::size_t d = 0; d != this->_dof_indices.size(); ++d)
591 phi_sum += this->_fe->get_phi()[d][qp];
594 satisfies_partition_of_unity =
false;
599 switch (this->_fe->get_family())
603 switch (this->_fe->get_order())
606 CPPUNIT_ASSERT(satisfies_partition_of_unity);
610 CPPUNIT_ASSERT(!satisfies_partition_of_unity);
619 switch (this->_fe->get_order())
622 CPPUNIT_ASSERT(satisfies_partition_of_unity);
626 CPPUNIT_ASSERT(!satisfies_partition_of_unity);
636 CPPUNIT_ASSERT(!satisfies_partition_of_unity);
645 CPPUNIT_ASSERT(satisfies_partition_of_unity);
650 CPPUNIT_FAIL(
"Uncovered FEFamily");
663 this->_fe->reinit(this->_elem);
665 const FEType fe_type = this->_sys->variable_type(0);
667 CPPUNIT_ASSERT_EQUAL(
669 this->_fe->n_shape_functions());
671 CPPUNIT_ASSERT_EQUAL(
673 this->_fe->get_continuity());
675 CPPUNIT_ASSERT_EQUAL(
677 this->_fe->is_hierarchic());
684 auto f = [
this](
Point p)
689 for (std::size_t d = 0; d != this->_dof_indices.size(); ++d)
690 u += this->_fe->get_phi()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
698 else if (FE_CAN_TEST_CUBIC)
701 true_u = p(0)*p(0) + 0.5*p(1)*p(1) + 0.25*p(2)*p(2) +
702 0.125*p(0)*p(1) + 0.0625*p(0)*p(2) + 0.03125*p(1)*p(2);
704 true_u = p(0) + 0.25*p(1) + 0.0625*p(2);
706 LIBMESH_ASSERT_FP_EQUAL
722 this->_fe->get_dual_phi();
725 this->_fe->reinit(this->_elem);
733 auto f = [
this](
Point p)
738 for (std::size_t d = 0; d != this->_dof_indices.size(); ++d)
739 grad_u += this->_fe->get_dphi()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
744 true_grad(0), this->_grad_tol);
747 true_grad(1), this->_grad_tol);
750 true_grad(2), this->_grad_tol);
760 auto f = [
this](
Point p)
764 Number grad_u_x = 0, grad_u_y = 0, grad_u_z = 0;
765 for (std::size_t d = 0; d != this->_dof_indices.size(); ++d)
767 grad_u_x += this->_fe->get_dphidx()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
769 grad_u_y += this->_fe->get_dphidy()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
772 grad_u_z += this->_fe->get_dphidz()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
779 true_grad(0), this->_grad_tol);
782 true_grad(1), this->_grad_tol);
785 true_grad(2), this->_grad_tol);
800 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 801 auto f = [
this](
Point p)
804 for (std::size_t d = 0; d != this->_dof_indices.size(); ++d)
805 hess_u += this->_fe->get_d2phi()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
813 LIBMESH_ASSERT_FP_EQUAL(true_hess(0,0),
libmesh_real(hess_u(0,0)),
819 LIBMESH_ASSERT_FP_EQUAL(true_hess(0,1),
libmesh_real(hess_u(0,1)),
821 LIBMESH_ASSERT_FP_EQUAL(true_hess(1,1),
libmesh_real(hess_u(1,1)),
830 LIBMESH_ASSERT_FP_EQUAL(true_hess(0,2),
libmesh_real(hess_u(0,2)),
832 LIBMESH_ASSERT_FP_EQUAL(true_hess(1,2),
libmesh_real(hess_u(1,2)),
834 LIBMESH_ASSERT_FP_EQUAL(true_hess(2,2),
libmesh_real(hess_u(2,2)),
840 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES 845 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 852 auto f = [
this](
Point p)
854 Number hess_u_xx = 0, hess_u_xy = 0, hess_u_yy = 0,
855 hess_u_xz = 0, hess_u_yz = 0, hess_u_zz = 0;
856 for (std::size_t d = 0; d != this->_dof_indices.size(); ++d)
858 hess_u_xx += this->_fe->get_d2phidx2()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
860 hess_u_xy += this->_fe->get_d2phidxdy()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
861 hess_u_yy += this->_fe->get_d2phidy2()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
864 hess_u_xz += this->_fe->get_d2phidxdz()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
865 hess_u_yz += this->_fe->get_d2phidydz()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
866 hess_u_zz += this->_fe->get_d2phidz2()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
876 LIBMESH_ASSERT_FP_EQUAL(true_hess(0,0),
libmesh_real(hess_u_xx),
880 LIBMESH_ASSERT_FP_EQUAL(true_hess(0,1),
libmesh_real(hess_u_xy),
882 LIBMESH_ASSERT_FP_EQUAL(true_hess(1,1),
libmesh_real(hess_u_yy),
887 LIBMESH_ASSERT_FP_EQUAL(true_hess(0,2),
libmesh_real(hess_u_xz),
889 LIBMESH_ASSERT_FP_EQUAL(true_hess(1,2),
libmesh_real(hess_u_yz),
891 LIBMESH_ASSERT_FP_EQUAL(true_hess(2,2),
libmesh_real(hess_u_zz),
897 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES 904 std::vector<Point> q_points;
905 std::vector<Real> weights;
906 q_points.resize(3); weights.resize(3);
907 q_points[0](0) = 0.0; q_points[0](1) = 0.0; weights[0] =
Real(1)/6;
908 q_points[1](0) = 1.0; q_points[1](1) = 0.0; weights[1] = weights[0];
909 q_points[2](0) = 0.0; q_points[2](1) = 1.0; weights[2] = weights[0];
911 FEType fe_type = this->_sys->variable_type(0);
912 std::unique_ptr<FEBase> fe (
FEBase::build(this->_dim, fe_type));
913 const int extraorder = 3;
915 fe->attach_quadrature_rule (qrule.get());
917 const std::vector<Point> & q_pos = fe->get_xyz();
919 for (
const auto & elem : this->_mesh->active_local_element_ptr_range()) {
920 fe->reinit (elem, &q_points, &weights);
921 CPPUNIT_ASSERT_EQUAL(q_points.size(), std::size_t(3));
922 CPPUNIT_ASSERT_EQUAL(q_pos.size(), std::size_t(3));
929 #define INSTANTIATE_FETEST(order, family, elemtype) \ 930 class FETest_##order##_##family##_##elemtype : public FETest<order, family, elemtype> { \ 932 FETest_##order##_##family##_##elemtype() : \ 933 FETest<order,family,elemtype>() { \ 934 if (unitlog->summarized_logs_enabled()) \ 935 this->libmesh_suite_name = "FETest"; \ 937 this->libmesh_suite_name = "FETest_" #order "_" #family "_" #elemtype; \ 939 CPPUNIT_TEST_SUITE( FETest_##order##_##family##_##elemtype ); \ 941 CPPUNIT_TEST_SUITE_END(); \ 944 CPPUNIT_TEST_SUITE_REGISTRATION( FETest_##order##_##family##_##elemtype );
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
std::unique_ptr< FunctionBase< Real > > clone() const override
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
void testDualDoesntScreamAndDie()
This class provides the ability to map between arbitrary, user-defined strings and several data types...
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
static constexpr Real TOLERANCE
std::unique_ptr< QGauss > _qrule
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)
Gradient quadratic_test_grad(const Point &p, const Parameters &, const std::string &, const std::string &)
void resize(const unsigned int n)
Resize the vector.
This is the base class from which all geometric element types are derived.
Gradient rational_test_grad(const Point &p, const Parameters &, const std::string &, const std::string &)
Order default_quadrature_order() const
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
Gradient fe_cubic_test_grad(const Point &p, const Parameters &, const std::string &, const std::string &)
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
The libMesh namespace provides an interface to certain functionality in the library.
Gradient fe_quartic_test_grad(const Point &p, const Parameters &, const std::string &, const std::string &)
Number fe_cubic_test(const Point &p, const Parameters &, const std::string &, const std::string &)
static const Real rational_w
std::unique_ptr< EquationSystems > _es
void project_solution(FunctionBase< Number > *f, FunctionBase< Gradient > *g=nullptr) const
Projects arbitrary functions onto the current solution.
std::unique_ptr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Number linear_test(const Point &p, const Parameters &, const std::string &, const std::string &)
Number rational_test(const Point &p, const Parameters &, const std::string &, const std::string &)
static unsigned int n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
static RealTensor true_hessian(Point p)
Manages consistently variables, degrees of freedom, and coefficient vectors.
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
NumberVectorValue Gradient
static bool is_hierarchic(const FEType &fe_type)
Returns whether or not the input FEType's higher-order shape functions are always hierarchic...
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
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.
Number quadratic_test(const Point &p, const Parameters &, const std::string &, const std::string &)
std::string libmesh_suite_name
std::unique_ptr< Mesh > _mesh
static FEContinuity get_continuity(const FEType &fe_type)
Returns the input FEType's FEContinuity based on the underlying FEFamily and potentially the Order...
const FEType & variable_type(const unsigned int i) const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::vector< dof_id_type > _dof_indices
static RealGradient true_gradient(Point p)
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...
std::unique_ptr< FEBase > _fe
Defines a dense vector for use in Finite Element-type computations.
Number fe_quartic_test(const Point &p, const Parameters &, const std::string &, const std::string &)
Gradient linear_test_grad(const Point &, const Parameters &, const std::string &, const std::string &)
Base class for functors that can be evaluated at a point and (optionally) time.
const DofMap & get_dof_map() const
A Point defines a location in LIBMESH_DIM dimensional Real space.
void testPartitionOfUnity()
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.