FETest< order, family, elem_type, CaseName > Class Template Reference

#include <fe_test.h>

Inheritance diagram for FETest< order, family, elem_type, CaseName >:

Public Member Functions
template<typename Functor >
void	testLoop (Functor f)
template<typename Functor >
void	testPastVerticesLoop (Functor f)
void	testRefspaceNodes ()
void	testNodalSoln ()
void	testPartitionOfUnity ()
void	testFEInterface ()
auto	checkU (const Real tol)
void	testU ()
void	testUPastVertices ()
void	testDualDoesntScreamAndDie ()
auto	checkGradU (const Real tol)
void	testGradU ()
void	testGradUPastVertices ()
void	testGradUComp ()
void	testHessU ()
void	testHessUComp ()
void	testCustomReinit ()
void	setUp ()
void	tearDown ()

Static Protected Member Functions
static RealGradient	true_gradient (Point p)
static RealTensor	true_hessian (Point p)

Protected Attributes
std::string	libmesh_suite_name
std::string	_case_name
	Name of the case, can be used in the test to switch cases.
unsigned int	_dim
unsigned int	_nx
unsigned int	_ny
unsigned int	_nz
Elem *	_elem
std::vector< dof_id_type >	_dof_indices
System *	_sys
std::unique_ptr< Mesh >	_mesh
std::unique_ptr< EquationSystems >	_es
std::unique_ptr< FEBase >	_fe
std::unique_ptr< QGauss >	_qrule
Real	_value_tol
Real	_grad_tol
Real	_hess_tol
Real	_past_vertex_offset
	The offset, as a fraction of an edge, of the points past the vertices.
Real	_past_vertex_tol_factor
	The factor by which the tolerances grow at the points past the vertices.

Detailed Description

template<Order order, FEFamily family, ElemType elem_type, typename CaseName>
class FETest< order, family, elem_type, CaseName >

Definition at line 655 of file fe_test.h.

Member Function Documentation

checkGradU()

template<Order order, FEFamily family, ElemType elem_type, typename CaseName >

auto FETest< order, family, elem_type, CaseName >::checkGradU ( const Real  tol )

inline

Definition at line 990 of file fe_test.h.

  {
    return [this, tol](Point p)
      {
        Parameters dummy;
 
        Gradient grad_u = 0;
        for (std::size_t d = 0; d != this->_dof_indices.size(); ++d)
          grad_u += this->_fe->get_dphi()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
 
        RealGradient true_grad = this->true_gradient(p);
 
        LIBMESH_ASSERT_NUMBERS_EQUAL
          (grad_u(0), true_grad(0), tol);
        if (this->_dim > 1)
          LIBMESH_ASSERT_NUMBERS_EQUAL
           (grad_u(1), true_grad(1), tol);
        if (this->_dim > 2)
          LIBMESH_ASSERT_NUMBERS_EQUAL
           (grad_u(2), true_grad(2), tol);
      };
  }

References FETestBase< order, family, elem_type, 1, CaseName >::_dim, FETestBase< order, family, elem_type, 1, CaseName >::_dof_indices, FETestBase< order, family, elem_type, 1, CaseName >::_fe, FETestBase< order, family, elem_type, 1, CaseName >::_sys, libMesh::System::current_local_solution, and FETestBase< order, family, elem_type, 1, CaseName >::true_gradient().

Referenced by FETest< order, family, elem_type, CaseName >::testGradU(), and FETest< order, family, elem_type, CaseName >::testGradUPastVertices().

checkU()

template<Order order, FEFamily family, ElemType elem_type, typename CaseName >

auto FETest< order, family, elem_type, CaseName >::checkU ( const Real  tol )

inline

Definition at line 932 of file fe_test.h.

  {
    return [this, tol](Point p)
      {
        Parameters dummy;
 
        Number u = 0;
        for (std::size_t d = 0; d != this->_dof_indices.size(); ++d)
          u += this->_fe->get_phi()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
 
        Number true_u;
 
        if (family == RATIONAL_BERNSTEIN && order > 1)
          true_u = rational_test(p, dummy, "", "");
        else if (order > 3)
          true_u = fe_quartic_test(p, dummy, "", "");
        else if (FE_CAN_TEST_CUBIC)
          true_u = fe_cubic_test(p, dummy, "", "");
        else if (order > 1)
          true_u = p(0)*p(0) + 0.5*p(1)*p(1) + 0.25*p(2)*p(2) +
            0.125*p(0)*p(1) + 0.0625*p(0)*p(2) + 0.03125*p(1)*p(2);
        else
          true_u = p(0) + 0.25*p(1) + 0.0625*p(2);
 
        LIBMESH_ASSERT_NUMBERS_EQUAL (true_u, u, tol);
      };
  }

References FETestBase< order, family, elem_type, 1, CaseName >::_dof_indices, FETestBase< order, family, elem_type, 1, CaseName >::_fe, FETestBase< order, family, elem_type, 1, CaseName >::_sys, libMesh::System::current_local_solution, fe_cubic_test(), fe_quartic_test(), libMesh::RATIONAL_BERNSTEIN, and rational_test().

Referenced by FETest< order, family, elem_type, CaseName >::testU(), and FETest< order, family, elem_type, CaseName >::testUPastVertices().

setUp()

void FETestBase< order, family, elem_type, build_nx, CaseName >::setUp ( )

inlineinherited

