FETestBase< order, family, elem_type, build_nx > Class Template Reference

#include <fe_test.h>

Inheritance diagram for FETestBase< order, family, elem_type, build_nx >:

Public Member Functions
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
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

Detailed Description

template<Order order, FEFamily family, ElemType elem_type, unsigned int build_nx>
class FETestBase< order, family, elem_type, build_nx >

Definition at line 266 of file fe_test.h.

Member Function Documentation

setUp()

template<Order order, FEFamily family, ElemType elem_type, unsigned int build_nx>

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

inline

Definition at line 350 of file fe_test.h.

Referenced by FESideTest< order, family, elem_type >::setUp().

  {
    _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);

        const std::vector<std::vector<unsigned int>> 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

        // 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<Elem> polyhedron = std::make_unique<C0Polyhedron>(sides);
        _mesh->add_elem(std::move(polyhedron));
        _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())
          {
            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)
              elem->node_ref(i).set_extra_datum<Real>(weight_index, 1.);
            for (unsigned int i = nv; i != nve; ++i)
              elem->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)
              elem->node_ref(i).set_extra_datum<Real>(weight_index, w2);
            const Real w3 = rational_w * w2;
            for (unsigned int i = nvef; i != nn; ++i)
              elem->node_ref(i).set_extra_datum<Real>(weight_index, w3);
          }
      }

    _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?!?

    // 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()

template<Order order, FEFamily family, ElemType elem_type, unsigned int build_nx>

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

inline

Definition at line 597 of file fe_test.h.

Referenced by FESideTest< order, family, elem_type >::tearDown().

{}

true_gradient()

template<Order order, FEFamily family, ElemType elem_type, unsigned int build_nx>

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

inlinestaticprotected

Definition at line 283 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()

template<Order order, FEFamily family, ElemType elem_type, unsigned int build_nx>

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

inlinestaticprotected

Definition at line 321 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

_dim

template<Order order, FEFamily family, ElemType elem_type, unsigned int build_nx>

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

protected

Definition at line 271 of file fe_test.h.

_dof_indices

template<Order order, FEFamily family, ElemType elem_type, unsigned int build_nx>

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

protected

Definition at line 273 of file fe_test.h.

_elem

template<Order order, FEFamily family, ElemType elem_type, unsigned int build_nx>

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

protected

Definition at line 272 of file fe_test.h.

_es

template<Order order, FEFamily family, ElemType elem_type, unsigned int build_nx>

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

protected

Definition at line 276 of file fe_test.h.

_fe

template<Order order, FEFamily family, ElemType elem_type, unsigned int build_nx>

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

protected

Definition at line 277 of file fe_test.h.

_grad_tol

template<Order order, FEFamily family, ElemType elem_type, unsigned int build_nx>

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

protected

Definition at line 280 of file fe_test.h.

_hess_tol

template<Order order, FEFamily family, ElemType elem_type, unsigned int build_nx>

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

protected

Definition at line 280 of file fe_test.h.

_mesh

template<Order order, FEFamily family, ElemType elem_type, unsigned int build_nx>

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

protected

Definition at line 275 of file fe_test.h.

_nx

template<Order order, FEFamily family, ElemType elem_type, unsigned int build_nx>

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

protected

Definition at line 271 of file fe_test.h.

_ny

template<Order order, FEFamily family, ElemType elem_type, unsigned int build_nx>

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

protected

Definition at line 271 of file fe_test.h.

_nz

template<Order order, FEFamily family, ElemType elem_type, unsigned int build_nx>

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

protected

Definition at line 271 of file fe_test.h.

_qrule

template<Order order, FEFamily family, ElemType elem_type, unsigned int build_nx>

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

protected

Definition at line 278 of file fe_test.h.

_sys

template<Order order, FEFamily family, ElemType elem_type, unsigned int build_nx>

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

protected

Definition at line 274 of file fe_test.h.

_value_tol

template<Order order, FEFamily family, ElemType elem_type, unsigned int build_nx>

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

protected

Definition at line 280 of file fe_test.h.

libmesh_suite_name

template<Order order, FEFamily family, ElemType elem_type, unsigned int build_nx>

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

protected

Definition at line 269 of file fe_test.h.

---

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

• tests/fe/fe_test.h