v0.14.0
operators_tests.cpp

Test operators in forms integrators

Date
2022-12-11

TODO: Add more operators.

/**
* @file operators_tests.cpp
* @example operators_tests.cpp
* @brief Test operators in forms integrators
* @date 2022-12-11
*
* @copyright Copyright (c) 2022
*
* TODO: Add more operators.
*
*/
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
template <int DIM> struct ElementsAndOps {};
template <> struct ElementsAndOps<2> {
};
template <> struct ElementsAndOps<3> {
};
//! [Define dimension]
constexpr int SPACE_DIM = 3;
constexpr AssemblyType A = AssemblyType::PETSC; //< selected assembly type
constexpr IntegrationType I =
IntegrationType::GAUSS; //< selected integration type
// Specializations for tested operators
template <int FIELD_DIM>
template <int FIELD_DIM>
template <int FIELD_DIM>
template <int FIELD_DIM>
template <int FIELD_DIM>
constexpr bool debug = false;
int main(int argc, char *argv[]) {
// 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";
// Add logging channel for example
auto core_log = logging::core::get();
core_log->add_sink(
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,
// 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) {
auto post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(m_field);
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(
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);
};
// 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 = [&]() {
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(
pip->getOpDomainRhsPipeline().push_back(
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");
}
};
// test operator for convection (part of material time directives)
auto TestOpConvection = [&]() {
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(
"VECTOR", vec_mat));
"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");
}
};
CHKERR TestOpGradGrad();
CHKERR TestOpConvection();
}
// finish work cleaning memory, getting statistics, etc.
return 0;
}
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
OpGradGrad
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpGradGrad< 1, FIELD_DIM, SPACE_DIM > OpGradGrad
Definition: operators_tests.cpp:44
OpConvectiveTermRhs
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpConvectiveTermRhs< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermRhs
Definition: operators_tests.cpp:52
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
OpConvectiveTermLhsDu
FormsIntegrators< DomainEleOp >::Assembly< A >::TriLinearForm< I >::OpConvectiveTermLhsDu< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermLhsDu
Definition: operators_tests.cpp:56
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
FIELD_DIM
constexpr int FIELD_DIM
Definition: child_and_parent.cpp:15
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1294
I
constexpr IntegrationType I
Definition: operators_tests.cpp:31
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:595
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
main
int main(int argc, char *argv[])
Definition: operators_tests.cpp:64
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEM::PipelineManager::FaceEle
MoFEM::FaceElementForcesAndSourcesCore FaceEle
Definition: PipelineManager.hpp:35
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:310
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
OpConvectiveTermLhsDy
FormsIntegrators< DomainEleOp >::Assembly< A >::TriLinearForm< I >::OpConvectiveTermLhsDy< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermLhsDy
Definition: operators_tests.cpp:60
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
BiLinearForm
OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
Definition: operators_tests.cpp:48
MoFEM::IntegrationType
IntegrationType
Form integrator integration types.
Definition: FormsIntegrators.hpp:136
MoFEM::OperatorsTester
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
Definition: OperatorsTester.hpp:21
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1536
EntData
EntitiesFieldData::EntData EntData
Definition: operators_tests.cpp:34
ElementsAndOps
Definition: child_and_parent.cpp:18
debug
constexpr bool debug
Definition: operators_tests.cpp:62
DomainEleOp
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
MoFEM::AssemblyType
AssemblyType
[Storage and set boundary conditions]
Definition: FormsIntegrators.hpp:104
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1102
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
MoFEM::SmartPetscObj< Vec >
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
SPACE_DIM
constexpr int SPACE_DIM
[Define dimension]
Definition: operators_tests.cpp:29
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
help
static char help[]
Definition: operators_tests.cpp:16