Definition at line 367 of file fe_test.h.

  {
    _mesh = std::make_unique<Mesh>(*TestCommWorld);
    _dim = Elem::type_to_dim_map[elem_type];
    const unsigned int build_ny = (_dim > 1) * build_nx;
    const unsigned int build_nz = (_dim > 2) * build_nx;
 
    unsigned char weight_index = 0;
 
    if (family == RATIONAL_BERNSTEIN)
      {
        // Make sure we can handle non-zero weight indices
        _mesh->add_node_integer("buffer integer");
 
        // Default weight to 1.0 so we don't get NaN from contains_point
        // checks with a default GhostPointNeighbors ghosting
        const Real default_weight = 1.0;
        weight_index = cast_int<unsigned char>
          (_mesh->add_node_datum<Real>("rational_weight", true,
                                       &default_weight));
        libmesh_assert_not_equal_to(weight_index, 0);
 
        // Set mapping data but not type, since here we're testing
        // rational basis functions in the FE space but testing then
        // with Lagrange bases for the mapping space.
        _mesh->set_default_mapping_data(weight_index);
      }
 
    if (elem_type == C0POLYGON)
      {
        // Build a pentagon by hand for testing purposes.  Make this
        // one non-skewed so a MONOMIAL basis will be able to exactly
        // represent the polynomials we're projecting.
 
        _mesh->add_point(Point(0, 0), 0);
        _mesh->add_point(Point(1, 0), 1);
        _mesh->add_point(Point(1+cos(2*libMesh::pi/5), sin(2*libMesh::pi/5)), 2);
        _mesh->add_point(Point(0.5, sin(2*libMesh::pi/5)+sin(libMesh::pi/5)), 3);
        _mesh->add_point(Point(-cos(2*libMesh::pi/5), sin(2*libMesh::pi/5)), 4);
 
        std::unique_ptr<Elem> polygon = std::make_unique<C0Polygon>(5);
        for (auto i : make_range(5))
          polygon->set_node(i, _mesh->node_ptr(i));
        polygon->set_id() = 0;
 
        _mesh->add_elem(std::move(polygon));
        _mesh->prepare_for_use();
      }
    else if (elem_type == C0POLYHEDRON)
      {
        // Build a plain cube by hand for testing purposes.  We could
        // upgrade this, but it should be to something non-skewed so a
        // MONOMIAL basis will be able to exactly represent the
        // polynomials we're projecting.
        //
        // See elem_test.h for the limitations here.
 
        _mesh->add_point(Point(0, 0, 0), 0);
        _mesh->add_point(Point(1, 0, 0), 1);
        _mesh->add_point(Point(1, 1, 0), 2);
        _mesh->add_point(Point(0, 1, 0), 3);
        _mesh->add_point(Point(0, 0, 1), 4);
        _mesh->add_point(Point(1, 0, 1), 5);
        _mesh->add_point(Point(1, 1, 1), 6);
        _mesh->add_point(Point(0, 1, 1), 7);
 
        // Pick case
        std::vector<std::vector<unsigned int>> nodes_on_side;
        if (_case_name == "Default")
          nodes_on_side = { {0, 1, 2, 3},   // min z
                            {0, 1, 5, 4},   // min y
                            {2, 6, 5, 1},   // max x
                            {2, 3, 7, 6},   // max y
                            {0, 4, 7, 3},   // min x
                            {5, 6, 7, 4} }; // max z
        else if (_case_name == "MidNode")
          nodes_on_side = { {0, 1, 2, 3},   // min z
                            {0, 1, 5, 4},   // min y
                            {1, 2, 6, 5},   // max x
                            {2, 3, 7, 6},   // max y
                            {3, 0, 4, 7},   // min x
                            {4, 5, 6, 7} }; // max z
        else
          libmesh_error_msg("Unknown case name: " + _case_name);
 
        // Build all the sides.
        std::vector<std::shared_ptr<Polygon>> sides(nodes_on_side.size());
 
        for (auto s : index_range(nodes_on_side))
          {
            const auto & nodes_on_s = nodes_on_side[s];
            sides[s] = std::make_shared<C0Polygon>(nodes_on_s.size());
            for (auto i : index_range(nodes_on_s))
              sides[s]->set_node(i, _mesh->node_ptr(nodes_on_s[i]));
          }
 
        std::unique_ptr<libMesh::Node> mid_elem_node;
        std::unique_ptr<Elem> polyhedron = std::make_unique<C0Polyhedron>(sides, mid_elem_node);
        _mesh->add_elem(std::move(polyhedron));
        if (mid_elem_node)
          _mesh->add_node(std::move(mid_elem_node));
        _mesh->prepare_for_use();
      }
    else
      {
        MeshTools::Generation::build_cube (*_mesh,
                                           build_nx, build_ny, build_nz,
                                           0., 1., 0., 1., 0., 1.,
                                           elem_type);
      }
 
    // For debugging purposes it can be helpful to only consider one
    // element even when we're using an element type that requires
    // more than one element to fill out a square or cube.
#if 0
    for (dof_id_type i = 0; i != _mesh->max_elem_id(); ++i)
      {
        Elem * elem = _mesh->query_elem_ptr(i);
        if (elem && elem->id())
          _mesh->delete_elem(elem);
      }
    _mesh->prepare_for_use();
    CPPUNIT_ASSERT_EQUAL(_mesh->n_elem(), dof_id_type(1));
#endif
 
    // Permute our elements randomly and rotate and skew our mesh so
    // we test all sorts of orientations ... except with Hermite
    // elements, which are only designed to support meshes with a
    // single orientation shared by all elements.  We're also not
    // rotating and/or skewing the rational elements, since our test
    // solution was designed for a specific weighted mesh.
    if (family != HERMITE &&
        family != RATIONAL_BERNSTEIN)
      {
        MeshTools::Modification::permute_elements(*_mesh);
 
        // Not yet testing manifolds embedded in higher-D space
        if (_dim > 1)
          MeshTools::Modification::rotate(*_mesh, 4,
                                          8*(_dim>2), 16*(_dim>2));
 
        SkewFunc skew_func;
        MeshTools::Modification::redistribute(*_mesh, skew_func);
      }
 
    // Set rational weights so we can exactly match our test solution
    if (family == RATIONAL_BERNSTEIN)
      {
        for (auto elem : _mesh->active_element_ptr_range())
          {
            // Workaround for nvc++ bug
            Elem & e = *elem;
 
            const unsigned int nv = elem->n_vertices();
            const unsigned int nn = elem->n_nodes();
            // We want interiors in lower dimensional elements treated
            // like edges or faces as appropriate.
            const unsigned int n_edges =
              (elem->type() == EDGE3) ? 1 : elem->n_edges();
            const unsigned int n_faces =
              (elem->type() == QUAD9) ? 1 : elem->n_faces();
            const unsigned int nve = std::min(nv + n_edges, nn);
            const unsigned int nvef = std::min(nve + n_faces, nn);
 
            for (unsigned int i = 0; i != nv; ++i)
              e.node_ref(i).set_extra_datum<Real>(weight_index, 1.);
            for (unsigned int i = nv; i != nve; ++i)
              e.node_ref(i).set_extra_datum<Real>(weight_index, rational_w);
            const Real w2 = rational_w * rational_w;
            for (unsigned int i = nve; i != nvef; ++i)
              e.node_ref(i).set_extra_datum<Real>(weight_index, w2);
            const Real w3 = rational_w * w2;
            for (unsigned int i = nvef; i != nn; ++i)
              e.node_ref(i).set_extra_datum<Real>(weight_index, w3);
          }
      }
 
    _mesh->complete_preparation();
    _es = std::make_unique<EquationSystems>(*_mesh);
    _sys = &(_es->add_system<System> ("SimpleSystem"));
    _sys->add_variable("u", order, family);
    _es->init();
 
    if (family == RATIONAL_BERNSTEIN && order > 1)
      {
        _sys->project_solution(rational_test, rational_test_grad, _es->parameters);
      }
    else if (order > 3)
      {
        _sys->project_solution(fe_quartic_test, fe_quartic_test_grad, _es->parameters);
      }
    // Lagrange "cubic" on Tet7 only supports a bubble function, not
    // all of P^3
    else if (FE_CAN_TEST_CUBIC)
      {
        _sys->project_solution(fe_cubic_test, fe_cubic_test_grad, _es->parameters);
      }
    else if (order > 1)
      {
        _sys->project_solution(quadratic_test, quadratic_test_grad, _es->parameters);
      }
    else
      {
        _sys->project_solution(linear_test, linear_test_grad, _es->parameters);
      }
 
    FEType fe_type = _sys->variable_type(0);
    _fe = FEBase::build(_dim, fe_type);
 
    // Create quadrature rule for use in computing dual shape coefficients
    _qrule = std::make_unique<QGauss>(_dim, fe_type.default_quadrature_order());
    _fe->attach_quadrature_rule(_qrule.get());
 
    auto rng = _mesh->active_local_element_ptr_range();
    this->_elem = rng.begin() == rng.end() ? nullptr : *(rng.begin());
 
    _sys->get_dof_map().dof_indices(this->_elem, _dof_indices);
 
    _nx = 10;
    _ny = (_dim > 1) ? _nx : 1;
    _nz = (_dim > 2) ? _nx : 1;
 
    this->_value_tol = TOLERANCE * sqrt(TOLERANCE);
 
    // We see 6.5*tol*sqrt(tol) errors on cubic Hermites with the fe_cubic
    // hermite test function
    // On Tri7 we see 10*tol*sqrt(tol) errors, even!
    // On Tet14 Monomial gives us at least 13*tol*sqrt(tol) in some
    // cases
    this->_grad_tol = 15 * TOLERANCE * sqrt(TOLERANCE);
 
    this->_hess_tol = sqrt(TOLERANCE); // FIXME: we see some ~1e-5 errors?!?
 
    // Reference vertices have coordinates of 0 and +/-1, so with a power of two offset the
    // barycentric coordinates of each point, and their sums and differences, are computed
    // exactly. A combination which vanishes on a point is then evaluated as exactly zero.
    this->_past_vertex_offset = 0.25;
 
    // Shape functions grow outside the element, and with them the rounding error in the
    // solution there, so evaluations past the vertices get looser tolerances. Along the line
    // through an edge the error is a polynomial of degree at most p, so an error within a
    // tolerance on the edge is within T_p(1 + 2*offset) times that tolerance at offset past
    // the edge's end, where T_p is the Chebyshev polynomial. Gradient errors have lower degree
    // and grow more slowly.
    const Real x = 1 + 2 * this->_past_vertex_offset;
    Real t_prev = 1, t = x;
    for (unsigned int k = 1; k < static_cast<unsigned int>(order); ++k)
      t_prev = std::exchange(t, 2 * x * t - t_prev);
    this->_past_vertex_tol_factor = (order == CONSTANT) ? 1 : t;
 
    // Prerequest everything we'll want to calculate later.
    _fe->get_phi();
    _fe->get_dphi();
    _fe->get_dphidx();
#if LIBMESH_DIM > 1
    _fe->get_dphidy();
#endif
#if LIBMESH_DIM > 2
    _fe->get_dphidz();
#endif
 
#if LIBMESH_ENABLE_SECOND_DERIVATIVES
 
    // Szabab elements don't have second derivatives yet
    if (family == SZABAB)
      return;
 
    _fe->get_d2phi();
    _fe->get_d2phidx2();
#if LIBMESH_DIM > 1
    _fe->get_d2phidxdy();
    _fe->get_d2phidy2();
#endif
#if LIBMESH_DIM > 2
    _fe->get_d2phidxdz();
    _fe->get_d2phidydz();
    _fe->get_d2phidz2();
#endif
 
#endif
  }

