v0.14.0
child_and_parent.cpp
Go to the documentation of this file.
1 /**
2  * \example child_and_parent.cpp
3  *
4  * Testing projection child and parent.
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 FIELD_DIM = 1;
16 constexpr int SPACE_DIM = 2;
17 
18 template <int DIM> struct ElementsAndOps {};
19 
20 template <> struct ElementsAndOps<2> {
21  using DomainEle = PipelineManager::FaceEle;
22  using DomainEleOp = DomainEle::UserDataOperator;
23  using DomianParentEle = FaceElementForcesAndSourcesCoreOnChildParent;
24  using BoundaryEle = PipelineManager::EdgeEle;
25  using BoundaryEleOp = BoundaryEle::UserDataOperator;
26  using BoundaryParentEle = EdgeElementForcesAndSourcesCoreOnChildParent;
27 };
28 
29 template <> struct ElementsAndOps<3> {
30  using DomainEle = VolumeElementForcesAndSourcesCore;
31  using DomainEleOp = DomainEle::UserDataOperator;
32 };
33 
34 using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle;
35 using DomainParentEle = ElementsAndOps<SPACE_DIM>::DomianParentEle;
36 using DomainEleOp = DomainEle::UserDataOperator;
37 using EntData = EntitiesFieldData::EntData;
38 
39 using BoundaryEle = ElementsAndOps<SPACE_DIM>::BoundaryEle;
40 using BoundaryEleOp = BoundaryEle::UserDataOperator;
41 using BoundaryParentEle = ElementsAndOps<SPACE_DIM>::BoundaryParentEle;
42 
43 template <int FIELD_DIM> struct ApproxFieldFunction;
44 
45 template <> struct ApproxFieldFunction<1> {
46  double operator()(const double x, const double y, const double z) {
47  return x * x + y * y + x * y + pow(x, 3) + pow(y, 3) + pow(x, 4) +
48  pow(y, 4);
49  }
50 };
51 
52 using OpDomainMass = FormsIntegrators<DomainEleOp>::Assembly<
53  PETSC>::BiLinearForm<GAUSS>::OpMass<1, FIELD_DIM>;
54 using OpDomainSource = FormsIntegrators<DomainEleOp>::Assembly<
55  PETSC>::LinearForm<GAUSS>::OpSource<1, FIELD_DIM>;
56 
57 struct AtomTest {
58 
59  AtomTest(MoFEM::Interface &m_field) : mField(m_field) {}
60 
61  MoFEMErrorCode runProblem();
62 
63 private:
64  MoFEM::Interface &mField;
65  Simple *simpleInterface;
66 
67  static ApproxFieldFunction<FIELD_DIM> approxFunction;
68 
69  MoFEMErrorCode readMesh();
70  MoFEMErrorCode setupProblem();
71  MoFEMErrorCode assembleSystem();
72  MoFEMErrorCode solveSystem();
73  MoFEMErrorCode
74  checkResults(boost::function<bool(FEMethod *fe_method_ptr)> test_bit);
75  MoFEMErrorCode refineResults();
76  struct CommonData {
77  boost::shared_ptr<VectorDouble> approxVals;
78  SmartPetscObj<Vec> L2Vec;
79  SmartPetscObj<Vec> resVec;
80  };
81 
82  template <int FIELD_DIM> struct OpError;
83 };
84 
85 ApproxFieldFunction<FIELD_DIM> AtomTest::approxFunction =
86  ApproxFieldFunction<FIELD_DIM>();
87 
88 template <> struct AtomTest::OpError<1> : public DomainEleOp {
89  boost::shared_ptr<CommonData> commonDataPtr;
90  OpError(boost::shared_ptr<CommonData> &common_data_ptr)
91  : DomainEleOp(FIELD_NAME, OPROW), commonDataPtr(common_data_ptr) {}
92  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
93  MoFEMFunctionBegin;
94 
95  if (const size_t nb_dofs = data.getIndices().size()) {
96 
97  const int nb_integration_pts = getGaussPts().size2();
98  auto t_w = getFTensor0IntegrationWeight();
99  auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
100  auto t_coords = getFTensor1CoordsAtGaussPts();
101 
102  VectorDouble nf(nb_dofs, false);
103  nf.clear();
104 
105  FTensor::Index<'i', 3> i;
106  const double volume = getMeasure();
107 
108  auto t_row_base = data.getFTensor0N();
109  double error = 0;
110  for (int gg = 0; gg != nb_integration_pts; ++gg) {
111 
112  const double alpha = t_w * volume;
113  double diff = t_val - AtomTest::approxFunction(t_coords(0), t_coords(1),
114  t_coords(2));
115  error += alpha * pow(diff, 2);
116 
117  for (size_t r = 0; r != nb_dofs; ++r) {
118  nf[r] += alpha * t_row_base * diff;
119  ++t_row_base;
120  }
121 
122  ++t_w;
123  ++t_val;
124  ++t_coords;
125  }
126 
127  const int index = 0;
128  CHKERR VecSetValue(commonDataPtr->L2Vec, index, error, ADD_VALUES);
129  CHKERR VecSetValues(commonDataPtr->resVec, data, &nf[0], ADD_VALUES);
130  }
131 
132  MoFEMFunctionReturn(0);
133  }
134 };
135 
136 template <typename ELE_OP, typename PARENT_ELE>
137 struct OpCheckGaussCoords : public ELE_OP {
138  OpCheckGaussCoords() : ELE_OP(NOSPACE, ELE_OP::OPSPACE) {}
139 
140  MoFEMErrorCode doWork(int side, EntityType type,
141  EntitiesFieldData::EntData &data) {
142  MoFEMFunctionBegin;
143 
144  MatrixDouble parent_coords;
145 
146  PARENT_ELE parent_fe(this->getPtrFE()->mField);
147  auto op = new ELE_OP(NOSPACE, ELE_OP::OPSPACE);
148  op->doWorkRhsHook = [&](DataOperator *op_ptr, int side,
149  EntityType type,
150  EntitiesFieldData::EntData &data) {
151  MoFEMFunctionBegin;
152 
153  MOFEM_LOG("SELF", Sev::noisy)
154  << "parent_coords in op "
155  << static_cast<ELE_OP *>(op_ptr)->getCoordsAtGaussPts();
156 
157  parent_coords = static_cast<ELE_OP *>(op_ptr)->getCoordsAtGaussPts();
158  MoFEMFunctionReturn(0);
159  };
160  parent_fe.getOpPtrVector().push_back(op);
161 
162  MOFEM_LOG("SELF", Sev::noisy) << "fe name " << this->getFEName();
163  CHKERR this->loopParent(this->getFEName(), &parent_fe);
164  MOFEM_LOG("SELF", Sev::noisy) << "parent_coords " << parent_coords;
165 
166  MatrixDouble child_coords = this->getCoordsAtGaussPts();
167  MOFEM_LOG("SELF", Sev::noisy) << "child_coords " << child_coords;
168 
169  child_coords -= parent_coords;
170 
171  MOFEM_LOG("SELF", Sev::noisy) << "Corrds diffs" << child_coords;
172 
173  double n = 0;
174  for (auto d : child_coords.data())
175  n += d * d;
176 
177  if (sqrt(n) > 1e-12)
178  SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
179  "Parent and child global coords at integration points are "
180  "diffrent norm = %3.2e",
181  sqrt(n));
182 
183  MoFEMFunctionReturn(0);
184  }
185 };
186 
187 //! [Run programme]
188 MoFEMErrorCode AtomTest::runProblem() {
189  MoFEMFunctionBegin;
190  CHKERR readMesh();
191  CHKERR setupProblem();
192  CHKERR assembleSystem();
193  CHKERR solveSystem();
194 
195  auto test_bit_child = [](FEMethod *fe_ptr) {
196  const auto &bit = fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel();
197  MOFEM_LOG("SELF", Sev::noisy) << bit << " " << bit.test(0);
198  return bit.test(1);
199  };
200 
201  CHKERR checkResults(test_bit_child);
202  CHKERR refineResults();
203  MoFEMFunctionReturn(0);
204 }
205 //! [Run programme]
206 
207 //! [Read mesh]
208 MoFEMErrorCode AtomTest::readMesh() {
209  MoFEMFunctionBegin;
210 
211  CHKERR mField.getInterface(simpleInterface);
212  CHKERR simpleInterface->getOptions();
213  CHKERR simpleInterface->loadFile();
214 
215  MOFEM_LOG("WORLD", Sev::verbose) << "Dim " << simpleInterface->getDim();
216 
217  auto bit_level0 = simpleInterface->getBitRefLevel();
218 
219  auto &moab = mField.get_moab();
220 
221  auto refine_mesh = [&](auto bit_level1) {
222  MoFEMFunctionBegin;
223 
224  auto refine = mField.getInterface<MeshRefinement>();
225 
226  auto meshset_level0_ptr = get_temp_meshset_ptr(moab);
227  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByRefLevel(
228  bit_level0, BitRefLevel().set(), *meshset_level0_ptr);
229 
230  // random mesh refinement
231  auto meshset_ref_edges_ptr = get_temp_meshset_ptr(moab);
232  Range edges_to_refine;
233  CHKERR moab.get_entities_by_type(*meshset_level0_ptr, MBEDGE,
234  edges_to_refine);
235  int ii = 0;
236  for (Range::iterator eit = edges_to_refine.begin();
237  eit != edges_to_refine.end(); eit++, ii++) {
238  int numb = ii % 2;
239  if (numb == 0) {
240  CHKERR moab.add_entities(*meshset_ref_edges_ptr, &*eit, 1);
241  }
242  }
243  CHKERR refine->addVerticesInTheMiddleOfEdges(*meshset_ref_edges_ptr,
244  bit_level1, false, VERBOSE);
245  if (simpleInterface->getDim() == 3) {
246  CHKERR refine->refineTets(*meshset_level0_ptr, bit_level1, VERBOSE);
247  } else if (simpleInterface->getDim() == 2) {
248  CHKERR refine->refineTris(*meshset_level0_ptr, bit_level1, VERBOSE);
249  } else {
250  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
251  "Dimension not handled by test");
252  }
253 
254  MoFEMFunctionReturn(0);
255  };
256 
257  BitRefLevel bit_level1;
258  bit_level1.set(1);
259  CHKERR refine_mesh(bit_level1);
260  simpleInterface->getBitRefLevel() = BitRefLevel().set();
261  simpleInterface->getBitRefLevelMask() = BitRefLevel().set();
262 
263  MoFEMFunctionReturn(0);
264 }
265 //! [Read mesh]
266 
267 //! [Set up problem]
268 MoFEMErrorCode AtomTest::setupProblem() {
269  MoFEMFunctionBegin;
270  // Add field
271  CHKERR simpleInterface->addDomainField(FIELD_NAME, H1,
272  AINSWORTH_LEGENDRE_BASE, FIELD_DIM);
273  CHKERR simpleInterface->addBoundaryField(FIELD_NAME, H1,
274  AINSWORTH_LEGENDRE_BASE, FIELD_DIM);
275  constexpr int order = 4;
276  CHKERR simpleInterface->setFieldOrder(FIELD_NAME, order);
277  CHKERR simpleInterface->setUp();
278 
279  CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
280  simpleInterface->getProblemName(), FIELD_NAME, BitRefLevel().set(0),
281  BitRefLevel().set(0));
282 
283  MoFEMFunctionReturn(0);
284 }
285 //! [Set up problem]
286 
287 boost::shared_ptr<DomainEle> domainChildLhs, domainChildRhs;
288 
289 //! [Push operators to pipeline]
290 MoFEMErrorCode AtomTest::assembleSystem() {
291  MoFEMFunctionBegin;
292  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
293 
294  auto rule = [](int, int, int p) -> int { return 2 * p; };
295 
296  CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
297  CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
298 
299  auto test_bit_parent = [](FEMethod *fe_ptr) {
300  const auto &bit = fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel();
301  MOFEM_LOG("SELF", Sev::noisy) << bit << " " << bit.test(0);
302  return bit.test(0);
303  };
304 
305  pipeline_mng->getDomainLhsFE()->exeTestHook = test_bit_parent;
306  pipeline_mng->getDomainRhsFE()->exeTestHook = test_bit_parent;
307 
308  auto beta = [](const double, const double, const double) { return 1; };
309 
310  // Make aliased shared pointer, and create child element
311  domainChildLhs = boost::make_shared<DomainEle>(mField);
312  domainChildLhs->getRuleHook = rule;
313  domainChildLhs->getOpPtrVector().push_back(
314  new OpDomainMass(FIELD_NAME, FIELD_NAME, beta));
315 
316  domainChildRhs = boost::make_shared<DomainEle>(mField);
317  domainChildLhs->getRuleHook = rule;
318  domainChildRhs->getOpPtrVector().push_back(
319  new OpDomainSource(FIELD_NAME, approxFunction));
320 
321  auto parent_op_lhs = new DomainEleOp(NOSPACE, DomainEleOp::OPSPACE);
322  parent_op_lhs->doWorkRhsHook = [&](DataOperator *op_ptr, int side,
323  EntityType type,
324  EntitiesFieldData::EntData &data) {
325  auto domain_op = static_cast<DomainEleOp *>(op_ptr);
326  MoFEMFunctionBegin;
327 
328  MOFEM_LOG("SELF", Sev::noisy) << "LHS Pipeline FE";
329 
330  if (!domainChildLhs)
331  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "FE not allocated");
332 
333  auto &bit =
334  domain_op->getFEMethod()->numeredEntFiniteElementPtr->getBitRefLevel();
335  if (bit == BitRefLevel().set(0)) {
336  CHKERR domain_op->loopChildren(domain_op->getFEName(),
337  domainChildLhs.get(), VERBOSE, Sev::noisy);
338  } else {
339  CHKERR domain_op->loopThis(domain_op->getFEName(), domainChildLhs.get(),
340  VERBOSE, Sev::noisy);
341  }
342  MoFEMFunctionReturn(0);
343  };
344 
345  auto parent_op_rhs = new DomainEleOp(NOSPACE, DomainEleOp::OPSPACE);
346  parent_op_rhs->doWorkRhsHook = [&](DataOperator *op_ptr, int side,
347  EntityType type,
348  EntitiesFieldData::EntData &data) {
349  auto domain_op = static_cast<DomainEleOp *>(op_ptr);
350  MoFEMFunctionBegin;
351 
352  MOFEM_LOG("SELF", Sev::noisy) << "RHS Pipeline FE";
353 
354  if (!domainChildRhs)
355  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "FE not allocated");
356 
357  auto &bit =
358  domain_op->getFEMethod()->numeredEntFiniteElementPtr->getBitRefLevel();
359  if (bit == BitRefLevel().set(0)) {
360  CHKERR domain_op->loopChildren(domain_op->getFEName(),
361  domainChildRhs.get(), VERBOSE, Sev::noisy);
362  } else if ((bit & BitRefLevel().set(0)).any()) {
363  CHKERR domain_op->loopThis(domain_op->getFEName(), domainChildRhs.get(),
364  VERBOSE, Sev::noisy);
365  }
366  MoFEMFunctionReturn(0);
367  };
368 
369  pipeline_mng->getOpDomainLhsPipeline().push_back(parent_op_lhs);
370  pipeline_mng->getOpDomainRhsPipeline().push_back(parent_op_rhs);
371 
372  MoFEMFunctionReturn(0);
373 }
374 //! [Push operators to pipeline]
375 
376 //! [Solve]
377 MoFEMErrorCode AtomTest::solveSystem() {
378  MoFEMFunctionBegin;
379  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
380 
381  MOFEM_LOG("WORLD", Sev::inform) << "Solve problem";
382 
383  auto solver = pipeline_mng->createKSP();
384  CHKERR KSPSetFromOptions(solver);
385  CHKERR KSPSetUp(solver);
386 
387  auto dm = simpleInterface->getDM();
388  auto D = createDMVector(dm);
389  auto F = vectorDuplicate(D);
390 
391  CHKERR KSPSolve(solver, F, D);
392  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
393  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
394  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
395  MoFEMFunctionReturn(0);
396 }
397 
398 //! [Solve]
399 MoFEMErrorCode AtomTest::refineResults() {
400  MoFEMFunctionBegin;
401 
402  auto &moab = mField.get_moab();
403 
404  auto bit_level0 = BitRefLevel().set(0);
405  auto bit_level1 = BitRefLevel().set(1);
406  auto bit_level2 = BitRefLevel().set(2);
407 
408  auto refine_mesh = [&]() {
409  MoFEMFunctionBegin;
410 
411  auto refine = mField.getInterface<MeshRefinement>();
412 
413  auto meshset_level1_ptr = get_temp_meshset_ptr(moab);
414  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByDimAndRefLevel(
415  bit_level1, BitRefLevel().set(), simpleInterface->getDim(),
416  *meshset_level1_ptr);
417 
418  // random mesh refinement
419  auto meshset_ref_edges_ptr = get_temp_meshset_ptr(moab);
420  Range edges_to_refine;
421  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
422  bit_level1, BitRefLevel().set(), MBEDGE, edges_to_refine);
423 
424  CHKERR refine->addVerticesInTheMiddleOfEdges(edges_to_refine, bit_level2,
425  VERBOSE);
426  if (simpleInterface->getDim() == 3) {
427  CHKERR refine->refineTets(*meshset_level1_ptr, bit_level2, VERBOSE);
428  } else if (simpleInterface->getDim() == 2) {
429  CHKERR refine->refineTris(*meshset_level1_ptr, bit_level2, VERBOSE);
430  } else {
431  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
432  "Dimension not handled by test");
433  }
434 
435  Range meshsets;
436  CHKERR moab.get_entities_by_type(0, MBENTITYSET, meshsets, true);
437  for (auto m : meshsets) {
438  CHKERR mField.getInterface<BitRefManager>()
439  ->updateMeshsetByEntitiesChildren(m, bit_level2, m, MBMAXTYPE, false);
440  }
441 
442  MoFEMFunctionReturn(0);
443  };
444 
445  CHKERR refine_mesh();
446 
447  simpleInterface->getBitRefLevel() = bit_level1 | bit_level2;
448  simpleInterface->getBitRefLevelMask() = BitRefLevel().set();
449 
450  CHKERR simpleInterface->reSetUp();
451 
452  CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
453  simpleInterface->getProblemName(), FIELD_NAME, BitRefLevel().set(),
454  bit_level0 | bit_level1);
455 
456  auto project_data = [&]() {
457  MoFEMFunctionBegin;
458 
459  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
460 
461  pipeline_mng->getDomainLhsFE().reset();
462  pipeline_mng->getDomainRhsFE().reset();
463  pipeline_mng->getBoundaryRhsFE().reset();
464 
465  auto rule = [](int, int, int p) -> int { return 2 * p; };
466 
467  CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
468  CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
469  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(rule);
470 
471  auto test_bit_ref = [](FEMethod *fe_ptr) {
472  const auto &bit = fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel();
473  MOFEM_LOG("SELF", Sev::noisy) << "ref : " << bit << " " << bit.test(2);
474  return bit.test(2);
475  };
476 
477  pipeline_mng->getDomainLhsFE()->exeTestHook = test_bit_ref;
478  pipeline_mng->getDomainRhsFE()->exeTestHook = test_bit_ref;
479  pipeline_mng->getBoundaryRhsFE()->exeTestHook = test_bit_ref;
480 
481  auto beta = [](const double, const double, const double) { return 1; };
482  auto field_vals_ptr = boost::make_shared<VectorDouble>();
483 
484  auto domainParentRhs = boost::make_shared<DomainParentEle>(mField);
485  domainParentRhs->getOpPtrVector().push_back(
486  new OpCalculateScalarFieldValues(FIELD_NAME, field_vals_ptr));
487 
488  pipeline_mng->getOpDomainLhsPipeline().push_back(
489  new OpDomainMass(FIELD_NAME, FIELD_NAME, beta));
490  pipeline_mng->getOpDomainRhsPipeline().push_back(
491  new OpRunParent(domainParentRhs, bit_level2, bit_level2,
492  domainParentRhs, bit_level2, BitRefLevel().set()));
493 
494  using OpDomainTimesScalarField = FormsIntegrators<DomainEleOp>::Assembly<
495  PETSC>::LinearForm<GAUSS>::OpBaseTimesScalar<1>;
496  pipeline_mng->getOpDomainRhsPipeline().push_back(
497  new OpDomainTimesScalarField(FIELD_NAME, field_vals_ptr, beta));
498 
499  pipeline_mng->getOpDomainRhsPipeline().push_back(
500  new OpCheckGaussCoords<DomainEleOp, DomainParentEle>());
501 
502  pipeline_mng->getOpBoundaryRhsPipeline().push_back(
503  new OpCheckGaussCoords<BoundaryEleOp, BoundaryParentEle>());
504 
505  CHKERR solveSystem();
506 
507  simpleInterface->getBitRefLevel() = bit_level2;
508  simpleInterface->getBitRefLevelMask() = BitRefLevel().set();
509  CHKERR simpleInterface->reSetUp();
510 
511  CHKERR checkResults([](FEMethod *fe_ptr) { return true; });
512 
513  MoFEMFunctionReturn(0);
514  };
515 
516  CHKERR project_data();
517 
518  MoFEMFunctionReturn(0);
519 }
520 //! [Postprocess results]
521 
522 //! [Check results]
523 MoFEMErrorCode AtomTest::checkResults(
524  boost::function<bool(FEMethod *fe_method_ptr)> test_bit) {
525  MoFEMFunctionBegin;
526  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
527  pipeline_mng->getDomainLhsFE().reset();
528  pipeline_mng->getDomainRhsFE().reset();
529  pipeline_mng->getBoundaryRhsFE().reset();
530 
531  auto rule = [](int, int, int p) -> int { return 2 * p + 1; };
532  CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
533  pipeline_mng->getDomainRhsFE()->exeTestHook = test_bit;
534 
535  auto common_data_ptr = boost::make_shared<CommonData>();
536  common_data_ptr->resVec = createDMVector(simpleInterface->getDM());
537  common_data_ptr->L2Vec = createVectorMPI(
538  mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
539  common_data_ptr->approxVals = boost::make_shared<VectorDouble>();
540 
541  pipeline_mng->getOpDomainRhsPipeline().push_back(
542  new OpCalculateScalarFieldValues(FIELD_NAME,
543  common_data_ptr->approxVals));
544  pipeline_mng->getOpDomainRhsPipeline().push_back(
545  new OpError<FIELD_DIM>(common_data_ptr));
546 
547  CHKERR VecZeroEntries(common_data_ptr->L2Vec);
548  CHKERR VecZeroEntries(common_data_ptr->resVec);
549 
550  CHKERR pipeline_mng->loopFiniteElements();
551 
552  CHKERR VecAssemblyBegin(common_data_ptr->L2Vec);
553  CHKERR VecAssemblyEnd(common_data_ptr->L2Vec);
554  CHKERR VecAssemblyBegin(common_data_ptr->resVec);
555  CHKERR VecAssemblyEnd(common_data_ptr->resVec);
556  double nrm2;
557  CHKERR VecNorm(common_data_ptr->resVec, NORM_2, &nrm2);
558  const double *array;
559  CHKERR VecGetArrayRead(common_data_ptr->L2Vec, &array);
560  MOFEM_LOG_C("WORLD", Sev::inform, "Error %6.4e Vec norm %6.4e\n",
561  std::sqrt(array[0]), nrm2);
562  CHKERR VecRestoreArrayRead(common_data_ptr->L2Vec, &array);
563 
564  constexpr double eps = 1e-8;
565  if (nrm2 > eps)
566  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
567  "Not converged solution err = %6.4e", nrm2);
568  MoFEMFunctionReturn(0);
569 }
570 //! [Check results]
571 
572 int main(int argc, char *argv[]) {
573 
574  // Initialisation of MoFEM/PETSc and MOAB data structures
575  MoFEM::Core::Initialize(&argc, &argv, NULL, help);
576 
577  try {
578 
579  //! [Register MoFEM discrete manager in PETSc]
580  DMType dm_name = "DMMOFEM";
581  CHKERR DMRegister_MoFEM(dm_name);
582  //! [Register MoFEM discrete manager in PETSc
583 
584  //! [Create MoAB]
585  moab::Core mb_instance; ///< mesh database
586  moab::Interface &moab = mb_instance; ///< mesh database interface
587  //! [Create MoAB]
588 
589  //! [Create MoFEM]
590  MoFEM::Core core(moab); ///< finite element database
591  MoFEM::Interface &m_field = core; ///< finite element database insterface
592  //! [Create MoFEM]
593 
594  //! [AtomTest]
595  AtomTest ex(m_field);
596  CHKERR ex.runProblem();
597  //! [AtomTest]
598  }
599  CATCH_ERRORS;
600 
601  CHKERR MoFEM::Core::Finalize();
602 }
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
OpCheckGaussCoords::OpCheckGaussCoords
OpCheckGaussCoords()
Definition: child_and_parent.cpp:138
MoFEM::EdgeElementForcesAndSourcesCore
Edge finite element.
Definition: EdgeElementForcesAndSourcesCore.hpp:30
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
MoFEM::PipelineManager::getDomainRhsFE
boost::shared_ptr< FEMethod > & getDomainRhsFE()
Definition: PipelineManager.hpp:455
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::DataOperator
base operator to do operations at Gauss Pt. level
Definition: DataOperators.hpp:24
help
static char help[]
Definition: child_and_parent.cpp:12
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:373
MoFEM::OpRunParent
Operator to execute finite element instance on parent element. This operator is typically used to pro...
Definition: MeshProjectionDataOperators.hpp:17
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
FIELD_NAME
constexpr char FIELD_NAME[]
Definition: child_and_parent.cpp:14
AtomTest::mField
MoFEM::Interface & mField
Definition: child_and_parent.cpp:64
MoFEM::PipelineManager::loopFiniteElements
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
Definition: PipelineManager.cpp:63
ApproxFieldFunction< 1 >::operator()
double operator()(const double x, const double y, const double z)
Definition: child_and_parent.cpp:46
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::PipelineManager::getOpDomainRhsPipeline
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
Definition: PipelineManager.hpp:848
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
AtomTest::simpleInterface
Simple * simpleInterface
Definition: child_and_parent.cpp:65
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
OpCheckGaussCoords::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: child_and_parent.cpp:140
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::PipelineManager::getBoundaryRhsFE
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
Definition: PipelineManager.hpp:463
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
OpCheckGaussCoords
Definition: child_and_parent.cpp:137
AtomTest::assembleSystem
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
Definition: child_and_parent.cpp:290
MoFEM::Simple::getBitRefLevelMask
BitRefLevel & getBitRefLevelMask()
Get the BitRefLevel.
Definition: Simple.hpp:380
test_bit_child
auto test_bit_child
lambda function used to select elements on which finite element pipelines are executed.
Definition: hanging_node_approx.cpp:173
ELE_OP
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:1589
sdf.r
int r
Definition: sdf.py:8
order
constexpr int order
Definition: dg_projection.cpp:18
AtomTest::CommonData::L2Vec
SmartPetscObj< Vec > L2Vec
Definition: child_and_parent.cpp:78
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
FIELD_DIM
constexpr int FIELD_DIM
Definition: child_and_parent.cpp:15
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1502
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEM::PipelineManager::EdgeEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Definition: PipelineManager.hpp:36
AtomTest::OpError< 1 >::OpError
OpError(boost::shared_ptr< CommonData > &common_data_ptr)
Definition: child_and_parent.cpp:90
AtomTest::OpError
Operator to evaluate errors.
Definition: child_and_parent.cpp:82
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:834
VERBOSE
@ VERBOSE
Definition: definitions.h:222
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
AtomTest::OpError< 1 >::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: child_and_parent.cpp:92
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
domainChildLhs
boost::shared_ptr< DomainEle > domainChildLhs
[Set up problem]
Definition: child_and_parent.cpp:287
AtomTest::OpError< 1 >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: child_and_parent.cpp:89
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
DomainEleOp
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
Definition: child_and_parent.cpp:36
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
OpDomainTimesScalarField
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpDomainTimesScalarField
Definition: initial_diffusion.cpp:32
BoundaryEleOp
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
MoFEM::MeshRefinement
Mesh refinement interface.
Definition: MeshRefinement.hpp:26
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
MoFEM::get_temp_meshset_ptr
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1886
double
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
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::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
DomainParentEle
ElementsAndOps< SPACE_DIM >::DomianParentEle DomainParentEle
Definition: child_and_parent.cpp:35
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:99
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1214
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::FaceElementForcesAndSourcesCoreOnChildParent
Base face element used to integrate on skeleton.
Definition: FaceElementForcesAndSourcesCoreOnParent.hpp:19
MoFEM::Simple::addBoundaryField
MoFEMErrorCode addBoundaryField(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on boundary.
Definition: Simple.cpp:358
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MoFEM::PipelineManager::setDomainRhsIntegrationRule
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:605
AtomTest::CommonData
Collected data use d by operator to evaluate errors for the test.
Definition: child_and_parent.cpp:76
AtomTest::refineResults
MoFEMErrorCode refineResults()
[Solve]
Definition: child_and_parent.cpp:399
AtomTest::AtomTest
AtomTest(MoFEM::Interface &m_field)
Definition: child_and_parent.cpp:59
AtomTest::approxFunction
static ApproxFieldFunction< FIELD_DIM > approxFunction
Definition: child_and_parent.cpp:67
MoFEM::PipelineManager::getOpDomainLhsPipeline
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
Definition: PipelineManager.hpp:824
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
BiLinearForm
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:109
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
FTensor::Index< 'i', 3 >
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:604
AtomTest::solveSystem
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
Definition: child_and_parent.cpp:377
MoFEM::EdgeElementForcesAndSourcesCore::UserDataOperator
default operator for EDGE element
Definition: EdgeElementForcesAndSourcesCore.hpp:68
convert.n
n
Definition: convert.py:82
MoFEM::PipelineManager::setDomainLhsIntegrationRule
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:578
ElementsAndOps
Definition: child_and_parent.cpp:18
Range
MoFEM::PipelineManager::getDomainLhsFE
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Definition: PipelineManager.hpp:451
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
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
main
int main(int argc, char *argv[])
[Check results]
Definition: child_and_parent.cpp:572
MoFEM::EdgeElementForcesAndSourcesCoreOnChildParent
Base face element used to integrate on skeleton.
Definition: EdgeElementForcesAndSourcesCoreOnParent.hpp:19
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
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
AtomTest::checkResults
MoFEMErrorCode checkResults()
[Check results]
Definition: dg_projection.cpp:213
MoFEM::Simple::getBitRefLevel
BitRefLevel & getBitRefLevel()
Get the BitRefLevel.
Definition: Simple.hpp:373
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
MoFEM::Simple::getDim
int getDim() const
Get the problem dimension.
Definition: Simple.hpp:331
AtomTest::CommonData::resVec
SmartPetscObj< Vec > resVec
Definition: child_and_parent.cpp:79
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
MoFEM::PipelineManager::getOpBoundaryRhsPipeline
boost::ptr_deque< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
Definition: PipelineManager.hpp:896
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
AtomTest::readMesh
MoFEMErrorCode readMesh()
[Run programme]
Definition: child_and_parent.cpp:208
MoFEM::Simple::reSetUp
MoFEMErrorCode reSetUp(bool only_dm=false)
Rebuild internal MoFEM data structures.
Definition: Simple.cpp:796
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
domainChildRhs
boost::shared_ptr< DomainEle > domainChildRhs
Definition: child_and_parent.cpp:287
MoFEM::PipelineManager::setBoundaryRhsIntegrationRule
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:659
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
AtomTest::CommonData::approxVals
boost::shared_ptr< VectorDouble > approxVals
Definition: child_and_parent.cpp:77
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
BoundaryParentEle
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
Definition: child_and_parent.cpp:41
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
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:770
MoFEM::Simple::getProblemName
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:408
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
MoFEM::Simple::addDomainField
MoFEMErrorCode addDomainField(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:264
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
ApproxFieldFunction
Definition: child_and_parent.cpp:43
F
@ F
Definition: free_surface.cpp:396
Generated by Doxygen 1.8.17 and hosted at University of Glasgow