v0.14.0
Loading...
Searching...
No Matches
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
*
*/
struct AtomTest {
AtomTest(MoFEM::Interface &m_field) : mField(m_field) {}
private:
/**
* @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]
}
//! [Run programme]
//! [Read mesh]
CHKERR simpleInterface->getOptions();
CHKERR simpleInterface->loadFile();
}
//! [Read mesh]
//! [Set up problem]
// Add field
CHKERR simpleInterface->addDomainField(FIELD_NAME, H1,
}
//! [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; };
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
}
//! [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 = 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 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]
}
}
static Index< 'p', 3 > p
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
static char help[]
int main()
Definition: adol-c_atom.cpp:46
static const double eps
constexpr char FIELD_NAME[]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
constexpr int SPACE_DIM
constexpr int FIELD_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
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
@ NOSPACE
Definition: definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
#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
auto fun
Function to approximate.
constexpr int BASE_DIM
constexpr int order
@ F
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:509
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
constexpr char FIELD_NAME[]
auto fun
Function to approximate.
constexpr int FIELD_DIM
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.
double D
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
boost::shared_ptr< VectorDouble > approxVals
SmartPetscObj< Vec > L2Vec
boost::shared_ptr< MatrixDouble > invJacPtr
boost::shared_ptr< MatrixDouble > approxHessianVals
boost::shared_ptr< MatrixDouble > approxGradVals
Operator to evaluate errors.
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< MatrixDouble > dataPtr
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]
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode readMesh()
[Run programme]
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
MoFEM::Interface & mField
MoFEMErrorCode runProblem()
[Run programme]
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
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
@ OPSPACE
operator do Work is execute on space data
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Execute "this" element in the operator.
PipelineManager interface.
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.