tearDown()

void FETestBase< order, family, elem_type, build_nx, CaseName >::tearDown ( )

inlineinherited

Definition at line 649 of file fe_test.h.

{}

testCustomReinit()

template<Order order, FEFamily family, ElemType elem_type, typename CaseName >

void FETest< order, family, elem_type, CaseName >::testCustomReinit ( )

inline

Definition at line 1169 of file fe_test.h.

  {
    LOG_UNIT_TEST;
 
    std::vector<Point> q_points;
    std::vector<Real> weights;
    q_points.resize(3); weights.resize(3);
    q_points[0](0) = 0.0; q_points[0](1) = 0.0; weights[0] = Real(1)/6;
    q_points[1](0) = 1.0; q_points[1](1) = 0.0; weights[1] = weights[0];
    q_points[2](0) = 0.0; q_points[2](1) = 1.0; weights[2] = weights[0];
 
    FEType fe_type = this->_sys->variable_type(0);
    std::unique_ptr<FEBase> fe (FEBase::build(this->_dim, fe_type));
    const int extraorder = 3;
    std::unique_ptr<QBase> qrule (fe_type.default_quadrature_rule (this->_dim, extraorder));
    fe->attach_quadrature_rule (qrule.get());
 
    const std::vector<Point> & q_pos = fe->get_xyz();
 
    for (const auto & elem : this->_mesh->active_local_element_ptr_range()) {
      fe->reinit (elem, &q_points, &weights);
      CPPUNIT_ASSERT_EQUAL(q_points.size(), std::size_t(3));
      CPPUNIT_ASSERT_EQUAL(q_pos.size(), std::size_t(3));  // 6? bug?
    }
  }

