v0.14.0
higher_derivatives.cpp

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
*
*/
PETSC>::LinearForm<GAUSS>::OpSource<BASE_DIM, FIELD_DIM>;
struct AtomTest {
AtomTest(MoFEM::Interface &m_field) : mField(m_field) {}
MoFEMErrorCode runProblem();
private:
Simple *simpleInterface;
MoFEMErrorCode setupProblem();
MoFEMErrorCode assembleSystem();
MoFEMErrorCode solveSystem();
MoFEMErrorCode checkResults();
/**
* @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> 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;
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
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
*(commonDataPtr
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
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));
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;
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"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
++t_w;
++t_val;
++t_hessian_val;
++t_inv_jac;
++t_coords;
}
// assemble error vector
std::array<int, 3> index = {0, 1, 2};
CHKERR VecSetValues(commonDataPtr->L2Vec, 3, index.data(), error.data(),
}
};
//! [Run programme]
CHKERR setupProblem();
CHKERR assembleSystem();
CHKERR solveSystem();
CHKERR checkResults();
}
//! [Run programme]
CHKERR mField.getInterface(simpleInterface);
CHKERR simpleInterface->getOptions();
}
//! [Set up problem]
constexpr int order = 4;
CHKERR simpleInterface->setFieldOrder(FIELD_NAME, order);
CHKERR simpleInterface->setUp();
}
//! [Set up problem]
//! [Push operators to pipeline]
auto rule = [](int, int, int p) -> int { return 2 * p; };
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
auto beta = [](const double, const double, const double) { return 1; };
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
}
//! [Push operators to pipeline]
//! [Solve]
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
MOFEM_LOG("WORLD", Sev::inform) << "Solve problem";
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSetUp(solver);
auto dm = simpleInterface->getDM();
auto D = createDMVector(dm);
auto F = vectorDuplicate(D);
CHKERR KSPSolve(solver, F, D);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
}
//! [Check results]
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getDomainRhsFE().reset();
pipeline_mng->getOpDomainRhsPipeline().clear();
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 = createVectorMPI(
mField.get_comm(), (!mField.get_comm_rank()) ? 3 : 0, 3);
common_data_ptr->approxVals = boost::make_shared<VectorDouble>();
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
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateHOJacForFace(jac_ptr));
// calculate incerse of jacobian at integration points
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
// calculate mass matrix to project derivatives
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpBaseDerivativesMass<BASE_DIM>(base_mass, data_l2,
// calculate second derivative of base functions, i.e. hessian
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpBaseDerivativesNext<BASE_DIM>(BaseDerivatives::SecondDerivative,
base_mass, data_l2,
// calculate third derivative
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpBaseDerivativesNext<BASE_DIM>(BaseDerivatives::ThirdDerivative,
base_mass, data_l2,
// calculate forth derivative
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpBaseDerivativesNext<BASE_DIM>(BaseDerivatives::ForthDerivative,
base_mass, data_l2,
// push first base directives tp physical element shape
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSetInvJacH1ForFace<1>(inv_jac_ptr));
// push second base directives tp physical element shape
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSetInvJacH1ForFace<2>(inv_jac_ptr));
// calculate value of function at integration points
pipeline_mng->getOpDomainRhsPipeline().push_back(
common_data_ptr->approxVals));
// calculate gradient at integration points
pipeline_mng->getOpDomainRhsPipeline().push_back(
// calculate hessian at integration points
pipeline_mng->getOpDomainRhsPipeline().push_back(
FIELD_NAME, common_data_ptr->approxHessianVals));
//FIXME: Note third and forth derivative is calculated but not properly tested
// push operator to evaluate errors
pipeline_mng->getOpDomainRhsPipeline().push_back(
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;
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)
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Wrong function value");
if (std::sqrt(array[1]) > eps)
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong function first direcative");
if (std::sqrt(array[2]) > eps)
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong function second direcative");
}
//! [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]
}
}
AtomTest
Definition: child_and_parent.cpp:57
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::PipelineManager::getDomainRhsFE
boost::shared_ptr< FEMethod > & getDomainRhsFE()
Definition: PipelineManager.hpp:405
FIELD_DIM
constexpr int FIELD_DIM
Definition: higher_derivatives.cpp:18
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
OpError
Definition: initial_diffusion.cpp:61
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
MoFEM::PipelineManager::loopFiniteElements
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
Definition: PipelineManager.cpp:19
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::PipelineManager::getOpDomainRhsPipeline
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
Definition: PipelineManager.hpp:773
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
AtomTest::assembleSystem
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
Definition: child_and_parent.cpp:290
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
Definition: child_and_parent.cpp:53
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1576
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1294
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MoFEM::OpCalculateScalarFieldHessian
Evaluate field gradient values for scalar field, i.e. gradient is tensor rank 1 (vector),...
Definition: UserDataOperators.hpp:1307
AtomTest::OpError
Operator to evaluate errors.
Definition: child_and_parent.cpp:82
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1067
diff2_fun
auto diff2_fun
Function second derivative.
Definition: higher_derivatives.cpp:55
BASE_DIM
constexpr int BASE_DIM
Definition: higher_derivatives.cpp:17
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
double
convert.type
type
Definition: convert.py:64
MoFEM::PipelineManager::FaceEle
MoFEM::FaceElementForcesAndSourcesCore FaceEle
Definition: PipelineManager.hpp:35
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:310
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::PipelineManager::createKSP
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
Definition: PipelineManager.cpp:49
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::PipelineManager::setDomainRhsIntegrationRule
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:530
MoFEM::PipelineManager::getOpDomainLhsPipeline
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
Definition: PipelineManager.hpp:749
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
BiLinearForm
FIELD_NAME
constexpr char FIELD_NAME[]
Definition: higher_derivatives.cpp:16
FTensor::Index< 'i', 2 >
AtomTest::solveSystem
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
Definition: child_and_parent.cpp:377
MoFEM::PipelineManager::setDomainLhsIntegrationRule
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:503
diff_fun
auto diff_fun
Function derivative.
Definition: higher_derivatives.cpp:45
ElementsAndOps
Definition: child_and_parent.cpp:18
MoFEM::PipelineManager::getDomainLhsFE
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Definition: PipelineManager.hpp:401
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
fun
auto fun
Function to approximate.
Definition: higher_derivatives.cpp:37
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
help
static char help[]
Definition: higher_derivatives.cpp:14
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1102
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
eps
static const double eps
OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, FIELD_DIM > OpDomainMass
OPerator to integrate mass matrix for least square approximation.
Definition: higher_derivatives.cpp:67
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
EntData
EntitiesFieldData::EntData EntData
Data on entities.
Definition: higher_derivatives.cpp:31
AtomTest::checkResults
MoFEMErrorCode checkResults()
[Check results]
Definition: dg_projection.cpp:213
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3254
MoFEM::OpCalculateHOJacForFace
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
Definition: HODataOperators.hpp:264
SPACE_DIM
constexpr int SPACE_DIM
Definition: higher_derivatives.cpp:19
[Run programme]
Definition: child_and_parent.cpp:208
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
main
int main(int argc, char *argv[])
[Check results]
Definition: higher_derivatives.cpp:397
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::SmartPetscObj< Vec >
OpDomainSource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
Definition: child_and_parent.cpp:55
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
OpDomainSource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, FIELD_DIM > OpDomainSource
Operator to integrate the right hand side matrix for the problem.
Definition: higher_derivatives.cpp:74
convert.int
int
Definition: convert.py:64
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
AtomTest::setupProblem
MoFEMErrorCode setupProblem()
Definition: child_and_parent.cpp:268
AtomTest::runProblem
MoFEMErrorCode runProblem()
[Run programme]
Definition: child_and_parent.cpp:188
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
F
@ F
Definition: free_surface.cpp:394