operators_tests.cpp File Reference

#include <MoFEM.hpp>

Go to the source code of this file.

Classes
struct	ElementsAndOps< DIM >
struct	ElementsAndOps< 2 >
struct	ElementsAndOps< 3 >

Typedefs
using	EntData = EntitiesFieldData::EntData
using	DomainEle = ElementsAndOps< SPACE_DIM >::DomainEle
using	DomainEleOp = DomainEle::UserDataOperator
using	PostProcEle = PostProcBrokenMeshInMoab< DomainEle >
template<int FIELD_DIM>
using	OpGradGrad = FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpGradGrad< 1, FIELD_DIM, SPACE_DIM >
template<int FIELD_DIM>
using	OpGradTimesTensor = FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM >
template<int FIELD_DIM>
using	OpConvectiveTermRhs = FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpConvectiveTermRhs< 1, FIELD_DIM, SPACE_DIM >
template<int FIELD_DIM>
using	OpConvectiveTermLhsDu = FormsIntegrators< DomainEleOp >::Assembly< A >::TriLinearForm< I >::OpConvectiveTermLhsDu< 1, FIELD_DIM, SPACE_DIM >
template<int FIELD_DIM>
using	OpConvectiveTermLhsDy = FormsIntegrators< DomainEleOp >::Assembly< A >::TriLinearForm< I >::OpConvectiveTermLhsDy< 1, FIELD_DIM, SPACE_DIM >

Functions
int	main (int argc, char *argv[])

Variables
static char	help [] = "...\n\n"
constexpr int	SPACE_DIM = 3
	[Define dimension] More...
constexpr AssemblyType	A = AssemblyType::PETSC
constexpr IntegrationType	I
constexpr bool	debug = false

Typedef Documentation

DomainEle

using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle

Definition at line 35 of file operators_tests.cpp.

DomainEleOp

using DomainEleOp = DomainEle::UserDataOperator

Definition at line 36 of file operators_tests.cpp.

EntData

using EntData = EntitiesFieldData::EntData

Definition at line 34 of file operators_tests.cpp.

OpConvectiveTermLhsDu

template<int FIELD_DIM>

using OpConvectiveTermLhsDu = FormsIntegrators<DomainEleOp>::Assembly< A>::TriLinearForm<I>::OpConvectiveTermLhsDu<1, FIELD_DIM, SPACE_DIM>

Definition at line 55 of file operators_tests.cpp.

OpConvectiveTermLhsDy

template<int FIELD_DIM>

using OpConvectiveTermLhsDy = FormsIntegrators<DomainEleOp>::Assembly< A>::TriLinearForm<I>::OpConvectiveTermLhsDy<1, FIELD_DIM, SPACE_DIM>

Definition at line 59 of file operators_tests.cpp.

OpConvectiveTermRhs

template<int FIELD_DIM>

using OpConvectiveTermRhs = FormsIntegrators<DomainEleOp>::Assembly< A>::LinearForm<I>::OpConvectiveTermRhs<1, FIELD_DIM, SPACE_DIM>

Definition at line 51 of file operators_tests.cpp.

OpGradGrad

template<int FIELD_DIM>

using OpGradGrad = FormsIntegrators<DomainEleOp>::Assembly<A>::BiLinearForm< I>::OpGradGrad<1, FIELD_DIM, SPACE_DIM>

Examples

heat_method.cpp.

Definition at line 43 of file operators_tests.cpp.

OpGradTimesTensor

template<int FIELD_DIM>

using OpGradTimesTensor = FormsIntegrators<DomainEleOp>::Assembly< A>::LinearForm<I>::OpGradTimesTensor<1, FIELD_DIM, SPACE_DIM>

Definition at line 47 of file operators_tests.cpp.

PostProcEle

using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>

Definition at line 37 of file operators_tests.cpp.

Function Documentation