References FETestBase< order, family, elem_type, 1, CaseName >::_dim, FETestBase< order, family, elem_type, 1, CaseName >::_mesh, FETestBase< order, family, elem_type, 1, CaseName >::_sys, libMesh::FEGenericBase< OutputType >::build(), libMesh::FEType::default_quadrature_rule(), libMesh::Real, and libMesh::System::variable_type().

testDualDoesntScreamAndDie()

template<Order order, FEFamily family, ElemType elem_type, typename CaseName >

void FETest< order, family, elem_type, CaseName >::testDualDoesntScreamAndDie ( )

inline

Definition at line 972 of file fe_test.h.

  {
    LOG_UNIT_TEST;
 
    // Handle the "more processors than elements" case
    if (!this->_elem)
      return;
 
    // Request dual calculations
    this->_fe->get_dual_phi();
 
    // reinit using the default quadrature rule in order to calculate the dual coefficients
    this->_fe->reinit(this->_elem);
  }

References FETestBase< order, family, elem_type, 1, CaseName >::_elem, and FETestBase< order, family, elem_type, 1, CaseName >::_fe.

testFEInterface()

template<Order order, FEFamily family, ElemType elem_type, typename CaseName >

void FETest< order, family, elem_type, CaseName >::testFEInterface ( )

inline

Definition at line 885 of file fe_test.h.

  {
    LOG_UNIT_TEST;
 
    // Handle the "more processors than elements" case
    if (!this->_elem)
      return;
 
    this->_fe->reinit(this->_elem);
 
    const FEType fe_type = this->_sys->variable_type(0);
 
    unsigned int my_n_dofs = 0;
 
    switch (this->_elem->dim())
    {
    case 0:
      my_n_dofs = FE<0,family>::n_dofs(this->_elem, order);
      break;
    case 1:
      my_n_dofs = FE<1,family>::n_dofs(this->_elem, order);
      break;
    case 2:
      my_n_dofs = FE<2,family>::n_dofs(this->_elem, order);
      break;
    case 3:
      my_n_dofs = FE<3,family>::n_dofs(this->_elem, order);
      break;
    default:
      libmesh_error();
    }
 
    CPPUNIT_ASSERT_EQUAL(
      FEInterface::n_shape_functions(fe_type, this->_elem),
      my_n_dofs);
 
    CPPUNIT_ASSERT_EQUAL(
      FEInterface::get_continuity(fe_type),
      this->_fe->get_continuity());
 
    CPPUNIT_ASSERT_EQUAL(
      FEInterface::is_hierarchic(fe_type),
      this->_fe->is_hierarchic());
  }

References FETestBase< order, family, elem_type, 1, CaseName >::_elem, FETestBase< order, family, elem_type, 1, CaseName >::_fe, FETestBase< order, family, elem_type, 1, CaseName >::_sys, libMesh::Elem::dim(), libMesh::FEInterface::get_continuity(), libMesh::FEInterface::is_hierarchic(), libMesh::FE< Dim, T >::n_dofs(), libMesh::FEInterface::n_shape_functions(), and libMesh::System::variable_type().

testGradU()

template<Order order, FEFamily family, ElemType elem_type, typename CaseName >

void FETest< order, family, elem_type, CaseName >::testGradU ( )

inline

Definition at line 1013 of file fe_test.h.

  {
    LOG_UNIT_TEST;
    testLoop(checkGradU(this->_grad_tol));
  }

References FETestBase< order, family, elem_type, 1, CaseName >::_grad_tol, FETest< order, family, elem_type, CaseName >::checkGradU(), and FETest< order, family, elem_type, CaseName >::testLoop().

testGradUComp()

template<Order order, FEFamily family, ElemType elem_type, typename CaseName >

void FETest< order, family, elem_type, CaseName >::testGradUComp ( )

inline

Definition at line 1025 of file fe_test.h.

  {
    LOG_UNIT_TEST;
 
    auto f = [this](Point p)
      {
        Parameters dummy;
 
        Number grad_u_x = 0, grad_u_y = 0, grad_u_z = 0;
        for (std::size_t d = 0; d != this->_dof_indices.size(); ++d)
          {
            grad_u_x += this->_fe->get_dphidx()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
#if LIBMESH_DIM > 1
            grad_u_y += this->_fe->get_dphidy()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
#endif
#if LIBMESH_DIM > 2
            grad_u_z += this->_fe->get_dphidz()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
#endif
          }
 
        RealGradient true_grad = this->true_gradient(p);
 
        LIBMESH_ASSERT_NUMBERS_EQUAL(grad_u_x,
                                true_grad(0), this->_grad_tol);
        if (this->_dim > 1)
          LIBMESH_ASSERT_NUMBERS_EQUAL
            (grad_u_y, true_grad(1), this->_grad_tol);
        if (this->_dim > 2)
          LIBMESH_ASSERT_NUMBERS_EQUAL
            (grad_u_z, true_grad(2), this->_grad_tol);
      };
 
    testLoop(f);
  }

References FETestBase< order, family, elem_type, 1, CaseName >::_dim, FETestBase< order, family, elem_type, 1, CaseName >::_dof_indices, FETestBase< order, family, elem_type, 1, CaseName >::_fe, FETestBase< order, family, elem_type, 1, CaseName >::_grad_tol, FETestBase< order, family, elem_type, 1, CaseName >::_sys, libMesh::System::current_local_solution, FETest< order, family, elem_type, CaseName >::testLoop(), and FETestBase< order, family, elem_type, 1, CaseName >::true_gradient().

testGradUPastVertices()

template<Order order, FEFamily family, ElemType elem_type, typename CaseName >

void FETest< order, family, elem_type, CaseName >::testGradUPastVertices ( )

inline

