dual_shape_verification_test.C

Go to the documentation of this file.

// libmesh includes
#include "libmesh/libmesh.h"
#include "libmesh/edge_edge2.h"
#include "libmesh/fe.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/mesh.h"

// unit test includes
#include "test_comm.h"

#include "libmesh_cppunit.h"

using namespace libMesh;

class DualShapeTest : public CppUnit::TestCase
{
public:
  LIBMESH_CPPUNIT_TEST_SUITE( DualShapeTest );
  CPPUNIT_TEST( testEdge2Lagrange );
  CPPUNIT_TEST_SUITE_END();

private:
  std::unique_ptr<Mesh> _mesh;
  std::unique_ptr<QGauss> _qrule;
  std::unique_ptr<FEBase> _fe;
  Elem * _elem;

public:

  void testEdge2Lagrange ()
  {
    LOG_UNIT_TEST;
    if (!_elem)
      return;

    _fe->reinit(_elem);

    const DenseMatrix<Real> & dual_coeff = _fe->get_dual_coeff();

    CPPUNIT_ASSERT_EQUAL(dual_coeff.m(), unsigned(2));
    CPPUNIT_ASSERT_EQUAL(dual_coeff.n(), unsigned(2));

    // TOLERANCE*TOLERANCE works with double but not float128
    Real my_tol = TOLERANCE*std::sqrt(TOLERANCE);

    LIBMESH_ASSERT_FP_EQUAL(2, dual_coeff(0,0), my_tol);
    LIBMESH_ASSERT_FP_EQUAL(-1, dual_coeff(0,1), my_tol);
    LIBMESH_ASSERT_FP_EQUAL(-1, dual_coeff(1,0), my_tol);
    LIBMESH_ASSERT_FP_EQUAL(2, dual_coeff(1,1), my_tol);

    const auto & dual_phi = _fe->get_dual_phi();

    CPPUNIT_ASSERT_EQUAL(std::size_t(2), dual_phi.size());

    const auto & qpoints = _qrule->get_points();

    CPPUNIT_ASSERT_EQUAL(qpoints.size(), dual_phi[0].size());

    for (auto qp : index_range(dual_phi[0]))
      LIBMESH_ASSERT_FP_EQUAL(1./2. * (1. - 3.*qpoints[qp](0)), dual_phi[0][qp],
        my_tol);

    CPPUNIT_ASSERT_EQUAL(qpoints.size(), dual_phi[1].size());

    for (auto qp : index_range(dual_phi[1]))
      LIBMESH_ASSERT_FP_EQUAL(1./2. * (1. + 3.*qpoints[qp](0)), dual_phi[1][qp],
        my_tol);
  }

  void setUp()
    {
      FEType fe_type(FIRST, LAGRANGE);
      _fe = FEBase::build(1, fe_type);
      _fe->get_phi();
      _fe->get_dual_phi();

      _mesh = std::make_unique<Mesh>(*TestCommWorld);

      MeshTools::Generation::build_line(*_mesh, 1, -1, 1, EDGE2);

      auto rng = _mesh->active_local_element_ptr_range();
      _elem = rng.begin() == rng.end() ? nullptr : *(rng.begin());

      _qrule = std::make_unique<QGauss>(1, fe_type.default_quadrature_order());
      _fe->attach_quadrature_rule(_qrule.get());
    }

  void tearDown() {}

};

CPPUNIT_TEST_SUITE_REGISTRATION( DualShapeTest );