No Matches
SCL-0: Least square approximaton
Prerequisites of this tutorial include MSH-1: Create a 2D mesh from Gmsh

Intended learning outcome:
  • general structure of a program developed using MoFEM
  • idea of Simple Interface in MoFEM and how to use it
  • idea of Domain element in MoFEM and how to use it
  • Use default Forms Integrals
  • how to push the developed UDOs to the Pipeline
  • utilisation of tools to convert outputs (MOAB) and visualise them (Paraview)
* \file approximation.cpp
* \example approximation.cpp
* Using PipelineManager interface calculate the divergence of base functions,
* and integral of flux on the boundary. Since the h-div space is used, volume
* integral and boundary integral should give the same result.
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
constexpr char FIELD_NAME[] = "U";
constexpr int FIELD_DIM = 1;
constexpr int SPACE_DIM = 2;
template <int DIM> struct ElementsAndOps {};
template <> struct ElementsAndOps<2> {
template <> struct ElementsAndOps<3> {
template <int FIELD_DIM> struct ApproxFieldFunction;
template <> struct ApproxFieldFunction<1> {
double operator()(const double x, const double y, const double z) {
return sin(x * 10.) * cos(y * 10.);
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
struct CommonData {
boost::shared_ptr<VectorDouble> approxVals;
boost::shared_ptr<CommonData> commonDataPtr;
template <int FIELD_DIM> struct OpError;
template <> struct Example::OpError<1> : public DomainEleOp {
boost::shared_ptr<CommonData> commonDataPtr;
OpError(boost::shared_ptr<CommonData> &common_data_ptr)
: DomainEleOp(FIELD_NAME, OPROW), commonDataPtr(common_data_ptr) {}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
if (const size_t nb_dofs = data.getIndices().size()) {
const int nb_integration_pts = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
auto t_coords = getFTensor1CoordsAtGaussPts();
VectorDouble nf(nb_dofs, false);
const double volume = getMeasure();
auto t_row_base = data.getFTensor0N();
double error = 0;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = t_w * volume;
double diff = t_val - Example::approxFunction(t_coords(0), t_coords(1),
error += alpha * pow(diff, 2);
for (size_t r = 0; r != nb_dofs; ++r) {
nf[r] += alpha * t_row_base * diff;
const int index = 0;
CHKERR VecSetValue(commonDataPtr->L2Vec, index, error, ADD_VALUES);
CHKERR VecSetValues(commonDataPtr->resVec, data, &nf[0], ADD_VALUES);
//! [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]
//! [Set integration rule]
auto rule = [](int, int, int p) -> int { return 2 * p; };
CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
//! [Set integration rule]
//! [Create common data]
commonDataPtr = boost::make_shared<CommonData>();
mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
//! [Create common data]
//! [Boundary condition]
//! [Boundary condition]
//! [Push operators to pipeline]
auto beta = [](const double, const double, const double) { return 1; };
//! [Push operators to pipeline]
//! [Solve]
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);
//! [Solve]
auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
auto u_ptr = boost::make_shared<VectorDouble>();
new OpPPMap(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{FIELD_NAME, u_ptr}},
pipeline_mng->getDomainRhsFE() = post_proc_fe;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR post_proc_fe->writeFile("out_approx.h5m");
//! [Postprocess results]
//! [Check results]
CHKERR pipeline_mng->loopFiniteElements();
CHKERR VecAssemblyBegin(commonDataPtr->L2Vec);
CHKERR VecAssemblyEnd(commonDataPtr->L2Vec);
CHKERR VecAssemblyBegin(commonDataPtr->resVec);
CHKERR VecAssemblyEnd(commonDataPtr->resVec);
double nrm2;
CHKERR VecNorm(commonDataPtr->resVec, NORM_2, &nrm2);
const double *array;
CHKERR VecGetArrayRead(commonDataPtr->L2Vec, &array);
if (mField.get_comm_rank() == 0)
PetscPrintf(PETSC_COMM_SELF, "Error %6.4e Vec norm %6.4e\n", std::sqrt(array[0]),
CHKERR VecRestoreArrayRead(commonDataPtr->L2Vec, &array);
constexpr double eps = 1e-8;
if (nrm2 > eps)
"Not converged solution");
//! [Check results]
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
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]
//! [Example]
Example ex(m_field);
CHKERR ex.runProblem();
//! [Example]
static Index< 'p', 3 > p
std::string param_file
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
Catch errors.
Definition: definitions.h:372
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ 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: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
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
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
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.
double D
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
int r
Definition: sdf.py:8
FTensor::Index< 'i', 3 > i
SmartPetscObj< Vec > L2Vec
SmartPetscObj< Vec > resVec
boost::shared_ptr< VectorDouble > approxVals
Definition: plastic.cpp:226
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
MoFEMErrorCode readMesh()
[Run problem]
MoFEMErrorCode setIntegrationRules()
[Set up problem]
Definition: plot_base.cpp:203
static ApproxFieldFunction< FIELD_DIM > approxFunction
boost::shared_ptr< CommonData > commonDataPtr
Definition: integration.cpp:42
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode solveSystem()
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:440
Simple * simpleInterface
Definition: helmholtz.cpp:41
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:262
MoFEM::Interface & mField
Definition: plastic.cpp:233
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:274
MoFEMErrorCode outputResults()
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
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
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:27
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.