Definition at line 1019 of file fe_test.h.

  {
    LOG_UNIT_TEST;
    testPastVerticesLoop(checkGradU(this->_past_vertex_tol_factor * this->_grad_tol));
  }

References FETestBase< order, family, elem_type, 1, CaseName >::_grad_tol, FETestBase< order, family, elem_type, 1, CaseName >::_past_vertex_tol_factor, FETest< order, family, elem_type, CaseName >::checkGradU(), and FETest< order, family, elem_type, CaseName >::testPastVerticesLoop().

testHessU()

template<Order order, FEFamily family, ElemType elem_type, typename CaseName >

void FETest< order, family, elem_type, CaseName >::testHessU ( )

inline

Definition at line 1061 of file fe_test.h.

  {
    LOG_UNIT_TEST;
 
    // Szabab elements don't have second derivatives yet
    if (family == SZABAB)
      return;
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
    auto f = [this](Point p)
      {
        Tensor hess_u;
        for (std::size_t d = 0; d != this->_dof_indices.size(); ++d)
          hess_u += this->_fe->get_d2phi()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
 
        // TODO: Yeah we'll test the ugly expressions later.
        if (family == RATIONAL_BERNSTEIN && order > 1)
          return;
 
        RealTensor true_hess = this->true_hessian(p);
 
        LIBMESH_ASSERT_NUMBERS_EQUAL
          (true_hess(0,0), hess_u(0,0), this->_hess_tol);
        if (this->_dim > 1)
          {
            LIBMESH_ASSERT_NUMBERS_EQUAL
              (hess_u(0,1), hess_u(1,0), this->_hess_tol);
            LIBMESH_ASSERT_NUMBERS_EQUAL
              (true_hess(0,1), hess_u(0,1), this->_hess_tol);
            LIBMESH_ASSERT_NUMBERS_EQUAL
              (true_hess(1,1), hess_u(1,1), this->_hess_tol);
          }
        if (this->_dim > 2)
          {
            LIBMESH_ASSERT_NUMBERS_EQUAL
              (hess_u(0,2), hess_u(2,0), this->_hess_tol);
            LIBMESH_ASSERT_NUMBERS_EQUAL
              (hess_u(1,2), hess_u(2,1), this->_hess_tol);
            LIBMESH_ASSERT_NUMBERS_EQUAL
              (true_hess(0,2), hess_u(0,2), this->_hess_tol);
            LIBMESH_ASSERT_NUMBERS_EQUAL
              (true_hess(1,2), hess_u(1,2), this->_hess_tol);
            LIBMESH_ASSERT_NUMBERS_EQUAL
              (true_hess(2,2), hess_u(2,2), this->_hess_tol);
          }
      };
 
    testLoop(f);
#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
  }

References FETestBase< order, family, elem_type, 1, CaseName >::_dim, FETestBase< order, family, elem_type, 1, CaseName >::_dof_indices, FETestBase< order, family, elem_type, 1, CaseName >::_fe, FETestBase< order, family, elem_type, 1, CaseName >::_hess_tol, FETestBase< order, family, elem_type, 1, CaseName >::_sys, libMesh::System::current_local_solution, libMesh::RATIONAL_BERNSTEIN, libMesh::SZABAB, FETest< order, family, elem_type, CaseName >::testLoop(), and FETestBase< order, family, elem_type, 1, CaseName >::true_hessian().

testHessUComp()

template<Order order, FEFamily family, ElemType elem_type, typename CaseName >

void FETest< order, family, elem_type, CaseName >::testHessUComp ( )

inline

Definition at line 1112 of file fe_test.h.

  {
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
    LOG_UNIT_TEST;
 
    // Szabab elements don't have second derivatives yet
    if (family == SZABAB)
      return;
 
    auto f = [this](Point p)
      {
        Number hess_u_xx = 0, hess_u_xy = 0, hess_u_yy = 0,
               hess_u_xz = 0, hess_u_yz = 0, hess_u_zz = 0;
        for (std::size_t d = 0; d != this->_dof_indices.size(); ++d)
          {
            hess_u_xx += this->_fe->get_d2phidx2()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
#if LIBMESH_DIM > 1
            hess_u_xy += this->_fe->get_d2phidxdy()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
            hess_u_yy += this->_fe->get_d2phidy2()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
#endif
#if LIBMESH_DIM > 2
            hess_u_xz += this->_fe->get_d2phidxdz()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
            hess_u_yz += this->_fe->get_d2phidydz()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
            hess_u_zz += this->_fe->get_d2phidz2()[d][0] * (*this->_sys->current_local_solution)(this->_dof_indices[d]);
#endif
          }
 
        // TODO: Yeah we'll test the ugly expressions later.
        if (family == RATIONAL_BERNSTEIN && order > 1)
          return;
 
        RealTensor true_hess = this->true_hessian(p);
 
        LIBMESH_ASSERT_NUMBERS_EQUAL
          (true_hess(0,0), hess_u_xx, this->_hess_tol);
        if (this->_dim > 1)
          {
            LIBMESH_ASSERT_NUMBERS_EQUAL
              (true_hess(0,1), hess_u_xy, this->_hess_tol);
            LIBMESH_ASSERT_NUMBERS_EQUAL
              (true_hess(1,1), hess_u_yy, this->_hess_tol);
          }
        if (this->_dim > 2)
          {
            LIBMESH_ASSERT_NUMBERS_EQUAL
              (true_hess(0,2), hess_u_xz, this->_hess_tol);
            LIBMESH_ASSERT_NUMBERS_EQUAL
              (true_hess(1,2), hess_u_yz, this->_hess_tol);
            LIBMESH_ASSERT_NUMBERS_EQUAL
              (true_hess(2,2), hess_u_zz, this->_hess_tol);
          }
      };
 
    testLoop(f);
#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
  }

References FETestBase< order, family, elem_type, 1, CaseName >::_dim, FETestBase< order, family, elem_type, 1, CaseName >::_dof_indices, FETestBase< order, family, elem_type, 1, CaseName >::_fe, FETestBase< order, family, elem_type, 1, CaseName >::_hess_tol, FETestBase< order, family, elem_type, 1, CaseName >::_sys, libMesh::System::current_local_solution, libMesh::RATIONAL_BERNSTEIN, libMesh::SZABAB, FETest< order, family, elem_type, CaseName >::testLoop(), and FETestBase< order, family, elem_type, 1, CaseName >::true_hessian().

