v0.14.0
mixed_poisson.cpp
Go to the documentation of this file.
1 /**
2  * \file mixed_poisson.cpp
3  * \example mixed_poisson.cpp
4  *
5  * MixedPoisson intended to show how to solve mixed formulation of the Dirichlet
6  * problem for the Poisson equation using error indicators and p-adaptivity
7  *
8  */
9 
10 #include <MoFEM.hpp>
11 
12 using namespace MoFEM;
13 
14 static char help[] = "...\n\n";
15 
16 #include <MoFEM.hpp>
17 
18 using DomainEle = PipelineManager::FaceEle;
19 using DomainEleOp = DomainEle::UserDataOperator;
20 using EntData = EntitiesFieldData::EntData;
21 
22 using OpHdivHdiv = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
23  GAUSS>::OpMass<3, 2>;
24 using OpHdivU = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
25  GAUSS>::OpMixDivTimesScalar<2>;
26 using OpDomainSource = FormsIntegrators<DomainEleOp>::Assembly<
27  PETSC>::LinearForm<GAUSS>::OpSource<1, 1>;
28 
29 inline double sqr(double x) { return x * x; }
30 
31 struct MixedPoisson {
32 
33  MixedPoisson(MoFEM::Interface &m_field) : mField(m_field) {}
34  MoFEMErrorCode runProblem();
35 
36 private:
37  MoFEM::Interface &mField;
38  Simple *simpleInterface;
39 
40  Range domainEntities;
41  double totErrorIndicator;
42  double maxErrorIndicator;
43 
44  double thetaParam;
45  double indicTolerance;
46  int initOrder;
47 
48  //! [Analytical function]
49  static double analyticalFunction(const double x, const double y,
50  const double z) {
51  return exp(-100. * (sqr(x) + sqr(y))) * cos(M_PI * x) * cos(M_PI * y);
52  }
53  //! [Analytical function]
54 
55  //! [Analytical function gradient]
56  static VectorDouble analyticalFunctionGrad(const double x, const double y,
57  const double z) {
58  VectorDouble res;
59  res.resize(2);
60  res[0] = -exp(-100. * (sqr(x) + sqr(y))) *
61  (200. * x * cos(M_PI * x) + M_PI * sin(M_PI * x)) * cos(M_PI * y);
62  res[1] = -exp(-100. * (sqr(x) + sqr(y))) *
63  (200. * y * cos(M_PI * y) + M_PI * sin(M_PI * y)) * cos(M_PI * x);
64  return res;
65  }
66  //! [Analytical function gradient]
67 
68  //! [Source function]
69  static double sourceFunction(const double x, const double y, const double z) {
70  return -exp(-100. * (sqr(x) + sqr(y))) *
71  (400. * M_PI *
72  (x * cos(M_PI * y) * sin(M_PI * x) +
73  y * cos(M_PI * x) * sin(M_PI * y)) +
74  2. * (20000. * (sqr(x) + sqr(y)) - 200. - sqr(M_PI)) *
75  cos(M_PI * x) * cos(M_PI * y));
76  }
77  //! [Source function]
78 
79  static MoFEMErrorCode getTagHandle(MoFEM::Interface &m_field,
80  const char *name, DataType type,
81  Tag &tag_handle) {
82  MoFEMFunctionBegin;
83  int int_val = 0;
84  double double_val = 0;
85  switch (type) {
86  case MB_TYPE_INTEGER:
87  CHKERR m_field.get_moab().tag_get_handle(
88  name, 1, type, tag_handle, MB_TAG_CREAT | MB_TAG_SPARSE, &int_val);
89  break;
90  case MB_TYPE_DOUBLE:
91  CHKERR m_field.get_moab().tag_get_handle(
92  name, 1, type, tag_handle, MB_TAG_CREAT | MB_TAG_SPARSE, &double_val);
93  break;
94  default:
95  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
96  "Wrong data type %d for tag", type);
97  }
98  MoFEMFunctionReturn(0);
99  }
100 
101  MoFEMErrorCode readMesh();
102  MoFEMErrorCode setupProblem();
103  MoFEMErrorCode setIntegrationRules();
104  MoFEMErrorCode createCommonData();
105  MoFEMErrorCode assembleSystem();
106  MoFEMErrorCode solveSystem();
107  MoFEMErrorCode checkError(int iter_num = 0);
108  MoFEMErrorCode refineOrder(int iter_num = 0);
109  MoFEMErrorCode solveRefineLoop();
110  MoFEMErrorCode outputResults(int iter_num = 0);
111 
112  struct CommonData {
113  boost::shared_ptr<VectorDouble> approxVals;
114  boost::shared_ptr<MatrixDouble> approxValsGrad;
115  boost::shared_ptr<MatrixDouble> approxFlux;
116  SmartPetscObj<Vec> petscVec;
117 
118  double maxErrorIndicator;
119 
120  enum VecElements {
121  ERROR_L2_NORM = 0,
122  ERROR_H1_SEMINORM,
123  ERROR_INDICATOR_TOTAL,
124  LAST_ELEMENT
125  };
126  };
127 
128  boost::shared_ptr<CommonData> commonDataPtr;
129  struct OpError : public DomainEleOp {
130  boost::shared_ptr<CommonData> commonDataPtr;
131  MoFEM::Interface &mField;
132  OpError(boost::shared_ptr<CommonData> &common_data_ptr,
133  MoFEM::Interface &m_field)
134  : DomainEleOp("U", OPROW), commonDataPtr(common_data_ptr),
135  mField(m_field) {
136  std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
137  doEntities[MBTRI] = doEntities[MBQUAD] = true;
138  }
139  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
140  };
141 };
142 
143 //! [Run programme]
144 MoFEMErrorCode MixedPoisson::runProblem() {
145  MoFEMFunctionBegin;
146  CHKERR readMesh();
147  CHKERR setupProblem();
148  CHKERR createCommonData();
149  CHKERR solveRefineLoop();
150  MoFEMFunctionReturn(0);
151 }
152 //! [Run programme]
153 
154 //! [Read mesh]
155 MoFEMErrorCode MixedPoisson::readMesh() {
156  MoFEMFunctionBegin;
157  CHKERR mField.getInterface(simpleInterface);
158  CHKERR simpleInterface->getOptions();
159  CHKERR simpleInterface->loadFile();
160  MoFEMFunctionReturn(0);
161 }
162 //! [Read mesh]
163 
164 //! [Set up problem]
165 MoFEMErrorCode MixedPoisson::setupProblem() {
166  MoFEMFunctionBegin;
167 
168  CHKERR simpleInterface->addDomainField("U", L2, AINSWORTH_LEGENDRE_BASE, 1);
169 
170  int nb_quads = 0;
171  CHKERR mField.get_moab().get_number_entities_by_type(0, MBQUAD, nb_quads);
172  auto base = AINSWORTH_LEGENDRE_BASE;
173  if (nb_quads) {
174  // AINSWORTH_LEGENDRE_BASE is not implemented for HDIV/HCURL space on quads
175  base = DEMKOWICZ_JACOBI_BASE;
176  }
177 
178  // Note that in 2D case HDIV and HCURL spaces are isomorphic, and therefore
179  // only base for HCURL has been implemented in 2D. Base vectors for HDIV space
180  // are be obtained after rotation of HCURL base vectors by a right angle
181  CHKERR simpleInterface->addDomainField("Q", HCURL, base, 1);
182 
183  thetaParam = 0.5;
184  CHKERR PetscOptionsGetReal(PETSC_NULL, "", "-theta", &thetaParam, PETSC_NULL);
185 
186  indicTolerance = 1e-3;
187  CHKERR PetscOptionsGetReal(PETSC_NULL, "", "-indic_tol", &indicTolerance,
188  PETSC_NULL);
189 
190  initOrder = 2;
191  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &initOrder, PETSC_NULL);
192  CHKERR simpleInterface->setFieldOrder("U", initOrder);
193  CHKERR simpleInterface->setFieldOrder("Q", initOrder + 1);
194  CHKERR simpleInterface->setUp();
195 
196  CHKERR mField.get_moab().get_entities_by_dimension(0, 2, domainEntities);
197  Tag th_order;
198  CHKERR getTagHandle(mField, "ORDER", MB_TYPE_INTEGER, th_order);
199  for (auto ent : domainEntities) {
200  CHKERR mField.get_moab().tag_set_data(th_order, &ent, 1, &initOrder);
201  }
202  MoFEMFunctionReturn(0);
203 }
204 //! [Set up problem]
205 
206 //! [Set integration rule]
207 MoFEMErrorCode MixedPoisson::setIntegrationRules() {
208  MoFEMFunctionBegin;
209 
210  auto rule = [](int, int, int p) -> int { return 2 * p; };
211 
212  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
213  CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
214  CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
215 
216  MoFEMFunctionReturn(0);
217 }
218 //! [Set integration rule]
219 
220 //! [Create common data]
221 MoFEMErrorCode MixedPoisson::createCommonData() {
222  MoFEMFunctionBegin;
223  commonDataPtr = boost::make_shared<CommonData>();
224  PetscInt ghosts[3] = {0, 1, 2};
225  if (!mField.get_comm_rank())
226  commonDataPtr->petscVec =
227  createGhostVector(mField.get_comm(), 3, 3, 0, ghosts);
228  else
229  commonDataPtr->petscVec =
230  createGhostVector(mField.get_comm(), 0, 3, 3, ghosts);
231  commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
232  commonDataPtr->approxValsGrad = boost::make_shared<MatrixDouble>();
233  commonDataPtr->approxFlux = boost::make_shared<MatrixDouble>();
234  MoFEMFunctionReturn(0);
235 }
236 //! [Create common data]
237 
238 //! [Assemble system]
239 MoFEMErrorCode MixedPoisson::assembleSystem() {
240  MoFEMFunctionBegin;
241  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
242  pipeline_mng->getDomainLhsFE().reset();
243  pipeline_mng->getDomainRhsFE().reset();
244  pipeline_mng->getOpDomainRhsPipeline().clear();
245  pipeline_mng->getOpDomainLhsPipeline().clear();
246 
247  CHKERR AddHOOps<2, 2, 2>::add(pipeline_mng->getOpDomainLhsPipeline(), {HDIV});
248 
249  auto beta = [](const double, const double, const double) { return 1; };
250  pipeline_mng->getOpDomainLhsPipeline().push_back(
251  new OpHdivHdiv("Q", "Q", beta));
252  auto unity = []() { return 1; };
253  pipeline_mng->getOpDomainLhsPipeline().push_back(
254  new OpHdivU("Q", "U", unity, true));
255  auto source = [&](const double x, const double y, const double z) {
256  return -sourceFunction(x, y, z);
257  };
258  pipeline_mng->getOpDomainRhsPipeline().push_back(
259  new OpDomainSource("U", source));
260  MoFEMFunctionReturn(0);
261 }
262 //! [Assemble system]
263 
264 //! [Solve]
265 MoFEMErrorCode MixedPoisson::solveSystem() {
266  MoFEMFunctionBegin;
267  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
268  auto solver = pipeline_mng->createKSP();
269  CHKERR KSPSetFromOptions(solver);
270  CHKERR KSPSetUp(solver);
271 
272  auto dm = simpleInterface->getDM();
273  auto D = createDMVector(dm);
274  auto F = vectorDuplicate(D);
275  CHKERR VecZeroEntries(D);
276 
277  CHKERR KSPSolve(solver, F, D);
278  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
279  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
280  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
281  MoFEMFunctionReturn(0);
282 }
283 //! [Solve]
284 
285 //! [Refine]
286 MoFEMErrorCode MixedPoisson::refineOrder(int iter_num) {
287  MoFEMFunctionBegin;
288  Tag th_error_ind, th_order;
289  CHKERR getTagHandle(mField, "ERROR_INDICATOR", MB_TYPE_DOUBLE, th_error_ind);
290  CHKERR getTagHandle(mField, "ORDER", MB_TYPE_INTEGER, th_order);
291 
292  std::vector<Range> refinement_levels;
293  refinement_levels.resize(iter_num + 1);
294  for (auto ent : domainEntities) {
295  double err_indic = 0;
296  CHKERR mField.get_moab().tag_get_data(th_error_ind, &ent, 1, &err_indic);
297 
298  int order, new_order;
299  CHKERR mField.get_moab().tag_get_data(th_order, &ent, 1, &order);
300  new_order = order + 1;
301  Range refined_ents;
302  if (err_indic > thetaParam * maxErrorIndicator) {
303  refined_ents.insert(ent);
304  Range adj;
305  CHKERR mField.get_moab().get_adjacencies(&ent, 1, 1, false, adj,
306  moab::Interface::UNION);
307  refined_ents.merge(adj);
308  refinement_levels[new_order - initOrder].merge(refined_ents);
309  CHKERR mField.get_moab().tag_set_data(th_order, &ent, 1, &new_order);
310  }
311  }
312 
313  for (int ll = 1; ll < refinement_levels.size(); ll++) {
314  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
315  refinement_levels[ll]);
316 
317  if (initOrder + ll > 8) {
318  MOFEM_LOG("EXAMPLE", Sev::warning)
319  << "setting approximation order higher than 8 is not currently "
320  "supported"
321  << endl;
322  } else {
323  CHKERR mField.set_field_order(refinement_levels[ll], "U", initOrder + ll);
324  CHKERR mField.set_field_order(refinement_levels[ll], "Q",
325  initOrder + ll + 1);
326  }
327  }
328 
329  CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities("Q");
330  CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities("U");
331  CHKERR mField.build_fields();
332  CHKERR mField.build_finite_elements();
333  CHKERR mField.build_adjacencies(simpleInterface->getBitRefLevel());
334  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
335  CHKERR DMSetUp_MoFEM(simpleInterface->getDM());
336  MoFEMFunctionReturn(0);
337 }
338 //! [Refine]
339 
340 //! [Solve and refine loop]
341 MoFEMErrorCode MixedPoisson::solveRefineLoop() {
342  MoFEMFunctionBegin;
343  CHKERR assembleSystem();
344  CHKERR setIntegrationRules();
345  CHKERR solveSystem();
346  CHKERR checkError();
347  CHKERR outputResults();
348 
349  int iter_num = 1;
350  while (fabs(indicTolerance) > DBL_EPSILON &&
351  totErrorIndicator > indicTolerance) {
352  MOFEM_LOG("EXAMPLE", Sev::inform) << "Refinement iteration " << iter_num;
353 
354  CHKERR refineOrder(iter_num);
355  CHKERR assembleSystem();
356  CHKERR setIntegrationRules();
357  CHKERR solveSystem();
358  CHKERR checkError(iter_num);
359  CHKERR outputResults(iter_num);
360 
361  iter_num++;
362  if (iter_num > 100)
363  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
364  "Too many refinement iterations");
365  }
366  MoFEMFunctionReturn(0);
367 }
368 //! [Solve and refine loop]
369 
370 //! [Check error]
371 MoFEMErrorCode MixedPoisson::checkError(int iter_num) {
372  MoFEMFunctionBegin;
373  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
374  pipeline_mng->getDomainLhsFE().reset();
375  pipeline_mng->getDomainRhsFE().reset();
376  pipeline_mng->getOpDomainRhsPipeline().clear();
377 
378  CHKERR AddHOOps<2, 2, 2>::add(pipeline_mng->getOpDomainRhsPipeline(),
379  {HDIV, L2});
380 
381  pipeline_mng->getOpDomainRhsPipeline().push_back(
382  new OpCalculateScalarFieldValues("U", commonDataPtr->approxVals));
383  pipeline_mng->getOpDomainRhsPipeline().push_back(
384  new OpCalculateScalarFieldGradient<2>("U",
385  commonDataPtr->approxValsGrad));
386  pipeline_mng->getOpDomainRhsPipeline().push_back(
387  new OpCalculateHVecVectorField<3>("Q", commonDataPtr->approxFlux));
388  pipeline_mng->getOpDomainRhsPipeline().push_back(
389  new OpError(commonDataPtr, mField));
390 
391  commonDataPtr->maxErrorIndicator = 0;
392  CHKERR VecZeroEntries(commonDataPtr->petscVec);
393  CHKERR pipeline_mng->loopFiniteElements();
394  CHKERR VecAssemblyBegin(commonDataPtr->petscVec);
395  CHKERR VecAssemblyEnd(commonDataPtr->petscVec);
396  CHKERR VecGhostUpdateBegin(commonDataPtr->petscVec, INSERT_VALUES,
397  SCATTER_FORWARD);
398  CHKERR VecGhostUpdateEnd(commonDataPtr->petscVec, INSERT_VALUES,
399  SCATTER_FORWARD);
400 
401  MPI_Allreduce(&commonDataPtr->maxErrorIndicator, &maxErrorIndicator, 1,
402  MPI_DOUBLE, MPI_MAX, mField.get_comm());
403 
404  const double *array;
405  CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
406  MOFEM_LOG("EXAMPLE", Sev::inform)
407  << "Global error indicator (max): " << commonDataPtr->maxErrorIndicator;
408  MOFEM_LOG("EXAMPLE", Sev::inform)
409  << "Global error indicator (total): "
410  << std::sqrt(array[CommonData::ERROR_INDICATOR_TOTAL]);
411  MOFEM_LOG("EXAMPLE", Sev::inform)
412  << "Global error L2 norm: "
413  << std::sqrt(array[CommonData::ERROR_L2_NORM]);
414  MOFEM_LOG("EXAMPLE", Sev::inform)
415  << "Global error H1 seminorm: "
416  << std::sqrt(array[CommonData::ERROR_H1_SEMINORM]);
417 
418  totErrorIndicator = std::sqrt(array[CommonData::ERROR_INDICATOR_TOTAL]);
419  CHKERR VecRestoreArrayRead(commonDataPtr->petscVec, &array);
420 
421  std::vector<Tag> tag_handles;
422  tag_handles.resize(4);
423  CHKERR getTagHandle(mField, "ERROR_L2_NORM", MB_TYPE_DOUBLE, tag_handles[0]);
424  CHKERR getTagHandle(mField, "ERROR_H1_SEMINORM", MB_TYPE_DOUBLE,
425  tag_handles[1]);
426  CHKERR getTagHandle(mField, "ERROR_INDICATOR", MB_TYPE_DOUBLE,
427  tag_handles[2]);
428  CHKERR getTagHandle(mField, "ORDER", MB_TYPE_INTEGER, tag_handles[3]);
429 
430  ParallelComm *pcomm =
431  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
432  if (pcomm == NULL)
433  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "Communicator not set");
434 
435  tag_handles.push_back(pcomm->part_tag());
436  std::ostringstream strm;
437  strm << "error_" << iter_num << ".h5m";
438  CHKERR mField.get_moab().write_file(strm.str().c_str(), "MOAB",
439  "PARALLEL=WRITE_PART", 0, 0,
440  tag_handles.data(), tag_handles.size());
441  MoFEMFunctionReturn(0);
442 }
443 //! [Check error]
444 
445 //! [Output results]
446 MoFEMErrorCode MixedPoisson::outputResults(int iter_num) {
447  MoFEMFunctionBegin;
448  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
449  pipeline_mng->getDomainLhsFE().reset();
450 
451  using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;
452 
453  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
454 
455  CHKERR AddHOOps<2, 2, 2>::add(post_proc_fe->getOpPtrVector(), {HDIV});
456 
457  auto u_ptr = boost::make_shared<VectorDouble>();
458  auto flux_ptr = boost::make_shared<MatrixDouble>();
459  post_proc_fe->getOpPtrVector().push_back(
460  new OpCalculateScalarFieldValues("U", u_ptr));
461  post_proc_fe->getOpPtrVector().push_back(
462  new OpCalculateHVecVectorField<3>("Q", flux_ptr));
463 
464  using OpPPMap = OpPostProcMapInMoab<3, 3>;
465 
466  post_proc_fe->getOpPtrVector().push_back(
467 
468  new OpPPMap(post_proc_fe->getPostProcMesh(),
469  post_proc_fe->getMapGaussPts(),
470 
471  OpPPMap::DataMapVec{{"U", u_ptr}},
472 
473  OpPPMap::DataMapMat{{"Q", flux_ptr}},
474 
475  OpPPMap::DataMapMat{},
476 
477  OpPPMap::DataMapMat{}
478 
479  )
480 
481  );
482 
483  pipeline_mng->getDomainRhsFE() = post_proc_fe;
484  CHKERR pipeline_mng->loopFiniteElements();
485 
486  std::ostringstream strm;
487  strm << "out_" << iter_num << ".h5m";
488  CHKERR post_proc_fe->writeFile(strm.str().c_str());
489  MoFEMFunctionReturn(0);
490 }
491 //! [Output results]
492 
493 //! [OpError]
494 MoFEMErrorCode MixedPoisson::OpError::doWork(int side, EntityType type,
495  EntData &data) {
496  MoFEMFunctionBegin;
497  const int nb_integration_pts = getGaussPts().size2();
498  const double area = getMeasure();
499  auto t_w = getFTensor0IntegrationWeight();
500  auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
501  auto t_val_grad = getFTensor1FromMat<2>(*(commonDataPtr->approxValsGrad));
502  auto t_flux = getFTensor1FromMat<3>(*(commonDataPtr->approxFlux));
503  auto t_coords = getFTensor1CoordsAtGaussPts();
504  FTensor::Tensor1<double, 2> t_diff;
505  FTensor::Index<'i', 2> i;
506 
507  double error_l2 = 0;
508  double error_h1 = 0;
509  double error_ind = 0;
510  for (int gg = 0; gg != nb_integration_pts; ++gg) {
511  const double alpha = t_w * area;
512 
513  double diff = t_val - MixedPoisson::analyticalFunction(
514  t_coords(0), t_coords(1), t_coords(2));
515  error_l2 += alpha * sqr(diff);
516 
517  VectorDouble vec = MixedPoisson::analyticalFunctionGrad(
518  t_coords(0), t_coords(1), t_coords(2));
519  auto t_fun_grad = getFTensor1FromArray<2, 2>(vec);
520  t_diff(i) = t_val_grad(i) - t_fun_grad(i);
521  error_h1 += alpha * t_diff(i) * t_diff(i);
522 
523  t_diff(i) = t_val_grad(i) - t_flux(i);
524  error_ind += alpha * t_diff(i) * t_diff(i);
525 
526  ++t_w;
527  ++t_val;
528  ++t_val_grad;
529  ++t_flux;
530  ++t_coords;
531  }
532 
533  const EntityHandle ent = getFEEntityHandle();
534  Tag th_error_l2, th_error_h1, th_error_ind;
535  CHKERR MixedPoisson::getTagHandle(mField, "ERROR_L2_NORM", MB_TYPE_DOUBLE,
536  th_error_l2);
537  CHKERR MixedPoisson::getTagHandle(mField, "ERROR_H1_SEMINORM", MB_TYPE_DOUBLE,
538  th_error_h1);
539  CHKERR MixedPoisson::getTagHandle(mField, "ERROR_INDICATOR", MB_TYPE_DOUBLE,
540  th_error_ind);
541 
542  double error_l2_norm = sqrt(error_l2);
543  double error_h1_seminorm = sqrt(error_h1);
544  double error_ind_local = sqrt(error_ind);
545  CHKERR mField.get_moab().tag_set_data(th_error_l2, &ent, 1, &error_l2_norm);
546  CHKERR mField.get_moab().tag_set_data(th_error_h1, &ent, 1,
547  &error_h1_seminorm);
548  CHKERR mField.get_moab().tag_set_data(th_error_ind, &ent, 1,
549  &error_ind_local);
550 
551  if (error_ind_local > commonDataPtr->maxErrorIndicator)
552  commonDataPtr->maxErrorIndicator = error_ind_local;
553 
554  int index = CommonData::ERROR_L2_NORM;
555  constexpr std::array<int, CommonData::LAST_ELEMENT> indices = {
556  CommonData::ERROR_L2_NORM, CommonData::ERROR_H1_SEMINORM,
557  CommonData::ERROR_INDICATOR_TOTAL};
558  std::array<double, CommonData::LAST_ELEMENT> values;
559  values[0] = error_l2;
560  values[1] = error_h1;
561  values[2] = error_ind;
562  CHKERR VecSetValues(commonDataPtr->petscVec, CommonData::LAST_ELEMENT,
563  indices.data(), values.data(), ADD_VALUES);
564  MoFEMFunctionReturn(0);
565 }
566 //! [OpError]
567 
568 int main(int argc, char *argv[]) {
569  // Initialisation of MoFEM/PETSc and MOAB data structures
570  const char param_file[] = "param_file.petsc";
571  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
572 
573  // Add logging channel for example problem
574  auto core_log = logging::core::get();
575  core_log->add_sink(
576  LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
577  LogManager::setLog("EXAMPLE");
578  MOFEM_LOG_TAG("EXAMPLE", "MixedPoisson");
579 
580  try {
581  //! [Register MoFEM discrete manager in PETSc]
582  DMType dm_name = "DMMOFEM";
583  CHKERR DMRegister_MoFEM(dm_name);
584  //! [Register MoFEM discrete manager in PETSc
585 
586  //! [Create MoAB]
587  moab::Core mb_instance; ///< mesh database
588  moab::Interface &moab = mb_instance; ///< mesh database interface
589  //! [Create MoAB]
590 
591  //! [Create MoFEM]
592  MoFEM::Core core(moab); ///< finite element database
593  MoFEM::Interface &m_field = core; ///< finite element database interface
594  //! [Create MoFEM]
595 
596  //! [MixedPoisson]
597  MixedPoisson ex(m_field);
598  CHKERR ex.runProblem();
599  //! [MixedPoisson]
600  }
601  CATCH_ERRORS;
602 
603  CHKERR MoFEM::Core::Finalize();
604 }
MixedPoisson::domainEntities
Range domainEntities
Definition: mixed_poisson.cpp:40
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MixedPoisson::OpError::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Output results]
Definition: mixed_poisson.cpp:494
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MoFEM::PipelineManager::getDomainRhsFE
boost::shared_ptr< FEMethod > & getDomainRhsFE()
Definition: PipelineManager.hpp:405
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MixedPoisson::CommonData::approxFlux
boost::shared_ptr< MatrixDouble > approxFlux
Definition: mixed_poisson.cpp:115
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
MixedPoisson::CommonData::petscVec
SmartPetscObj< Vec > petscVec
Definition: mixed_poisson.cpp:116
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
OpHdivHdiv
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, SPACE_DIM > OpHdivHdiv
[Linear elastic problem]
Definition: thermo_elastic.cpp:47
EntityHandle
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
MoFEM::OpPostProcMapInMoab::DataMapVec
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
Definition: PostProcBrokenMeshInMoabBase.hpp:700
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
MoFEM::PipelineManager::loopFiniteElements
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
Definition: PipelineManager.cpp:19
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MixedPoisson::analyticalFunctionGrad
static VectorDouble analyticalFunctionGrad(const double x, const double y, const double z)
[Analytical function]
Definition: mixed_poisson.cpp:56
MoFEM::PipelineManager::getOpDomainRhsPipeline
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
Definition: PipelineManager.hpp:773
MixedPoisson::CommonData
Definition: mixed_poisson.cpp:112
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM::DMSetUp_MoFEM
PetscErrorCode DMSetUp_MoFEM(DM dm)
Definition: DMMoFEM.cpp:1303
MoFEM.hpp
MixedPoisson::readMesh
MoFEMErrorCode readMesh()
[Run programme]
Definition: mixed_poisson.cpp:155
MixedPoisson::CommonData::ERROR_INDICATOR_TOTAL
@ ERROR_INDICATOR_TOTAL
Definition: mixed_poisson.cpp:123
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
MixedPoisson::CommonData::VecElements
VecElements
Definition: mixed_poisson.cpp:120
MixedPoisson::CommonData::approxVals
boost::shared_ptr< VectorDouble > approxVals
Definition: mixed_poisson.cpp:113
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
MixedPoisson::initOrder
int initOrder
Definition: mixed_poisson.cpp:46
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MixedPoisson::outputResults
MoFEMErrorCode outputResults(int iter_num=0)
[Check error]
Definition: mixed_poisson.cpp:446
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1293
MixedPoisson::CommonData::maxErrorIndicator
double maxErrorIndicator
Definition: mixed_poisson.cpp:118
MixedPoisson::OpError
Definition: mixed_poisson.cpp:129
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MixedPoisson::CommonData::ERROR_L2_NORM
@ ERROR_L2_NORM
Definition: mixed_poisson.cpp:121
MixedPoisson::mField
MoFEM::Interface & mField
Definition: mixed_poisson.cpp:37
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
MixedPoisson::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: mixed_poisson.cpp:128
MixedPoisson::CommonData::LAST_ELEMENT
@ LAST_ELEMENT
Definition: mixed_poisson.cpp:124
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::createGhostVector
auto createGhostVector(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[])
Create smart ghost vector.
Definition: PetscSmartObj.hpp:179
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MixedPoisson::OpError::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: mixed_poisson.cpp:130
MixedPoisson::sourceFunction
static double sourceFunction(const double x, const double y, const double z)
[Analytical function gradient]
Definition: mixed_poisson.cpp:69
MixedPoisson::indicTolerance
double indicTolerance
Definition: mixed_poisson.cpp:45
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MixedPoisson::refineOrder
MoFEMErrorCode refineOrder(int iter_num=0)
[Solve]
Definition: mixed_poisson.cpp:286
MixedPoisson::setupProblem
MoFEMErrorCode setupProblem()
[Read mesh]
Definition: mixed_poisson.cpp:165
MixedPoisson::setIntegrationRules
MoFEMErrorCode setIntegrationRules()
[Set up problem]
Definition: mixed_poisson.cpp:207
OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: seepage.cpp:57
MoFEM::PetscOptionsGetReal
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:152
MixedPoisson::solveSystem
MoFEMErrorCode solveSystem()
[Assemble system]
Definition: mixed_poisson.cpp:265
MoFEM::OpCalculateHVecVectorField
Get vector field for H-div approximation.
Definition: UserDataOperators.hpp:2115
MixedPoisson::OpError::OpError
OpError(boost::shared_ptr< CommonData > &common_data_ptr, MoFEM::Interface &m_field)
Definition: mixed_poisson.cpp:132
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
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
main
int main(int argc, char *argv[])
[OpError]
Definition: mixed_poisson.cpp:568
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:49
help
static char help[]
Definition: mixed_poisson.cpp:14
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
MoFEM::PipelineManager::setDomainRhsIntegrationRule
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:530
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MoFEM::PipelineManager::getOpDomainLhsPipeline
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
Definition: PipelineManager.hpp:749
MixedPoisson::MixedPoisson
MixedPoisson(MoFEM::Interface &m_field)
Definition: mixed_poisson.cpp:33
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::OpPostProcMapInMoab::DataMapMat
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Definition: PostProcBrokenMeshInMoabBase.hpp:701
MixedPoisson::OpError::mField
MoFEM::Interface & mField
Definition: mixed_poisson.cpp:131
BiLinearForm
EntData
EntitiesFieldData::EntData EntData
Definition: mixed_poisson.cpp:20
FTensor::Index< 'i', 2 >
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:413
MixedPoisson::analyticalFunction
static double analyticalFunction(const double x, const double y, const double z)
[Analytical function]
Definition: mixed_poisson.cpp:49
MoFEM::PipelineManager::setDomainLhsIntegrationRule
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:503
MixedPoisson::maxErrorIndicator
double maxErrorIndicator
Definition: mixed_poisson.cpp:42
Range
MoFEM::PipelineManager::getDomainLhsFE
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Definition: PipelineManager.hpp:401
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
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
OpHdivHdiv
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, 2 > OpHdivHdiv
Definition: mixed_poisson.cpp:23
MixedPoisson::simpleInterface
Simple * simpleInterface
Definition: mixed_poisson.cpp:38
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
OpHdivU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< 2 > OpHdivU
Definition: mixed_poisson.cpp:25
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
MixedPoisson
Definition: mixed_poisson.cpp:31
MixedPoisson::thetaParam
double thetaParam
Definition: mixed_poisson.cpp:44
MixedPoisson::totErrorIndicator
double totErrorIndicator
Definition: mixed_poisson.cpp:41
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MixedPoisson::solveRefineLoop
MoFEMErrorCode solveRefineLoop()
[Refine]
Definition: mixed_poisson.cpp:341
MixedPoisson::checkError
MoFEMErrorCode checkError(int iter_num=0)
[Solve and refine loop]
Definition: mixed_poisson.cpp:371
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
MixedPoisson::runProblem
MoFEMErrorCode runProblem()
[Run programme]
Definition: mixed_poisson.cpp:144
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
sqr
double sqr(double x)
Definition: mixed_poisson.cpp:29
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MixedPoisson::assembleSystem
MoFEMErrorCode assembleSystem()
[Create common data]
Definition: mixed_poisson.cpp:239
MoFEM::SmartPetscObj< Vec >
source
auto source
Definition: poisson_2d_dis_galerkin.cpp:62
OpDomainSource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
Definition: child_and_parent.cpp:55
MixedPoisson::CommonData::ERROR_H1_SEMINORM
@ ERROR_H1_SEMINORM
Definition: mixed_poisson.cpp:122
MixedPoisson::createCommonData
MoFEMErrorCode createCommonData()
[Set integration rule]
Definition: mixed_poisson.cpp:221
convert.int
int
Definition: convert.py:64
OpDomainSource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpDomainSource
Definition: mixed_poisson.cpp:27
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MixedPoisson::CommonData::approxValsGrad
boost::shared_ptr< MatrixDouble > approxValsGrad
Definition: mixed_poisson.cpp:114
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MixedPoisson::getTagHandle
static MoFEMErrorCode getTagHandle(MoFEM::Interface &m_field, const char *name, DataType type, Tag &tag_handle)
[Source function]
Definition: mixed_poisson.cpp:79
F
@ F
Definition: free_surface.cpp:394
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
Generated by Doxygen 1.8.17 and hosted at University of Glasgow