Testing higher derivatives of base functions

* \example higher_derivatives.cpp
* Testing higher derivatives of base functions
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
constexpr char FIELD_NAME[] = "U";
constexpr int BASE_DIM = 1;
constexpr int FIELD_DIM = 1;
constexpr int SPACE_DIM = 2;
template <int DIM> struct ElementsAndOps {};
template <> struct ElementsAndOps<2> {
using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle; ///< Finite elenent type
using DomainEleOp =
DomainEle::UserDataOperator; ///< Finire element operator type
using EntData = EntitiesFieldData::EntData; ///< Data on entities
* @brief Function to approximate
auto fun = [](const double x, const double y, const double z) {
return x * x + y * y + x * y + pow(x, 3) + pow(y, 3) + pow(x, 4) + pow(y, 4);
* @brief Function derivative
auto diff_fun = [](const double x, const double y, const double z) {
2 * x + y + 3 * pow(x, 2) + 4 * pow(x, 3),
2 * y + x + 3 * pow(y, 2) + 4 * pow(y, 3)};
* @brief Function second derivative
auto diff2_fun = [](const double x, const double y, const double z) {
2 + 6 * x + 12 * pow(x, 2), 1.,
1., 2 + 6 * y + 12 * pow(y, 2)};
* @brief OPerator to integrate mass matrix for least square approximation
* @brief Operator to integrate the right hand side matrix for the problem
struct AtomTest {
AtomTest(MoFEM::Interface &m_field) : mField(m_field) {}
* @brief Collected data use d by operator to evaluate errors for the test
struct CommonData {
boost::shared_ptr<MatrixDouble> invJacPtr;
boost::shared_ptr<VectorDouble> approxVals;
boost::shared_ptr<MatrixDouble> approxGradVals;
boost::shared_ptr<MatrixDouble> approxHessianVals;
* @brief Operator to evaluate errors
struct OpError;
* @brief Operator to evaluate errors
struct AtomTest::OpError : public DomainEleOp {
boost::shared_ptr<CommonData> commonDataPtr;
OpError(boost::shared_ptr<CommonData> &common_data_ptr)
: DomainEleOp(FIELD_NAME, OPROW), commonDataPtr(common_data_ptr) {
std::fill(doEntities.begin(), doEntities.end(), false);
doEntities[MBVERTEX] = true;
const int nb_integration_pts = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight(); // ger integration weights
auto t_val = getFTensor0FromVec(*(
commonDataPtr->approxVals)); // get function approximation at gauss pts
auto t_grad_val = getFTensor1FromMat<SPACE_DIM>(
->approxGradVals)); // get gradient of approximation at gauss pts
auto t_hessian_val = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
*(commonDataPtr)->approxHessianVals); // get hessian of approximation
// at integration pts
auto t_inv_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
*(commonDataPtr->invJacPtr)); // get inverse of element jacobian
auto t_coords = getFTensor1CoordsAtGaussPts(); // get coordinates of
// integration points
// Indices used for tensor operations
const double volume = getMeasure(); // get finite element area
auto t_row_base = data.getFTensor0N();
auto t_diff_row_base = data.getFTensor1DiffN<2>();
std::array<double, 3> error = {0, 0,
0}; // array for storing operator errors
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = t_w * volume;
// Calculate errors
double diff = t_val - fun(t_coords(0), t_coords(1), t_coords(2));
error[0] += alpha * pow(diff, 2);
auto t_diff_grad = diff_fun(t_coords(0), t_coords(1), t_coords(2));
t_diff_grad(i) -= t_grad_val(i);
error[1] += alpha * t_diff_grad(i) *
t_diff_grad(i); // note push forward derivatives
MOFEM_LOG("SELF", Sev::noisy) << "t_hessian_val " << t_hessian_push;
// hessian expected to have symmetry
if (std::abs(t_hessian_val(0, 1) - t_hessian_val(1, 0)) >
std::numeric_limits<float>::epsilon()) {
MOFEM_LOG("SELF", Sev::error) << "t_hessian_val " << t_hessian_val;
"Hessian should be symmetric");
auto t_diff_hessian = diff2_fun(t_coords(0), t_coords(1), t_coords(2));
t_diff_hessian(i, j) -= t_hessian_val(i, j);
error[2] = t_diff_hessian(i, j) * t_diff_hessian(i, j);
// move data to next integration point
// assemble error vector
std::array<int, 3> index = {0, 1, 2};
CHKERR VecSetValues(commonDataPtr->L2Vec, 3, index.data(), error.data(),
//! [Run programme]
//! [Run programme]
//! [Read mesh]
CHKERR simpleInterface->getOptions();
CHKERR simpleInterface->loadFile();
//! [Read mesh]
//! [Set up problem]
// Add field
CHKERR simpleInterface->addDomainField(FIELD_NAME, H1,
constexpr int order = 4;
//! [Set up problem]
//! [Push operators to pipeline]
auto rule = [](int, int, int p) -> int { return 2 * p; };
CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
auto beta = [](const double, const double, const double) { return 1; };
//! [Push operators to pipeline]
//! [Solve]
MOFEM_LOG("WORLD", Sev::inform) << "Solve problem";
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSetUp(solver);
auto dm = simpleInterface->getDM();
auto D = smartCreateDMVector(dm);
CHKERR KSPSolve(solver, F, D);
//! [Check results]
auto rule = [](int, int, int p) -> int { return 2 * p; };
rule); // set integration rule
// create data structures for operator
auto common_data_ptr = boost::make_shared<CommonData>();
common_data_ptr->L2Vec = createSmartVectorMPI(
mField.get_comm(), (!mField.get_comm_rank()) ? 3 : 0, 3);
common_data_ptr->approxVals = boost::make_shared<VectorDouble>();
common_data_ptr->approxGradVals = boost::make_shared<MatrixDouble>();
common_data_ptr->approxHessianVals = boost::make_shared<MatrixDouble>();
// create data strutires for evaluation of higher order derivatives
auto base_mass = boost::make_shared<MatrixDouble>();
auto data_l2 = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
common_data_ptr->invJacPtr = inv_jac_ptr;
// calculate jacobian at integration points
new OpCalculateHOJacForFace(jac_ptr));
// calculate incerse of jacobian at integration points
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
// calculate mass matrix to project derivatives
new OpBaseDerivativesMass<BASE_DIM>(base_mass, data_l2,
// calculate second derivative of base functions, i.e. hessian
new OpBaseDerivativesNext<BASE_DIM>(BaseDerivatives::SecondDerivative,
base_mass, data_l2,
// calculate third derivative
new OpBaseDerivativesNext<BASE_DIM>(BaseDerivatives::ThirdDerivative,
base_mass, data_l2,
// calculate forth derivative
new OpBaseDerivativesNext<BASE_DIM>(BaseDerivatives::ForthDerivative,
base_mass, data_l2,
// push first base directives tp physical element shape
new OpSetInvJacH1ForFace<1>(inv_jac_ptr));
// push second base directives tp physical element shape
new OpSetInvJacH1ForFace<2>(inv_jac_ptr));
// calculate value of function at integration points
// calculate gradient at integration points
FIELD_NAME, common_data_ptr->approxGradVals));
// calculate hessian at integration points
FIELD_NAME, common_data_ptr->approxHessianVals));
//FIXME: Note third and forth derivative is calculated but not properly tested
// push operator to evaluate errors
new OpError(common_data_ptr));
// zero error vector and iterate over all elements on the mesh
CHKERR VecZeroEntries(common_data_ptr->L2Vec);
CHKERR pipeline_mng->loopFiniteElements();
// assemble error vector
CHKERR VecAssemblyBegin(common_data_ptr->L2Vec);
CHKERR VecAssemblyEnd(common_data_ptr->L2Vec);
// check if errors are small and print results
constexpr double eps = 1e-8;
const double *array;
CHKERR VecGetArrayRead(common_data_ptr->L2Vec, &array);
MOFEM_LOG_C("WORLD", Sev::inform,
"Error %6.4e Diff Error %6.4e Diff2 Error %6.4e\n",
std::sqrt(array[0]), std::sqrt(array[1]), std::sqrt(array[2]));
if (std::sqrt(array[0]) > eps)
if (std::sqrt(array[1]) > eps)
"Wrong function first direcative");
if (std::sqrt(array[2]) > eps)
"Wrong function second direcative");
CHKERR VecRestoreArrayRead(common_data_ptr->L2Vec, &array);
//! [Check results]
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
MoFEM::Core::Initialize(&argc, &argv, NULL, help);
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database insterface
//! [Create MoFEM]
//! [AtomTest]
AtomTest ex(m_field);
CHKERR ex.runProblem();
//! [AtomTest]
static Index< 'p', 3 > p
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:304
MoFEM::FaceElementForcesAndSourcesCore FaceEle
static const double eps
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
Catch errors.
Definition: definitions.h:372
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
Definition: definitions.h:40
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:470
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
#define MOFEM_LOG(channel, severity)
Definition: LogManager.hpp:301
FTensor::Index< 'i', SPACE_DIM > i
auto diff2_fun
Function second derivative.
int main(int argc, char *argv[])
[Check results]
EntitiesFieldData::EntData EntData
Data on entities.
static char help[]
constexpr char FIELD_NAME[]
auto fun
Function to approximate.
constexpr int SPACE_DIM
constexpr int FIELD_DIM
constexpr int BASE_DIM
auto diff_fun
Function derivative.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, FIELD_DIM > OpDomainMass
OPerator to integrate mass matrix for least square approximation.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, FIELD_DIM > OpDomainSource
Operator to integrate the right hand side matrix for the problem.
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
auto createSmartVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1086
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
const double D
boost::shared_ptr< VectorDouble > approxVals
boost::shared_ptr< MatrixDouble > approxGradVals
SmartPetscObj< Vec > L2Vec
boost::shared_ptr< MatrixDouble > invJacPtr
boost::shared_ptr< MatrixDouble > approxHessianVals
Operator to evaluate errors.
OpError(boost::shared_ptr< CommonData > &common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Simple * simpleInterface
MoFEMErrorCode checkResults()
[Check results]
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
AtomTest(MoFEM::Interface &m_field)
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode readMesh()
[Run programme]
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
MoFEM::Interface & mField
MoFEMErrorCode runProblem()
[Run programme]
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Core (interface) class.
Definition: Core.hpp:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
operator doWork function is executed on FE rows
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Evaluate field gradient values for scalar field, i.e. gradient is tensor rank 1 (vector),...
Get value at integration points for scalar field.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.