v0.14.0
dg_projection.cpp
Go to the documentation of this file.
1 /**
2  * \example dg_projection.cpp
3  *
4  * Testing DG projection operators
5  *
6  */
7 
8 #include <MoFEM.hpp>
9 
10 using namespace MoFEM;
11 
12 static char help[] = "...\n\n";
13 
14 constexpr char FIELD_NAME[] = "U";
15 constexpr int BASE_DIM = 1;
16 constexpr int FIELD_DIM = 1;
17 constexpr int SPACE_DIM = 2;
18 constexpr int order = 2;
19 
20 template <int DIM> struct ElementsAndOps {};
21 
22 template <> struct ElementsAndOps<2> {
25 };
26 
27 using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle; ///< Finite elenent type
28 using DomainEleOp =
29  DomainEle::UserDataOperator; ///< Finire element operator type
30 using EntData = EntitiesFieldData::EntData; ///< Data on entities
31 
32 /**
33  * @brief Function to approximate
34  *
35  */
36 auto fun = [](const double x, const double y, const double z) {
37  return x + y + x * x + y * y;
38 };
39 
40 /**
41  * @brief OPerator to integrate mass matrix for least square approximation
42  *
43  */
46 
47 /**
48  * @brief Operator to integrate the right hand side matrix for the problem
49  *
50  */
52  PETSC>::LinearForm<GAUSS>::OpSource<BASE_DIM, FIELD_DIM>;
53 
54 struct AtomTest {
55 
56  AtomTest(MoFEM::Interface &m_field) : mField(m_field) {}
57 
58  MoFEMErrorCode runProblem();
59 
60 private:
61  MoFEM::Interface &mField;
62  Simple *simpleInterface;
63 
64  MoFEMErrorCode readMesh();
65  MoFEMErrorCode setupProblem();
66  MoFEMErrorCode assembleSystem();
67  MoFEMErrorCode solveSystem();
68  MoFEMErrorCode checkResults();
69 
70  /**
71  * @brief Collected data use d by operator to evaluate errors for the test
72  *
73  */
74  struct CommonData {
75  boost::shared_ptr<MatrixDouble> invJacPtr;
76  boost::shared_ptr<VectorDouble> approxVals;
77  boost::shared_ptr<MatrixDouble> approxGradVals;
78  boost::shared_ptr<MatrixDouble> approxHessianVals;
79  SmartPetscObj<Vec> L2Vec;
80  };
81 
82  /**
83  * @brief Operator to evaluate errors
84  *
85  */
86  struct OpError;
87 };
88 
89 /**
90  * @brief Operator to evaluate errors
91  *
92  */
93 struct AtomTest::OpError : public DomainEleOp {
94  boost::shared_ptr<CommonData> commonDataPtr;
95  OpError(boost::shared_ptr<MatrixDouble> data_ptr)
96  : DomainEleOp(NOSPACE, OPSPACE), dataPtr(data_ptr) {}
97 
98  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
100 
101  const int nb_integration_pts = getGaussPts().size2();
102  auto t_val = getFTensor1FromMat<1>(
103  *(dataPtr)); // get function approximation at gauss pts
104  auto t_coords = getFTensor1CoordsAtGaussPts(); // get coordinates of
105  // integration points
106 
107  for (int gg = 0; gg != nb_integration_pts; ++gg) {
108 
109  // Calculate errors
110 
111  double diff = t_val(0) - fun(t_coords(0), t_coords(1), t_coords(2));
112  constexpr double eps = 1e-8;
113  if (std::abs(diff) > eps) {
114  MOFEM_LOG("SELF", Sev::error) << "Wrong function value " << diff;
115  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
116  "Wrong function value");
117  }
118 
119  // move data to next integration point
120  ++t_val;
121  ++t_coords;
122  }
123 
124  MOFEM_LOG("SELF", Sev::noisy) << "All is OK";
125 
127  }
128 
129  private:
130  boost::shared_ptr<MatrixDouble> dataPtr;
131 };
132 
133 //! [Run programme]
136  CHKERR readMesh();
137  CHKERR setupProblem();
138  CHKERR assembleSystem();
139  CHKERR solveSystem();
140  CHKERR checkResults();
142 }
143 //! [Run programme]
144 
145 //! [Read mesh]
148 
149  CHKERR mField.getInterface(simpleInterface);
150  CHKERR simpleInterface->getOptions();
151  CHKERR simpleInterface->loadFile();
152 
154 }
155 //! [Read mesh]
156 
157 //! [Set up problem]
160  // Add field
161  CHKERR simpleInterface->addDomainField(FIELD_NAME, H1,
163  CHKERR simpleInterface->setFieldOrder(FIELD_NAME, order);
164  CHKERR simpleInterface->setUp();
165 
167 }
168 //! [Set up problem]
169 
170 //! [Push operators to pipeline]
173 
174  auto rule = [](int, int, int p) -> int { return 2 * p; };
175 
176  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
177  CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
178  CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
179 
180  auto beta = [](const double, const double, const double) { return 1; };
181  pipeline_mng->getOpDomainLhsPipeline().push_back(
182  new OpDomainMass(FIELD_NAME, FIELD_NAME, beta));
183  pipeline_mng->getOpDomainRhsPipeline().push_back(
185 
187 }
188 //! [Push operators to pipeline]
189 
190 //! [Solve]
193  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
194 
195  MOFEM_LOG("WORLD", Sev::inform) << "Solve problem";
196 
197  auto solver = pipeline_mng->createKSP();
198  CHKERR KSPSetFromOptions(solver);
199  CHKERR KSPSetUp(solver);
200 
201  auto dm = simpleInterface->getDM();
202  auto D = createDMVector(dm);
203  auto F = vectorDuplicate(D);
204 
205  CHKERR KSPSolve(solver, F, D);
206  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
207  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
208  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
210 }
211 
212 //! [Check results]
215  auto simple = mField.getInterface<Simple>();
216  auto pipeline_mng = mField.getInterface<PipelineManager>();
217  pipeline_mng->getDomainLhsFE().reset();
218  pipeline_mng->getDomainRhsFE().reset();
219  pipeline_mng->getOpDomainRhsPipeline().clear();
220 
221  auto rule = [](int, int, int p) -> int { return 2 * p + 1; };
222  CHKERR pipeline_mng->setDomainRhsIntegrationRule(
223  rule); // set integration rule
224 
225  auto entity_data_l2 = boost::make_shared<EntitiesFieldData>(
226  MBENTITYSET); // entity data shared between
227  // physical and post proc
228  // elements
229  auto mass_ptr = boost::make_shared<MatrixDouble>(); // integrated mass matrix
230  // of post proc element
231  auto coeffs_ptr =
232  boost::make_shared<MatrixDouble>(); // vector of coeffs shared between
233  // physical and post proc elements
234  auto data_ptr =
235  boost::make_shared<MatrixDouble>(); // data stored at integration points
236  // of the physical element and
237  // evaluated at integration points of
238  // the post proc element
239 
240  auto op_this =
241  new OpLoopThis<DomainEle>(mField, simple->getDomainFEName(), Sev::noisy);
242  pipeline_mng->getOpDomainRhsPipeline().push_back(op_this); // 1
243  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpDGProjectionEvaluation(
244  data_ptr, coeffs_ptr, entity_data_l2, AINSWORTH_LEGENDRE_BASE,
245  L2)); // 5
246  pipeline_mng->getOpDomainRhsPipeline().push_back(
247  new OpError(data_ptr)); // 6
248 
249  auto fe_physics_ptr = op_this->getThisFEPtr();
250  fe_physics_ptr->getRuleHook = [](int, int, int p) { return 2 * p; };
251 
252  fe_physics_ptr->getOpPtrVector().push_back(new OpDGProjectionMassMatrix(
253  order, mass_ptr, entity_data_l2, AINSWORTH_LEGENDRE_BASE, L2)); // 2
254  fe_physics_ptr->getOpPtrVector().push_back(
256  data_ptr)); // 3
257  fe_physics_ptr->getOpPtrVector().push_back(
258  new OpDGProjectionCoefficients(data_ptr, coeffs_ptr, mass_ptr,
259  entity_data_l2, AINSWORTH_LEGENDRE_BASE,
260  L2)); // 4
261 
262  CHKERR pipeline_mng->loopFiniteElements();
263 
265 }
266 //! [Check results]
267 
268 int main(int argc, char *argv[]) {
269 
270  // Initialisation of MoFEM/PETSc and MOAB data structures
271  MoFEM::Core::Initialize(&argc, &argv, NULL, help);
272 
273  try {
274 
275  //! [Register MoFEM discrete manager in PETSc]
276  DMType dm_name = "DMMOFEM";
277  CHKERR DMRegister_MoFEM(dm_name);
278  //! [Register MoFEM discrete manager in PETSc
279 
280  //! [Create MoAB]
281  moab::Core mb_instance; ///< mesh database
282  moab::Interface &moab = mb_instance; ///< mesh database interface
283  //! [Create MoAB]
284 
285  //! [Create MoFEM]
286  MoFEM::Core core(moab); ///< finite element database
287  MoFEM::Interface &m_field = core; ///< finite element database insterface
288  //! [Create MoFEM]
289 
290  //! [AtomTest]
291  AtomTest ex(m_field);
292  CHKERR ex.runProblem();
293  //! [AtomTest]
294  }
295  CATCH_ERRORS;
296 
298 }
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:127
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
MoFEM::OpLoopThis
Execute "this" element in the operator.
Definition: DGProjection.hpp:68
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
AtomTest::CommonData::approxHessianVals
boost::shared_ptr< MatrixDouble > approxHessianVals
Definition: dg_projection.cpp:78
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:1975
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:535
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1067
AtomTest::CommonData::invJacPtr
boost::shared_ptr< MatrixDouble > invJacPtr
Definition: dg_projection.cpp:75
AtomTest::CommonData::approxGradVals
boost::shared_ptr< MatrixDouble > approxGradVals
Definition: dg_projection.cpp:77
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
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
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:310
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::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MoFEM::PipelineManager::setDomainRhsIntegrationRule
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:530
AtomTest::AtomTest
AtomTest(MoFEM::Interface &m_field)
Definition: dg_projection.cpp:56
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
AtomTest::OpError::OpError
OpError(boost::shared_ptr< MatrixDouble > data_ptr)
Definition: dg_projection.cpp:95
MoFEM::PipelineManager::setDomainLhsIntegrationRule
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:503
AtomTest::OpError::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: dg_projection.cpp:98
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:372
AtomTest::OpError::dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
Definition: dg_projection.cpp:130
MoFEM::OpDGProjectionEvaluation
Definition: DGProjection.hpp:136
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1102
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
AtomTest::OpError::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: dg_projection.cpp:94
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:416
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:346
F
@ F
Definition: free_surface.cpp:394
MoFEM::OpDGProjectionMassMatrix
Definition: DGProjection.hpp:106