testLoop()

template<Order order, FEFamily family, ElemType elem_type, typename CaseName >

template<typename Functor >

void FETest< order, family, elem_type, CaseName >::testLoop ( Functor  f )

inline

Definition at line 660 of file fe_test.h.

  {
    // Handle the "more processors than elements" case
    if (!this->_elem)
      return;
 
    // These tests require exceptions to be enabled because a
    // TypeTensor::solve() call down in Elem::contains_point()
    // actually throws a non-fatal exception for a certain Point which
    // is not in the Elem. When exceptions are not enabled, this test
    // simply aborts.
#ifdef LIBMESH_ENABLE_EXCEPTIONS
    for (unsigned int i=0; i != this->_nx; ++i)
      for (unsigned int j=0; j != this->_ny; ++j)
        for (unsigned int k=0; k != this->_nz; ++k)
          {
            Point p = Real(i)/this->_nx;
            if (j > 0)
              p(1) = Real(j)/this->_ny;
            if (k > 0)
              p(2) = Real(k)/this->_nz;
            if (!this->_elem->contains_point(p))
              continue;
 
            // If at a singular node, cannot use FEMap::map
            if (this->_elem->local_singular_node(p) != invalid_uint)
              continue;
 
            std::vector<Point> master_points
              (1, FEMap::inverse_map(this->_dim, this->_elem, p));
 
            // Reinit at point to test against analytic solution
            this->_fe->reinit(this->_elem, &master_points);
 
            f(p);
          }
#endif // LIBMESH_ENABLE_EXCEPTIONS
  }

References FETestBase< order, family, elem_type, 1, CaseName >::_dim, FETestBase< order, family, elem_type, 1, CaseName >::_elem, FETestBase< order, family, elem_type, 1, CaseName >::_fe, FETestBase< order, family, elem_type, 1, CaseName >::_nx, FETestBase< order, family, elem_type, 1, CaseName >::_ny, FETestBase< order, family, elem_type, 1, CaseName >::_nz, libMesh::Elem::contains_point(), libMesh::invalid_uint, libMesh::FEMap::inverse_map(), libMesh::Elem::local_singular_node(), and libMesh::Real.

Referenced by FETest< order, family, elem_type, CaseName >::testGradU(), FETest< order, family, elem_type, CaseName >::testGradUComp(), FETest< order, family, elem_type, CaseName >::testHessU(), FETest< order, family, elem_type, CaseName >::testHessUComp(), and FETest< order, family, elem_type, CaseName >::testU().

testNodalSoln()

template<Order order, FEFamily family, ElemType elem_type, typename CaseName >

void FETest< order, family, elem_type, CaseName >::testNodalSoln ( )

inline

Definition at line 770 of file fe_test.h.

    {
      LOG_UNIT_TEST;
 
      if (!this->_elem)
        return;
 
      // The Lagrange nodal solution is the interpolant evaluated at
      // each node, including nodes which carry no degree of freedom
      if (family != LAGRANGE && family != L2_LAGRANGE)
        return;
 
      if (this->_elem->runtime_topology())
        return;
 
      const FEType fe_type = this->_sys->variable_type(0);
      const unsigned int n_dofs = FEInterface::n_dofs(fe_type, this->_elem);
 
      std::vector<Number> elem_soln(n_dofs, 0), nodal_soln;
      for (auto i : make_range(n_dofs))
        {
          elem_soln[i] = 1;
          FEInterface::nodal_soln(this->_dim, fe_type, this->_elem,
                                  elem_soln, nodal_soln);
          elem_soln[i] = 0;
 
          CPPUNIT_ASSERT_EQUAL(std::size_t(this->_elem->n_nodes()),
                               nodal_soln.size());
          for (auto n : this->_elem->node_index_range())
            LIBMESH_ASSERT_NUMBERS_EQUAL
              (FEInterface::shape(fe_type, this->_elem, i,
                                  this->_elem->master_point(n)),
               nodal_soln[n], TOLERANCE*TOLERANCE);
        }
    }

References FETestBase< order, family, elem_type, 1, CaseName >::_dim, FETestBase< order, family, elem_type, 1, CaseName >::_elem, FETestBase< order, family, elem_type, 1, CaseName >::_sys, libMesh::L2_LAGRANGE, libMesh::LAGRANGE, libMesh::make_range(), libMesh::Elem::master_point(), libMesh::FEInterface::n_dofs(), libMesh::Elem::n_nodes(), libMesh::FEInterface::nodal_soln(), libMesh::Elem::node_index_range(), libMesh::Elem::runtime_topology(), libMesh::FEInterface::shape(), libMesh::TOLERANCE, and libMesh::System::variable_type().

testPartitionOfUnity()

template<Order order, FEFamily family, ElemType elem_type, typename CaseName >

void FETest< order, family, elem_type, CaseName >::testPartitionOfUnity ( )

inline

Definition at line 806 of file fe_test.h.

    {
      LOG_UNIT_TEST;
 
      if (!this->_elem)
        return;
 
      this->_fe->reinit(this->_elem);
 
      bool satisfies_partition_of_unity = true;
      for (const auto qp : make_range(this->_qrule->n_points()))
      {
        Real phi_sum = 0;
        for (std::size_t d = 0; d != this->_dof_indices.size(); ++d)
          phi_sum += this->_fe->get_phi()[d][qp];
        if (phi_sum < (1 - TOLERANCE) || phi_sum > (1 + TOLERANCE))
        {
          satisfies_partition_of_unity = false;
          break;
        }
      }
 
      switch (this->_fe->get_family())
      {
        case MONOMIAL:
        {
          switch (this->_fe->get_order())
          {
            case CONSTANT:
              CPPUNIT_ASSERT(satisfies_partition_of_unity);
              break;
 
            default:
              CPPUNIT_ASSERT(!satisfies_partition_of_unity);
              break;
          }
          break;
        }
        case SZABAB:
        case HIERARCHIC:
        case L2_HIERARCHIC:
        {
          switch (this->_fe->get_order())
          {
            case FIRST:
              CPPUNIT_ASSERT(satisfies_partition_of_unity);
              break;
 
            default:
              CPPUNIT_ASSERT(!satisfies_partition_of_unity);
              break;
          }
          break;
        }
 
        case XYZ:
        case CLOUGH:
        case HERMITE:
        {
          CPPUNIT_ASSERT(!satisfies_partition_of_unity);
          break;
        }
 
        case LAGRANGE:
        case L2_LAGRANGE:
        case L2_LAGRANGE_GLL:
        case BERNSTEIN:
        case RATIONAL_BERNSTEIN:
        {
          CPPUNIT_ASSERT(satisfies_partition_of_unity);
          break;
        }
 
        default:
          CPPUNIT_FAIL("Uncovered FEFamily");
      }
    }