main()

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 64 of file operators_tests.cpp.

                                 {
 
  // initialize petsc
  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
 
  try {
 
    // Create MoAB database
    moab::Core moab_core;
    moab::Interface &moab = moab_core;
 
    // Create MoFEM database and link it to MoAB
    MoFEM::Core mofem_core(moab);
    MoFEM::Interface &m_field = mofem_core;
 
    // Register DM Manager
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
 
    // Add logging channel for example
    auto core_log = logging::core::get();
    core_log->add_sink(
        LogManager::createSink(LogManager::getStrmWorld(), "OpTester"));
    LogManager::setLog("OpTester");
    MOFEM_LOG_TAG("OpTester", "OpTester");
 
    // Simple interface
    auto simple = m_field.getInterface<Simple>();
 
    // get options from command line
    CHKERR simple->getOptions();
    // load mesh file
    CHKERR simple->loadFile();
 
    // Scalar fields and vector field is tested. Add more fields, i.e. vector
    // field if needed.
    CHKERR simple->addDomainField("SCALAR", H1, AINSWORTH_LEGENDRE_BASE, 1);
    CHKERR simple->addDomainField("VECTOR", H1, AINSWORTH_LEGENDRE_BASE,
                                  SPACE_DIM);
 
    // set fields order, i.e. for most first cases order is sufficient.
    CHKERR simple->setFieldOrder("SCALAR", 1);
    CHKERR simple->setFieldOrder("VECTOR", 1);
 
    // setup problem
    CHKERR simple->setUp();
 
    // get operators tester
    auto opt = m_field.getInterface<OperatorsTester>(); // get interface to
                                                        // OperatorsTester
    auto pip = m_field.getInterface<PipelineManager>(); // get interface to
                                                        // pipeline manager
 
    auto post_proc = [&](auto dm, auto f_res, auto out_name) {
      MoFEMFunctionBegin;
      auto post_proc_fe =
          boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(m_field);
 
      using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
      auto scal_vec = boost::make_shared<VectorDouble>();
      auto vec_mat = boost::make_shared<MatrixDouble>();
 
      post_proc_fe->getOpPtrVector().push_back(
          new OpCalculateScalarFieldValues("SCALAR", scal_vec, f_res));
      post_proc_fe->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("VECTOR", vec_mat,
                                                      f_res));
 
      post_proc_fe->getOpPtrVector().push_back(
 
          new OpPPMap(
 
              post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
              {{"SCALAR", scal_vec}},
 
              {{"VECTOR", vec_mat}},
 
              {}, {})
 
      );
 
      CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
                                      post_proc_fe);
      post_proc_fe->writeFile(out_name);
      MoFEMFunctionReturn(0);
    };
 
    // Note: I do not push OPs transferring base to physical coordinates. That
    // is not needed to test OPs, and not required, since we like to test only
    // operators.
 
    // Test grad-grad operator for scalar and vector field specialisation
    auto TestOpGradGrad = [&]() {
      MoFEMFunctionBegin;
      MOFEM_LOG("OpTester", Sev::verbose) << "TestOpGradGrad";
 
      pip->getOpDomainLhsPipeline().clear();
      pip->getOpDomainRhsPipeline().clear();
 
      pip->setDomainLhsIntegrationRule([](int, int, int o) { return 2 * o; });
      pip->setDomainLhsIntegrationRule([](int, int, int o) { return 2 * o; });
 
      pip->getOpDomainLhsPipeline().push_back(new OpGradGrad<1>(
          "SCALAR", "SCALAR", [](double, double, double) { return 1; }));
      pip->getOpDomainLhsPipeline().push_back(new OpGradGrad<SPACE_DIM>(
          "VECTOR", "VECTOR", [](double, double, double) { return 1; }));
 
      auto scl_mat = boost::make_shared<MatrixDouble>();
      auto vec_mat = boost::make_shared<MatrixDouble>();
 
      pip->getOpDomainRhsPipeline().push_back(
          new OpCalculateScalarFieldGradient<SPACE_DIM>("SCALAR", scl_mat));
      pip->getOpDomainRhsPipeline().push_back(
          new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>("VECTOR",
                                                                   vec_mat));
 
      pip->getOpDomainRhsPipeline().push_back(new OpGradTimesTensor<1>(
          "SCALAR", scl_mat, [](double, double, double) { return 1; }));
      pip->getOpDomainRhsPipeline().push_back(new OpGradTimesTensor<SPACE_DIM>(
          "VECTOR", vec_mat, [](double, double, double) { return 1; }));
 
      constexpr double eps = 1e-6;
 
      auto x = opt->setRandomFields(simple->getDM(),
                                    {{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
      auto diff_x = opt->setRandomFields(
          simple->getDM(), {{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
 
      auto diff_res = opt->checkCentralFiniteDifference(
          simple->getDM(), simple->getDomainFEName(), pip->getDomainRhsFE(),
          pip->getDomainLhsFE(), x, SmartPetscObj<Vec>(), SmartPetscObj<Vec>(),
          diff_x, 0, 1, eps);
 
      // Calculate norm of difference between directive calculated from finite
      // difference, and tangent matrix.
      double fnorm;
      CHKERR VecNorm(diff_res, NORM_2, &fnorm);
      MOFEM_LOG_C("OpTester", Sev::inform, "TestOpGradGrad %3.4e", fnorm);
 
      constexpr double err = 1e-9;
      if (fnorm > err)
        SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                "Norm of directional derivative too large");
 
      if (debug) {
        // Example how to plot direction in direction diff_x. If instead
        // directionalCentralFiniteDifference(...) diff_res is used, then error
        // on directive is plotted.
        CHKERR post_proc(simple->getDM(),
                         //
                         opt->directionalCentralFiniteDifference(
                             simple->getDM(), simple->getDomainFEName(),
                             pip->getDomainRhsFE(), x, SmartPetscObj<Vec>(),
                             SmartPetscObj<Vec>(), diff_res, 0, 1, eps),
                         //
                         "out_directional_directive_gradgrad.h5m");
      }
 
      MoFEMFunctionReturn(0);
    };
 
    // test operator for convection (part of material time directives)
    auto TestOpConvection = [&]() {
      MoFEMFunctionBegin;
      MOFEM_LOG("OpTester", Sev::verbose) << "TestOpConvection";
 
      pip->getOpDomainLhsPipeline().clear();
      pip->getOpDomainRhsPipeline().clear();
 
      pip->setDomainLhsIntegrationRule([](int, int, int o) { return 2 * o; });
      pip->setDomainLhsIntegrationRule([](int, int, int o) { return 2 * o; });
 
      auto evaluate_field = [&](auto &pipe) {
        auto scl_mat = boost::make_shared<MatrixDouble>();
        auto vec_mat = boost::make_shared<MatrixDouble>();
        auto dot_vec_vel = boost::make_shared<MatrixDouble>();
        pipe.push_back(
            new OpCalculateScalarFieldGradient<SPACE_DIM>("SCALAR", scl_mat));
        pipe.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
            "VECTOR", vec_mat));
        pipe.push_back(new OpCalculateVectorFieldValuesDot<SPACE_DIM>(
            "VECTOR", dot_vec_vel));
        return std::make_tuple(scl_mat, vec_mat, dot_vec_vel);
      };
 
      auto [rhs_scl_mat, rhs_vec_mat, rhs_dot_vec_vel] =
          evaluate_field(pip->getOpDomainRhsPipeline());
      pip->getOpDomainRhsPipeline().push_back(
          new OpConvectiveTermRhs<1>("SCALAR", rhs_dot_vec_vel, rhs_scl_mat));
      pip->getOpDomainRhsPipeline().push_back(
          new OpConvectiveTermRhs<SPACE_DIM>("VECTOR", rhs_dot_vec_vel,
                                             rhs_vec_mat));
 
      auto [lhs_scl_mat, lhs_vec_mat, lhs_dot_vec_vel] =
          evaluate_field(pip->getOpDomainLhsPipeline());
 
      // I create op, set scaling function to calculate time directive, and add
      // operator pinter to pipeline
      auto op_convective_term_lhs_du_scalar =
          new OpConvectiveTermLhsDu<1>("SCALAR", "VECTOR", lhs_scl_mat);
      op_convective_term_lhs_du_scalar->feScalingFun =
          [](const FEMethod *fe_ptr) { return fe_ptr->ts_a; };
      pip->getOpDomainLhsPipeline().push_back(op_convective_term_lhs_du_scalar);
      pip->getOpDomainLhsPipeline().push_back(
          new OpConvectiveTermLhsDy<1>("SCALAR", "SCALAR", lhs_dot_vec_vel));
 
      auto op_convective_term_lhs_du_vec =
          new OpConvectiveTermLhsDu<SPACE_DIM>("VECTOR", "VECTOR", lhs_vec_mat);
      op_convective_term_lhs_du_vec->feScalingFun = [](const FEMethod *fe_ptr) {
        return fe_ptr->ts_a;
      };
      pip->getOpDomainLhsPipeline().push_back(op_convective_term_lhs_du_vec);
      pip->getOpDomainLhsPipeline().push_back(
          new OpConvectiveTermLhsDy<SPACE_DIM>("VECTOR", "VECTOR",
                                               lhs_dot_vec_vel));
 
      constexpr double eps = 1e-6;
 
      auto x = opt->setRandomFields(simple->getDM(),
                                    {{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
      auto x_t = opt->setRandomFields(
          simple->getDM(), {{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
      auto diff_x = opt->setRandomFields(
          simple->getDM(), {{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
 
      auto diff_res = opt->checkCentralFiniteDifference(
          simple->getDM(), simple->getDomainFEName(), pip->getDomainRhsFE(),
          pip->getDomainLhsFE(), x, x_t, SmartPetscObj<Vec>(), diff_x, 0, 1,
          eps);
 
      // Calculate norm of difference between directive calculated from finite
      // difference, and tangent matrix.
      double fnorm;
      CHKERR VecNorm(diff_res, NORM_2, &fnorm);
      MOFEM_LOG_C("OpTester", Sev::inform, "TestOpConvection %3.4e", fnorm);
 
      constexpr double err = 1e-9;
      if (fnorm > err) {
        SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                "Norm of directional derivative too large");
      }
 
      if (debug) {
        // Example how to plot direction in direction diff_x. If instead
        // directionalCentralFiniteDifference(...) diff_res is used, then error
        // on directive is plotted.
        CHKERR post_proc(simple->getDM(),
                         //
                         opt->directionalCentralFiniteDifference(
                             simple->getDM(), simple->getDomainFEName(),
                             pip->getDomainRhsFE(), x, x_t,
                             SmartPetscObj<Vec>(), diff_res, 0, 1, eps),
                         //
                         "out_directional_directive_convection.h5m");
      }
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR TestOpGradGrad();
    CHKERR TestOpConvection();
  }
  CATCH_ERRORS;
 
  // finish work cleaning memory, getting statistics, etc.
  MoFEM::Core::Finalize();
 
  return 0;
}

Variable Documentation

A

constexpr AssemblyType A = AssemblyType::PETSC

constexpr

Examples

ContactOps.hpp, EshelbianPlasticity.cpp, PlasticOps.hpp, PlasticOpsGeneric.hpp, Remodeling.cpp, ThermoElasticOps.hpp, analytical_poisson_field_split.cpp, approx_sphere.cpp, build_large_problem.cpp, cell_forces.cpp, dm_create_subdm.cpp, eigen_elastic.cpp, free_surface.cpp, heat_equation.cpp, level_set.cpp, minimal_surface_area.cpp, operators_tests.cpp, photon_diffusion.cpp, plate.cpp, prism_polynomial_approximation.cpp, quad_polynomial_approximation.cpp, simple_elasticity.cpp, simple_interface.cpp, tensor_divergence_operator.cpp, test_jacobian_of_simple_contact_element.cpp, testing_jacobian_of_hook_element.cpp, testing_jacobian_of_hook_scaled_with_density_element.cpp, and thermo_elastic.cpp.

Definition at line 30 of file operators_tests.cpp.

debug

constexpr bool debug = false

constexpr

Definition at line 62 of file operators_tests.cpp.

help

char help[] = "...\n\n"

static

Definition at line 16 of file operators_tests.cpp.

I

constexpr IntegrationType I

constexpr

Initial value:

=
    IntegrationType::GAUSS

Examples

Remodeling.cpp, Remodeling.hpp, ThermoElasticOps.hpp, free_surface.cpp, level_set.cpp, operators_tests.cpp, and tensor_divergence_operator.cpp.

Definition at line 31 of file operators_tests.cpp.

SPACE_DIM

constexpr int SPACE_DIM = 3

constexpr

[Define dimension]

Definition at line 29 of file operators_tests.cpp.