#include <MoFEM.hpp>

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

Examples: operators_tests.cpp.

Definition at line 34 of file operators_tests.cpp.

◆ OpConvectiveTermLhsDu

template<int FIELD_DIM>

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

Examples: operators_tests.cpp, and shallow_wave.cpp.

Definition at line 56 of file operators_tests.cpp.

◆ OpConvectiveTermLhsDy

template<int FIELD_DIM>

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

Examples: operators_tests.cpp, and shallow_wave.cpp.

Definition at line 60 of file operators_tests.cpp.

◆ OpConvectiveTermRhs

template<int FIELD_DIM>

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

Examples: operators_tests.cpp, and shallow_wave.cpp.

Definition at line 52 of file operators_tests.cpp.

◆ OpGradGrad

template<int FIELD_DIM>

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

Definition at line 44 of file operators_tests.cpp.

◆ OpGradTimesTensor

template<int FIELD_DIM>

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

Examples: heat_equation.cpp, operators_tests.cpp, photon_diffusion.cpp, poisson_2d_homogeneous.cpp, seepage.cpp, tensor_divergence_operator.cpp, and wave_equation.cpp.

Definition at line 48 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[]
	)

Examples: operators_tests.cpp.

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;
 }