References FETestBase< order, family, elem_type, 1, CaseName >::_dof_indices, FETestBase< order, family, elem_type, 1, CaseName >::_elem, FETestBase< order, family, elem_type, 1, CaseName >::_fe, FETestBase< order, family, elem_type, 1, CaseName >::_qrule, libMesh::BERNSTEIN, libMesh::CLOUGH, libMesh::CONSTANT, libMesh::FIRST, libMesh::HERMITE, libMesh::HIERARCHIC, libMesh::L2_HIERARCHIC, libMesh::L2_LAGRANGE, libMesh::L2_LAGRANGE_GLL, libMesh::LAGRANGE, libMesh::make_range(), libMesh::MONOMIAL, libMesh::RATIONAL_BERNSTEIN, libMesh::Real, libMesh::SZABAB, libMesh::TOLERANCE, and libMesh::XYZ.

testPastVerticesLoop()

template<Order order, FEFamily family, ElemType elem_type, typename CaseName >

template<typename Functor >

void FETest< order, family, elem_type, CaseName >::testPastVerticesLoop ( Functor  f )

inline

Definition at line 707 of file fe_test.h.

  {
    // Handle the "more processors than elements" case
    if (!this->_elem)
      return;
 
    // We don't have reference elements for these
    if (this->_elem->infinite() || this->_elem->runtime_topology())
      return;
 
    const Real delta = this->_past_vertex_offset;
 
    // An edge element is its own edge
    std::vector<std::pair<unsigned int, unsigned int>> edges;
    if (this->_dim == 1)
      edges.emplace_back(0, 1);
    for (const auto e : this->_elem->edge_index_range())
      edges.emplace_back(this->_elem->local_edge_node(e, 0),
                         this->_elem->local_edge_node(e, 1));
 
    for (const auto v : make_range(this->_elem->n_vertices()))
      {
        // Shape functions are singular at, and on a plane through, a singular node
        if (this->_elem->is_singular_node(v))
          continue;
 
        for (const auto & [a, b] : edges)
          for (const Real sign : {Real(-1), Real(1)})
            {
              const Point master_p = this->_elem->master_point(v) + sign * delta *
                (this->_elem->master_point(b) - this->_elem->master_point(a));
 
              std::vector<Point> master_points(1, master_p);
              this->_fe->reinit(this->_elem, &master_points);
 
              f(FEMap::map(this->_dim, this->_elem, master_p));
            }
      }
  }

References FETestBase< order, family, elem_type, 1, CaseName >::_dim, FETestBase< order, family, elem_type, 1, CaseName >::_elem, FETestBase< order, family, elem_type, 1, CaseName >::_fe, FETestBase< order, family, elem_type, 1, CaseName >::_past_vertex_offset, b, libMesh::Elem::edge_index_range(), libMesh::Elem::infinite(), libMesh::Elem::is_singular_node(), libMesh::Elem::local_edge_node(), libMesh::make_range(), libMesh::FEMap::map(), libMesh::Elem::master_point(), libMesh::Elem::n_vertices(), and libMesh::Real.

Referenced by FETest< order, family, elem_type, CaseName >::testGradUPastVertices(), and FETest< order, family, elem_type, CaseName >::testUPastVertices().

testRefspaceNodes()

template<Order order, FEFamily family, ElemType elem_type, typename CaseName >

void FETest< order, family, elem_type, CaseName >::testRefspaceNodes ( )

inline

Definition at line 747 of file fe_test.h.

    {
      LOG_UNIT_TEST;
 
      if (!this->_elem)
        return;
 
      // We don't currently support reference infinite elements
      if (this->_elem->infinite())
        return;
 
      // We'll probably never support reference polytopes
      if (this->_elem->runtime_topology())
        return;
 
      std::vector<Point> nodes;
      this->_fe->get_refspace_nodes(this->_elem->type(), nodes);
      for (auto n : index_range(nodes))
        LIBMESH_ASSERT_REALVEC_EQUAL(nodes[n],
                                     this->_elem->master_point(n),
                                     TOLERANCE*TOLERANCE);
    }

References FETestBase< order, family, elem_type, 1, CaseName >::_elem, FETestBase< order, family, elem_type, 1, CaseName >::_fe, libMesh::index_range(), libMesh::Elem::infinite(), libMesh::Elem::master_point(), libMesh::Elem::runtime_topology(), libMesh::TOLERANCE, and libMesh::Elem::type().

testU()

template<Order order, FEFamily family, ElemType elem_type, typename CaseName >

void FETest< order, family, elem_type, CaseName >::testU ( )

inline

Definition at line 960 of file fe_test.h.

  {
    LOG_UNIT_TEST;
    testLoop(checkU(this->_value_tol));
  }

References FETestBase< order, family, elem_type, 1, CaseName >::_value_tol, FETest< order, family, elem_type, CaseName >::checkU(), and FETest< order, family, elem_type, CaseName >::testLoop().

testUPastVertices()

template<Order order, FEFamily family, ElemType elem_type, typename CaseName >

void FETest< order, family, elem_type, CaseName >::testUPastVertices ( )

inline

Definition at line 966 of file fe_test.h.

  {
    LOG_UNIT_TEST;
    testPastVerticesLoop(checkU(this->_past_vertex_tol_factor * this->_value_tol));
  }

