v0.14.0
dg_projection.cpp

Testing DG projection operators

/**
* \example dg_projection.cpp
*
* Testing DG projection operators
*
*/
#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;
constexpr int order = 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 + y + x * x + y * y;
};
/**
* @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 readMesh();
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> 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<MatrixDouble> data_ptr)
: DomainEleOp(NOSPACE, OPSPACE), dataPtr(data_ptr) {}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
const int nb_integration_pts = getGaussPts().size2();
auto t_val = getFTensor1FromMat<1>(
*(dataPtr)); // get function approximation at gauss pts
auto t_coords = getFTensor1CoordsAtGaussPts(); // get coordinates of
// integration points
for (int gg = 0; gg != nb_integration_pts; ++gg) {
// Calculate errors
double diff = t_val(0) - fun(t_coords(0), t_coords(1), t_coords(2));
constexpr double eps = 1e-8;
if (std::abs(diff) > eps) {
MOFEM_LOG("SELF", Sev::error) << "Wrong function value " << diff;
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong function value");
}
// move data to next integration point
++t_val;
++t_coords;
}
MOFEM_LOG("SELF", Sev::noisy) << "All is OK";
}
private:
boost::shared_ptr<MatrixDouble> dataPtr;
};
//! [Run programme]
CHKERR readMesh();
CHKERR setupProblem();
CHKERR assembleSystem();
CHKERR solveSystem();
CHKERR checkResults();
}
//! [Run programme]
//! [Read mesh]
CHKERR mField.getInterface(simpleInterface);
CHKERR simpleInterface->getOptions();
CHKERR simpleInterface->loadFile();
}
//! [Read mesh]
//! [Set up problem]
// Add field
CHKERR simpleInterface->addDomainField(FIELD_NAME, H1,
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]
auto simple = mField.getInterface<Simple>();
auto 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 + 1; };
CHKERR pipeline_mng->setDomainRhsIntegrationRule(
rule); // set integration rule
auto entity_data_l2 = boost::make_shared<EntitiesFieldData>(
MBENTITYSET); // entity data shared between
// physical and post proc
// elements
auto mass_ptr = boost::make_shared<MatrixDouble>(); // integrated mass matrix
// of post proc element
auto coeffs_ptr =
boost::make_shared<MatrixDouble>(); // vector of coeffs shared between
// physical and post proc elements
auto data_ptr =
boost::make_shared<MatrixDouble>(); // data stored at integration points
// of the physical element and
// evaluated at integration points of
// the post proc element
auto op_this =
new OpLoopThis<DomainEle>(mField, simple->getDomainFEName(), Sev::noisy);
pipeline_mng->getOpDomainRhsPipeline().push_back(op_this); // 1
pipeline_mng->getOpDomainRhsPipeline().push_back(new OpDGProjectionEvaluation(
data_ptr, coeffs_ptr, entity_data_l2, AINSWORTH_LEGENDRE_BASE,
L2)); // 5
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpError(data_ptr)); // 6
auto fe_physics_ptr = op_this->getThisFEPtr();
fe_physics_ptr->getRuleHook = [](int, int, int p) { return 2 * p; };
fe_physics_ptr->getOpPtrVector().push_back(new OpDGProjectionMassMatrix(
order, mass_ptr, entity_data_l2, AINSWORTH_LEGENDRE_BASE, L2)); // 2
fe_physics_ptr->getOpPtrVector().push_back(
data_ptr)); // 3
fe_physics_ptr->getOpPtrVector().push_back(
new OpDGProjectionCoefficients(data_ptr, coeffs_ptr, mass_ptr,
entity_data_l2, AINSWORTH_LEGENDRE_BASE,
L2)); // 4
CHKERR pipeline_mng->loopFiniteElements();
}
//! [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
NOSPACE
@ NOSPACE
Definition: definitions.h:83
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:128
EntData
EntitiesFieldData::EntData EntData
Data on entities.
Definition: dg_projection.cpp:30
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
OpError
Definition: initial_diffusion.cpp:61
MoFEM::OpLoopThis
Execute "this" element in the operator.
Definition: DGProjection.hpp:68
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
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
main
int main(int argc, char *argv[])
[Check results]
Definition: dg_projection.cpp:268
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
MoFEM::OpDGProjectionCoefficients
Definition: DGProjection.hpp:119
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
BASE_DIM
constexpr int BASE_DIM
Definition: dg_projection.cpp:15
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
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
AtomTest::OpError
Operator to evaluate errors.
Definition: child_and_parent.cpp:82
SPACE_DIM
constexpr int SPACE_DIM
Definition: dg_projection.cpp:17
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, FIELD_DIM > OpDomainMass
OPerator to integrate mass matrix for least square approximation.
Definition: dg_projection.cpp:45
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: dg_projection.cpp:52
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
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:317
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
FIELD_DIM
constexpr int FIELD_DIM
Definition: dg_projection.cpp:16
BiLinearForm
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
ElementsAndOps
Definition: child_and_parent.cpp:18
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
FIELD_NAME
constexpr char FIELD_NAME[]
Definition: dg_projection.cpp:14
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MoFEM::OpDGProjectionEvaluation
Definition: DGProjection.hpp:136
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
CommonData
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:22
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
AtomTest::checkResults
MoFEMErrorCode checkResults()
[Check results]
Definition: dg_projection.cpp:213
AtomTest::readMesh
MoFEMErrorCode readMesh()
[Run programme]
Definition: child_and_parent.cpp:208
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
help
static char help[]
Definition: dg_projection.cpp:12
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
fun
auto fun
Function to approximate.
Definition: dg_projection.cpp:36
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:429
AtomTest::setupProblem
MoFEMErrorCode setupProblem()
[Read mesh]
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:359
F
@ F
Definition: free_surface.cpp:394
MoFEM::OpDGProjectionMassMatrix
Definition: DGProjection.hpp:106