FETest< order, family, elem_type > Class Template Reference

#include <fe_test.h>

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

Public Member Functions
void	setUp ()
void	tearDown ()
void	testU ()
void	testGradU ()
void	testGradUComp ()

Private Attributes
unsigned int	_dim
unsigned int	_nx
unsigned int	_ny
unsigned int	_nz
Elem *	_elem
std::vector< dof_id_type >	_dof_indices
FEBase *	_fe
Mesh *	_mesh
System *	_sys
EquationSystems *	_es

Detailed Description

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

Definition at line 120 of file fe_test.h.

Member Function Documentation

setUp()

template<Order order, FEFamily family, ElemType elem_type>

void FETest< order, family, elem_type >::setUp ( )

inline

Definition at line 132 of file fe_test.h.

  {
    _mesh = new Mesh(*TestCommWorld);
    const std::unique_ptr<Elem> test_elem = Elem::build(elem_type);
    _dim = test_elem->dim();
    const unsigned int ny = _dim > 1;
    const unsigned int nz = _dim > 2;
 
    unsigned char weight_index = 0;
 
    if (family == RATIONAL_BERNSTEIN)
      {
        // Make sure we can handle non-zero weight indices
        _mesh->add_node_integer("buffer integer");
        weight_index = cast_int<unsigned char>
          (_mesh->add_node_datum<Real>("rational_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);
      }
 
    MeshTools::Generation::build_cube (*_mesh,
                                       1, ny, nz,
                                       0., 1., 0., ny, 0., nz,
                                       elem_type);
 
    // 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 = new EquationSystems(*_mesh);
    _sys = &(_es->add_system<System> ("SimpleSystem"));
    _sys->add_variable("u", order, family);
    _es->init();
 
    // Clough-Tocher elements still don't work multithreaded
    if (family == CLOUGH && libMesh::n_threads() > 1)
      return;
 
    if (family == RATIONAL_BERNSTEIN && order > 1)
      {
        _sys->project_solution(rational_test, rational_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).release();
    _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
 
    auto rng = _mesh->active_local_element_ptr_range();
    _elem = rng.begin() == rng.end() ? nullptr : *(rng.begin());
 
    _sys->get_dof_map().dof_indices(_elem, _dof_indices);
 
    _nx = 10;
    _ny = (_dim > 1) ? _nx : 1;
    _nz = (_dim > 2) ? _nx : 1;
  }

References libMesh::DistributedMesh::active_element_ptr_range(), libMesh::DistributedMesh::active_local_element_ptr_range(), libMesh::MeshBase::add_node_datum(), libMesh::MeshBase::add_node_integer(), libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), libMesh::FEGenericBase< OutputType >::build(), libMesh::Elem::build(), libMesh::MeshTools::Generation::build_cube(), libMesh::CLOUGH, libMesh::DofMap::dof_indices(), libMesh::EDGE3, libMesh::System::get_dof_map(), libMesh::FEGenericBase< OutputType >::get_dphi(), libMesh::FEGenericBase< OutputType >::get_dphidx(), libMesh::FEGenericBase< OutputType >::get_dphidy(), libMesh::FEGenericBase< OutputType >::get_dphidz(), libMesh::FEGenericBase< OutputType >::get_phi(), libMesh::EquationSystems::init(), linear_test(), linear_test_grad(), libMesh::n_threads(), libMesh::EquationSystems::parameters, libMesh::System::project_solution(), libMesh::QUAD9, libMesh::RATIONAL_BERNSTEIN, rational_test(), rational_test_grad(), rational_w, libMesh::Real, libMesh::MeshBase::set_default_mapping_data(), TestCommWorld, and libMesh::System::variable_type().

tearDown()

template<Order order, FEFamily family, ElemType elem_type>

void FETest< order, family, elem_type >::tearDown ( )

inline

Definition at line 230 of file fe_test.h.

  {
    delete _fe;
    delete _es;
    delete _mesh;
  }

testGradU()

template<Order order, FEFamily family, ElemType elem_type>

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

inline

Definition at line 291 of file fe_test.h.

  {
    // Clough-Tocher elements still don't work multithreaded
    if (family == CLOUGH && libMesh::n_threads() > 1)
      return;
 
    // Handle the "more processors than elements" case
    if (!_elem)
      return;
 
    Parameters dummy;
 
    // 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 != _nx; ++i)
      for (unsigned int j=0; j != _ny; ++j)
        for (unsigned int k=0; k != _nz; ++k)
          {
            Point p(Real(i)/_nx);
            if (j > 0)
              p(1) = Real(j)/_ny;
            if (k > 0)
              p(2) = Real(k)/_ny;
            if (!_elem->contains_point(p))
              continue;
 
            std::vector<Point> master_points
              (1, FEMap::inverse_map(_dim, _elem, p));
 
            _fe->reinit(_elem, &master_points);
 
            Gradient grad_u = 0;
            for (std::size_t d = 0; d != _dof_indices.size(); ++d)
              grad_u += _fe->get_dphi()[d][0] * (*_sys->current_local_solution)(_dof_indices[d]);
 
            if (family == RATIONAL_BERNSTEIN && order > 1)
              {
                const Gradient rat_grad =
                  rational_test_grad(p, dummy, "", "");
 
                LIBMESH_ASSERT_FP_EQUAL(libmesh_real(grad_u(0)),
                                        libmesh_real(rat_grad(0)),
                                        TOLERANCE*sqrt(TOLERANCE));
                if (_dim > 1)
                  LIBMESH_ASSERT_FP_EQUAL(libmesh_real(grad_u(1)),
                                          libmesh_real(rat_grad(1)),
                                          TOLERANCE*sqrt(TOLERANCE));
                if (_dim > 2)
                  LIBMESH_ASSERT_FP_EQUAL(libmesh_real(grad_u(2)),
                                          libmesh_real(rat_grad(2)),
                                          TOLERANCE*sqrt(TOLERANCE));
              }
            else
              {
                LIBMESH_ASSERT_FP_EQUAL(libmesh_real(grad_u(0)), 1.0,
                                        TOLERANCE*sqrt(TOLERANCE));
                if (_dim > 1)
                  LIBMESH_ASSERT_FP_EQUAL(libmesh_real(grad_u(1)), 0.25,
                                          TOLERANCE*sqrt(TOLERANCE));
                if (_dim > 2)
                  LIBMESH_ASSERT_FP_EQUAL(libmesh_real(grad_u(2)), 0.0625,
                                          TOLERANCE*sqrt(TOLERANCE));
              }
          }
#endif
  }

References libMesh::CLOUGH, libMesh::Elem::contains_point(), libMesh::System::current_local_solution, libMesh::FEGenericBase< OutputType >::get_dphi(), libMesh::FEMap::inverse_map(), libMesh::libmesh_real(), libMesh::n_threads(), libMesh::RATIONAL_BERNSTEIN, rational_test_grad(), libMesh::Real, libMesh::FEAbstract::reinit(), std::sqrt(), and libMesh::TOLERANCE.

testGradUComp()

template<Order order, FEFamily family, ElemType elem_type>

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

inline

Definition at line 362 of file fe_test.h.

  {
    // Clough-Tocher elements still don't work multithreaded
    if (family == CLOUGH && libMesh::n_threads() > 1)
      return;
 
    // Handle the "more processors than elements" case
    if (!_elem)
      return;
 
    Parameters dummy;
 
    // 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 != _nx; ++i)
      for (unsigned int j=0; j != _ny; ++j)
        for (unsigned int k=0; k != _nz; ++k)
          {
            Point p(Real(i)/_nx);
            if (j > 0)
              p(1) = Real(j)/_ny;
            if (k > 0)
              p(2) = Real(k)/_ny;
            if (!_elem->contains_point(p))
              continue;
 
            std::vector<Point> master_points
              (1, FEMap::inverse_map(_dim, _elem, p));
 
            _fe->reinit(_elem, &master_points);
 
            Number grad_u_x = 0, grad_u_y = 0, grad_u_z = 0;
            for (std::size_t d = 0; d != _dof_indices.size(); ++d)
              {
                grad_u_x += _fe->get_dphidx()[d][0] * (*_sys->current_local_solution)(_dof_indices[d]);
#if LIBMESH_DIM > 1
                grad_u_y += _fe->get_dphidy()[d][0] * (*_sys->current_local_solution)(_dof_indices[d]);
#endif
#if LIBMESH_DIM > 2
                grad_u_z += _fe->get_dphidz()[d][0] * (*_sys->current_local_solution)(_dof_indices[d]);
#endif
              }
 
            if (family == RATIONAL_BERNSTEIN && order > 1)
              {
                const Gradient rat_grad =
                  rational_test_grad(p, dummy, "", "");
 
                LIBMESH_ASSERT_FP_EQUAL(libmesh_real(grad_u_x),
                                        libmesh_real(rat_grad(0)),
                                        TOLERANCE*sqrt(TOLERANCE));
                if (_dim > 1)
                  LIBMESH_ASSERT_FP_EQUAL(libmesh_real(grad_u_y),
                                          libmesh_real(rat_grad(1)),
                                          TOLERANCE*sqrt(TOLERANCE));
                if (_dim > 2)
                  LIBMESH_ASSERT_FP_EQUAL(libmesh_real(grad_u_z),
                                          libmesh_real(rat_grad(2)),
                                          TOLERANCE*sqrt(TOLERANCE));
              }
            else
              {
                LIBMESH_ASSERT_FP_EQUAL(libmesh_real(grad_u_x), 1.0,
                                        TOLERANCE*sqrt(TOLERANCE));
                if (_dim > 1)
                  LIBMESH_ASSERT_FP_EQUAL(libmesh_real(grad_u_y), 0.25,
                                          TOLERANCE*sqrt(TOLERANCE));
                if (_dim > 2)
                  LIBMESH_ASSERT_FP_EQUAL(libmesh_real(grad_u_z), 0.0625,
                                          TOLERANCE*sqrt(TOLERANCE));
              }
          }
#endif
  }

References libMesh::CLOUGH, libMesh::Elem::contains_point(), libMesh::System::current_local_solution, libMesh::FEGenericBase< OutputType >::get_dphidx(), libMesh::FEGenericBase< OutputType >::get_dphidy(), libMesh::FEGenericBase< OutputType >::get_dphidz(), libMesh::FEMap::inverse_map(), libMesh::libmesh_real(), libMesh::n_threads(), libMesh::RATIONAL_BERNSTEIN, rational_test_grad(), libMesh::Real, libMesh::FEAbstract::reinit(), std::sqrt(), and libMesh::TOLERANCE.

testU()

template<Order order, FEFamily family, ElemType elem_type>

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

inline

Definition at line 237 of file fe_test.h.

  {
    // Clough-Tocher elements still don't work multithreaded
    if (family == CLOUGH && libMesh::n_threads() > 1)
      return;
 
    // Handle the "more processors than elements" case
    if (!_elem)
      return;
 
    Parameters dummy;
 
    // 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 != _nx; ++i)
      for (unsigned int j=0; j != _ny; ++j)
        for (unsigned int k=0; k != _nz; ++k)
          {
            Real x = (Real(i)/_nx), y = 0, z = 0;
            Point p = x;
            if (j > 0)
              p(1) = y = (Real(j)/_ny);
            if (k > 0)
              p(2) = z = (Real(k)/_nz);
            if (!_elem->contains_point(p))
              continue;
 
            std::vector<Point> master_points
              (1, FEMap::inverse_map(_dim, _elem, p));
 
            _fe->reinit(_elem, &master_points);
 
            Number u = 0;
            for (std::size_t d = 0; d != _dof_indices.size(); ++d)
              u += _fe->get_phi()[d][0] * (*_sys->current_local_solution)(_dof_indices[d]);
 
            if (family == RATIONAL_BERNSTEIN && order > 1)
              LIBMESH_ASSERT_FP_EQUAL
                (libmesh_real(u),
                 libmesh_real(rational_test(p, dummy, "", "")),
                 TOLERANCE*TOLERANCE);
            else
              LIBMESH_ASSERT_FP_EQUAL
                (libmesh_real(u),
                 libmesh_real(x + 0.25*y + 0.0625*z),
                 TOLERANCE*TOLERANCE);
          }
#endif
  }

References libMesh::CLOUGH, libMesh::Elem::contains_point(), libMesh::System::current_local_solution, libMesh::FEGenericBase< OutputType >::get_phi(), libMesh::FEMap::inverse_map(), libMesh::libmesh_real(), libMesh::n_threads(), libMesh::RATIONAL_BERNSTEIN, rational_test(), libMesh::Real, libMesh::FEAbstract::reinit(), and libMesh::TOLERANCE.

Member Data Documentation

_dim

template<Order order, FEFamily family, ElemType elem_type>

unsigned int FETest< order, family, elem_type >::_dim

private

Definition at line 123 of file fe_test.h.

_dof_indices

template<Order order, FEFamily family, ElemType elem_type>

std::vector<dof_id_type> FETest< order, family, elem_type >::_dof_indices

private

Definition at line 125 of file fe_test.h.

_elem

template<Order order, FEFamily family, ElemType elem_type>

Elem* FETest< order, family, elem_type >::_elem

private

Definition at line 124 of file fe_test.h.

_es

template<Order order, FEFamily family, ElemType elem_type>

EquationSystems* FETest< order, family, elem_type >::_es

private

Definition at line 129 of file fe_test.h.

_fe

template<Order order, FEFamily family, ElemType elem_type>

FEBase* FETest< order, family, elem_type >::_fe

private

Definition at line 126 of file fe_test.h.

_mesh

template<Order order, FEFamily family, ElemType elem_type>

Mesh* FETest< order, family, elem_type >::_mesh

private

Definition at line 127 of file fe_test.h.

_nx

template<Order order, FEFamily family, ElemType elem_type>

unsigned int FETest< order, family, elem_type >::_nx

private

Definition at line 123 of file fe_test.h.

_ny

template<Order order, FEFamily family, ElemType elem_type>

unsigned int FETest< order, family, elem_type >::_ny

private

Definition at line 123 of file fe_test.h.

_nz

template<Order order, FEFamily family, ElemType elem_type>

unsigned int FETest< order, family, elem_type >::_nz

private

Definition at line 123 of file fe_test.h.

_sys

template<Order order, FEFamily family, ElemType elem_type>

System* FETest< order, family, elem_type >::_sys

private

Definition at line 128 of file fe_test.h.

---

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

• tests/fe/fe_test.h