References FETestBase< order, family, elem_type, 1, CaseName >::_past_vertex_tol_factor, FETestBase< order, family, elem_type, 1, CaseName >::_value_tol, FETest< order, family, elem_type, CaseName >::checkU(), and FETest< order, family, elem_type, CaseName >::testPastVerticesLoop().

true_gradient()

static RealGradient FETestBase< order, family, elem_type, build_nx, CaseName >::true_gradient ( Point  p )

inlinestaticprotectedinherited

Definition at line 300 of file fe_test.h.

  {
    Parameters dummy;
 
    Gradient true_grad;
    RealGradient returnval;
 
    if (family == RATIONAL_BERNSTEIN && order > 1)
      true_grad = rational_test_grad(p, dummy, "", "");
    else if (order > 3)
      true_grad = fe_quartic_test_grad(p, dummy, "", "");
    else if (FE_CAN_TEST_CUBIC)
      true_grad = fe_cubic_test_grad(p, dummy, "", "");
    else if (order > 1)
      {
        const Real & x = p(0);
        const Real & y = (LIBMESH_DIM > 1) ? p(1) : 0;
        const Real & z = (LIBMESH_DIM > 2) ? p(2) : 0;
 
        true_grad = Gradient(2*x+0.125*y+0.0625*z,
                             y+0.125*x+0.03125*z,
                             0.5*z+0.0625*x+0.03125*y);
      }
    else
      true_grad = Gradient(1.0, 0.25, 0.0625);
 
    for (unsigned int d=0; d != LIBMESH_DIM; ++d)
      {
        CPPUNIT_ASSERT(true_grad(d) ==
                       Number(libmesh_real(true_grad(d))));
 
        returnval(d) = libmesh_real(true_grad(d));
      }
 
    return returnval;
  }

true_hessian()

static RealTensor FETestBase< order, family, elem_type, build_nx, CaseName >::true_hessian ( Point  p )

inlinestaticprotectedinherited

Definition at line 338 of file fe_test.h.

  {
    const Real & x = p(0);
    const Real & y = LIBMESH_DIM > 1 ? p(1) : 0;
    const Real & z = LIBMESH_DIM > 2 ? p(2) : 0;
 
    if (order > 3)
      return RealTensor
        { 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),
                 -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),
                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 };
    else if (FE_CAN_TEST_CUBIC)
      return RealTensor
        { 6*x-4+2*(1-y), -2*x+z-1,     y-1,
               -2*x+z-1,     -2*z, x+1-2*y,
                    y-1,  x+1-2*y, 6*z-4 };
    else if (order > 1)
      return RealTensor
        { 2, 0.125, 0.0625,
          0.125, 1, 0.03125,
          0.0625, 0.03125, 0.5 };
 
    return RealTensor
      { 0, 0, 0,
        0, 0, 0,
        0, 0, 0 };
  }

Member Data Documentation

_case_name

std::string FETestBase< order, family, elem_type, build_nx, CaseName >::_case_name

protectedinherited

Name of the case, can be used in the test to switch cases.

Definition at line 280 of file fe_test.h.

_dim

unsigned int FETestBase< order, family, elem_type, build_nx, CaseName >::_dim

protectedinherited

Definition at line 282 of file fe_test.h.

_dof_indices

std::vector<dof_id_type> FETestBase< order, family, elem_type, build_nx, CaseName >::_dof_indices

protectedinherited

Definition at line 284 of file fe_test.h.

_elem

Elem* FETestBase< order, family, elem_type, build_nx, CaseName >::_elem

protectedinherited

Definition at line 283 of file fe_test.h.

_es

std::unique_ptr<EquationSystems> FETestBase< order, family, elem_type, build_nx, CaseName >::_es

protectedinherited

Definition at line 287 of file fe_test.h.

_fe

std::unique_ptr<FEBase> FETestBase< order, family, elem_type, build_nx, CaseName >::_fe

protectedinherited

Definition at line 288 of file fe_test.h.

_grad_tol

Real FETestBase< order, family, elem_type, build_nx, CaseName >::_grad_tol

protectedinherited

Definition at line 291 of file fe_test.h.

_hess_tol

Real FETestBase< order, family, elem_type, build_nx, CaseName >::_hess_tol

protectedinherited

Definition at line 291 of file fe_test.h.

_mesh

std::unique_ptr<Mesh> FETestBase< order, family, elem_type, build_nx, CaseName >::_mesh

protectedinherited

Definition at line 286 of file fe_test.h.

_nx

unsigned int FETestBase< order, family, elem_type, build_nx, CaseName >::_nx

protectedinherited

Definition at line 282 of file fe_test.h.

_ny

unsigned int FETestBase< order, family, elem_type, build_nx, CaseName >::_ny

protectedinherited

Definition at line 282 of file fe_test.h.

_nz

unsigned int FETestBase< order, family, elem_type, build_nx, CaseName >::_nz

protectedinherited

Definition at line 282 of file fe_test.h.

_past_vertex_offset

Real FETestBase< order, family, elem_type, build_nx, CaseName >::_past_vertex_offset

protectedinherited

The offset, as a fraction of an edge, of the points past the vertices.

Definition at line 294 of file fe_test.h.

_past_vertex_tol_factor

Real FETestBase< order, family, elem_type, build_nx, CaseName >::_past_vertex_tol_factor

protectedinherited

The factor by which the tolerances grow at the points past the vertices.

Definition at line 297 of file fe_test.h.

_qrule

std::unique_ptr<QGauss> FETestBase< order, family, elem_type, build_nx, CaseName >::_qrule

protectedinherited

Definition at line 289 of file fe_test.h.

_sys

System* FETestBase< order, family, elem_type, build_nx, CaseName >::_sys

protectedinherited

Definition at line 285 of file fe_test.h.

_value_tol

Real FETestBase< order, family, elem_type, build_nx, CaseName >::_value_tol

protectedinherited

Definition at line 291 of file fe_test.h.

libmesh_suite_name

std::string FETestBase< order, family, elem_type, build_nx, CaseName >::libmesh_suite_name

protectedinherited

Definition at line 277 of file fe_test.h.

---

The documentation for this class was generated from the following file:

• tests/fe/fe_test.h