v0.13.2
Loading...
Searching...
No Matches
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
12using namespace MoFEM;
13
14static char help[] = "...\n\n";
15
17
21
25 GAUSS>::OpMixDivTimesScalar<2>;
28
29inline double sqr(double x) { return x * x; }
30
32
33 MixedPoisson(MoFEM::Interface &m_field) : mField(m_field) {}
35
36private:
39
45
46 //! [Analytical function]
47 static double analyticalFunction(const double x, const double y,
48 const double z) {
49 return exp(-100. * (sqr(x) + sqr(y))) * cos(M_PI * x) * cos(M_PI * y);
50 }
51 //! [Analytical function]
52
53 //! [Analytical function gradient]
54 static VectorDouble analyticalFunctionGrad(const double x, const double y,
55 const double z) {
56 VectorDouble res;
57 res.resize(2);
58 res[0] = -exp(-100. * (sqr(x) + sqr(y))) *
59 (200. * x * cos(M_PI * x) + M_PI * sin(M_PI * x)) * cos(M_PI * y);
60 res[1] = -exp(-100. * (sqr(x) + sqr(y))) *
61 (200. * y * cos(M_PI * y) + M_PI * sin(M_PI * y)) * cos(M_PI * x);
62 return res;
63 }
64 //! [Analytical function gradient]
65
66 //! [Source function]
67 static double sourceFunction(const double x, const double y, const double z) {
68 return -exp(-100. * (sqr(x) + sqr(y))) *
69 (400. * M_PI *
70 (x * cos(M_PI * y) * sin(M_PI * x) +
71 y * cos(M_PI * x) * sin(M_PI * y)) +
72 2. * (20000. * (sqr(x) + sqr(y)) - 200. - sqr(M_PI)) *
73 cos(M_PI * x) * cos(M_PI * y));
74 }
75 //! [Source function]
76
78 const char *name, DataType type,
79 Tag &tag_handle) {
81 int int_val = 0;
82 double double_val = 0;
83 switch (type) {
84 case MB_TYPE_INTEGER:
85 CHKERR m_field.get_moab().tag_get_handle(
86 name, 1, type, tag_handle, MB_TAG_CREAT | MB_TAG_SPARSE, &int_val);
87 break;
88 case MB_TYPE_DOUBLE:
89 CHKERR m_field.get_moab().tag_get_handle(
90 name, 1, type, tag_handle, MB_TAG_CREAT | MB_TAG_SPARSE, &double_val);
91 break;
92 default:
93 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
94 "Wrong data type %d for tag", type);
95 }
97 }
98
105 MoFEMErrorCode checkError(int iter_num = 0);
108 MoFEMErrorCode outputResults(int iter_num = 0);
109
110 struct CommonData {
111 boost::shared_ptr<VectorDouble> approxVals;
112 boost::shared_ptr<MatrixDouble> approxValsGrad;
113 boost::shared_ptr<MatrixDouble> approxFlux;
115
122 };
123 };
124
125 boost::shared_ptr<CommonData> commonDataPtr;
126 struct OpError : public DomainEleOp {
127 boost::shared_ptr<CommonData> commonDataPtr;
129 OpError(boost::shared_ptr<CommonData> &common_data_ptr,
130 MoFEM::Interface &m_field)
131 : DomainEleOp("U", OPROW), commonDataPtr(common_data_ptr),
132 mField(m_field) {
133 std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
134 doEntities[MBTRI] = doEntities[MBQUAD] = true;
135 }
136 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
137 };
138};
139
140//! [Run programme]
149}
150//! [Run programme]
151
152//! [Read mesh]
159}
160//! [Read mesh]
161
162//! [Set up problem]
165 // Note that in 2D case HDIV and HCURL spaces are isomorphic, and therefore
166 // only base for HCURL has been implemented in 2D. Base vectors for HDIV space
167 // are be obtained after rotation of HCURL base vectors by a right angle
169 1);
170 // We use AINSWORTH_LEGENDRE_BASE since DEMKOWICZ_JACOBI_BASE for triangle
171 // is not yet implemented for L2 space. For quads DEMKOWICZ_JACOBI_BASE and
172 // AINSWORTH_LEGENDRE_BASE are construcreed in the same way
174
175 refIterNum = 0;
176 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-ref_iter_num", &refIterNum,
177 PETSC_NULL);
178 baseOrder = 2;
179 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-base_order", &baseOrder,
180 PETSC_NULL);
184
185 CHKERR mField.get_moab().get_entities_by_dimension(0, 2, domainEntities,
186 false);
187 Tag th_order;
188 CHKERR getTagHandle(mField, "ORDER", MB_TYPE_INTEGER, th_order);
189 for (auto ent : domainEntities) {
190 CHKERR mField.get_moab().tag_set_data(th_order, &ent, 1, &baseOrder);
191 }
193}
194//! [Set up problem]
195
196//! [Set integration rule]
199
200 auto rule = [](int, int, int p) -> int { return 2 * p + 1; };
201
203 CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
204 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
205
207}
208//! [Set integration rule]
209
210//! [Create common data]
213 commonDataPtr = boost::make_shared<CommonData>();
214 PetscInt ghosts[4] = {0, 1, 2, 3};
215 if (!mField.get_comm_rank())
216 commonDataPtr->petscVec =
217 createGhostVector(mField.get_comm(), 4, 4, 0, ghosts);
218 else
219 commonDataPtr->petscVec =
220 createGhostVector(mField.get_comm(), 0, 4, 4, ghosts);
221 commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
222 commonDataPtr->approxValsGrad = boost::make_shared<MatrixDouble>();
223 commonDataPtr->approxFlux = boost::make_shared<MatrixDouble>();
225}
226//! [Create common data]
227
228//! [Assemble system]
232 pipeline_mng->getDomainLhsFE().reset();
233 pipeline_mng->getDomainRhsFE().reset();
234 pipeline_mng->getOpDomainRhsPipeline().clear();
235 pipeline_mng->getOpDomainLhsPipeline().clear();
236
237 auto det_ptr = boost::make_shared<VectorDouble>();
238 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
239 auto jac_ptr = boost::make_shared<MatrixDouble>();
240 pipeline_mng->getOpDomainLhsPipeline().push_back(
241 new OpCalculateHOJacForFace(jac_ptr));
242 pipeline_mng->getOpDomainLhsPipeline().push_back(
243 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
244 pipeline_mng->getOpDomainLhsPipeline().push_back(new OpMakeHdivFromHcurl());
245 pipeline_mng->getOpDomainLhsPipeline().push_back(
247 pipeline_mng->getOpDomainLhsPipeline().push_back(
248 new OpSetInvJacHcurlFace(inv_jac_ptr));
249 pipeline_mng->getOpDomainLhsPipeline().push_back(new OpSetHOWeightsOnFace());
250
251 auto beta = [](const double, const double, const double) { return 1; };
252 pipeline_mng->getOpDomainLhsPipeline().push_back(
253 new OpHdivHdiv("FLUX", "FLUX", beta));
254 auto unity = []() { return 1; };
255 pipeline_mng->getOpDomainLhsPipeline().push_back(
256 new OpHdivU("FLUX", "U", unity, true));
257 auto source = [&](const double x, const double y, const double z) {
258 return -sourceFunction(x, y, z);
259 };
260 pipeline_mng->getOpDomainRhsPipeline().push_back(
261 new OpDomainSource("U", source));
263}
264//! [Assemble system]
265
266//! [Solve]
270 auto solver = pipeline_mng->createKSP();
271 CHKERR KSPSetFromOptions(solver);
272 CHKERR KSPSetUp(solver);
273
274 auto dm = simpleInterface->getDM();
275 auto D = createDMVector(dm);
276 auto F = vectorDuplicate(D);
277 CHKERR VecZeroEntries(D);
278
279 CHKERR KSPSolve(solver, F, D);
280 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
281 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
282 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
284}
285//! [Solve]
286
287//! [Refine]
290 Tag th_error_ind, th_order;
291 CHKERR getTagHandle(mField, "ERROR_INDICATOR", MB_TYPE_DOUBLE, th_error_ind);
292 CHKERR getTagHandle(mField, "ORDER", MB_TYPE_INTEGER, th_order);
293
294 std::vector<Range> refinement_levels;
295 refinement_levels.resize(refIterNum + 1);
296 for (auto ent : domainEntities) {
297 double err_indic = 0;
298 CHKERR mField.get_moab().tag_get_data(th_error_ind, &ent, 1, &err_indic);
299 int order, new_order;
300 CHKERR mField.get_moab().tag_get_data(th_order, &ent, 1, &order);
301 new_order = order + 1;
302 Range refined_ents;
303 if (err_indic > errorIndicatorIntegral / totalElementNumber) {
304 refined_ents.insert(ent);
305 Range adj;
306 CHKERR mField.get_moab().get_adjacencies(&ent, 1, 1, false, adj,
307 moab::Interface::UNION);
308 refined_ents.merge(adj);
309 refinement_levels[new_order - baseOrder].merge(refined_ents);
310 CHKERR mField.get_moab().tag_set_data(th_order, &ent, 1, &new_order);
311 }
312 }
313
314 for (int ll = 1; ll < refinement_levels.size(); ll++) {
315 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
316 refinement_levels[ll]);
317 CHKERR mField.set_field_order(refinement_levels[ll], "FLUX",
318 baseOrder + ll);
319 CHKERR mField.set_field_order(refinement_levels[ll], "U",
320 baseOrder + ll - 1);
321 }
322
323 CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities("FLUX");
324 CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities("U");
328 mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
331}
332//! [Refine]
333
334//! [Solve and refine loop]
341
342 for (int ii = 0; ii < refIterNum; ii++) {
344
347 CHKERR checkError(ii + 1);
348 CHKERR outputResults(ii + 1);
349 }
351}
352//! [Solve and refine loop]
353
354//! [Check error]
358 pipeline_mng->getDomainLhsFE().reset();
359 pipeline_mng->getDomainRhsFE().reset();
360 pipeline_mng->getOpDomainRhsPipeline().clear();
361
362 auto det_ptr = boost::make_shared<VectorDouble>();
363 auto jac_ptr = boost::make_shared<MatrixDouble>();
364 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
365
366 pipeline_mng->getOpDomainRhsPipeline().push_back(
367 new OpCalculateHOJacForFace(jac_ptr));
368 pipeline_mng->getOpDomainRhsPipeline().push_back(
369 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
370 pipeline_mng->getOpDomainRhsPipeline().push_back(new OpMakeHdivFromHcurl());
371 pipeline_mng->getOpDomainRhsPipeline().push_back(
373 pipeline_mng->getOpDomainRhsPipeline().push_back(
374 new OpSetInvJacHcurlFace(inv_jac_ptr));
375 pipeline_mng->getOpDomainRhsPipeline().push_back(
376 new OpSetInvJacL2ForFace(inv_jac_ptr));
377 pipeline_mng->getOpDomainRhsPipeline().push_back(new OpSetHOWeightsOnFace());
378
379 pipeline_mng->getOpDomainRhsPipeline().push_back(
380 new OpCalculateScalarFieldValues("U", commonDataPtr->approxVals));
381 pipeline_mng->getOpDomainRhsPipeline().push_back(
383 commonDataPtr->approxValsGrad));
384 pipeline_mng->getOpDomainRhsPipeline().push_back(
385 new OpCalculateHVecVectorField<3>("FLUX", commonDataPtr->approxFlux));
386 pipeline_mng->getOpDomainRhsPipeline().push_back(
388
389 CHKERR VecZeroEntries(commonDataPtr->petscVec);
390 CHKERR VecGhostUpdateBegin(commonDataPtr->petscVec, INSERT_VALUES,
391 SCATTER_FORWARD);
392 CHKERR VecGhostUpdateEnd(commonDataPtr->petscVec, INSERT_VALUES,
393 SCATTER_FORWARD);
394 CHKERR pipeline_mng->loopFiniteElements();
395 CHKERR VecAssemblyBegin(commonDataPtr->petscVec);
396 CHKERR VecAssemblyEnd(commonDataPtr->petscVec);
397 CHKERR VecGhostUpdateBegin(commonDataPtr->petscVec, ADD_VALUES,
398 SCATTER_REVERSE);
399 CHKERR VecGhostUpdateEnd(commonDataPtr->petscVec, ADD_VALUES,
400 SCATTER_REVERSE);
401 CHKERR VecGhostUpdateBegin(commonDataPtr->petscVec, INSERT_VALUES,
402 SCATTER_FORWARD);
403 CHKERR VecGhostUpdateEnd(commonDataPtr->petscVec, INSERT_VALUES,
404 SCATTER_FORWARD);
405 const double *array;
406 CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
407 MOFEM_LOG("EXAMPLE", Sev::inform)
408 << "Global error L2 norm: " << std::sqrt(array[0]);
409 MOFEM_LOG("EXAMPLE", Sev::inform)
410 << "Global error H1 seminorm: " << std::sqrt(array[1]);
411 MOFEM_LOG("EXAMPLE", Sev::inform)
412 << "Global error indicator: " << std::sqrt(array[2]);
413 MOFEM_LOG("EXAMPLE", Sev::inform)
414 << "Total number of elements: " << (int)array[3];
415
416 errorIndicatorIntegral = array[2];
417 totalElementNumber = (int)array[3];
418 CHKERR VecRestoreArrayRead(commonDataPtr->petscVec, &array);
419
420 std::vector<Tag> tag_handles;
421 tag_handles.resize(4);
422 CHKERR getTagHandle(mField, "ERROR_L2_NORM", MB_TYPE_DOUBLE, tag_handles[0]);
423 CHKERR getTagHandle(mField, "ERROR_H1_SEMINORM", MB_TYPE_DOUBLE,
424 tag_handles[1]);
425 CHKERR getTagHandle(mField, "ERROR_INDICATOR", MB_TYPE_DOUBLE,
426 tag_handles[2]);
427 CHKERR getTagHandle(mField, "ORDER", MB_TYPE_INTEGER, tag_handles[3]);
428
429 ParallelComm *pcomm =
430 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
431 if (pcomm == NULL)
432 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "Communicator not set");
433
434 tag_handles.push_back(pcomm->part_tag());
435 std::ostringstream strm;
436 strm << "error_" << iter_num << ".h5m";
437 CHKERR mField.get_moab().write_file(strm.str().c_str(), "MOAB",
438 "PARALLEL=WRITE_PART", 0, 0,
439 tag_handles.data(), tag_handles.size());
441}
442//! [Check error]
443
444//! [Output results]
448 pipeline_mng->getDomainLhsFE().reset();
449
451
452 auto post_proc_fe =
453 boost::make_shared<PostProcEle>(mField);
454
455 auto det_ptr = boost::make_shared<VectorDouble>();
456 auto jac_ptr = boost::make_shared<MatrixDouble>();
457 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
458
459 post_proc_fe->getOpPtrVector().push_back(
460 new OpCalculateHOJacForFace(jac_ptr));
461 post_proc_fe->getOpPtrVector().push_back(
462 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
463 post_proc_fe->getOpPtrVector().push_back(new OpMakeHdivFromHcurl());
464 post_proc_fe->getOpPtrVector().push_back(
466 post_proc_fe->getOpPtrVector().push_back(
467 new OpSetInvJacHcurlFace(inv_jac_ptr));
468
469 auto u_ptr = boost::make_shared<VectorDouble>();
470 auto flux_ptr = boost::make_shared<MatrixDouble>();
471 post_proc_fe->getOpPtrVector().push_back(
472 new OpCalculateScalarFieldValues("U", u_ptr));
473 post_proc_fe->getOpPtrVector().push_back(
474 new OpCalculateHVecVectorField<3>("FLUX", flux_ptr));
475
477
478 post_proc_fe->getOpPtrVector().push_back(
479
480 new OpPPMap(post_proc_fe->getPostProcMesh(),
481 post_proc_fe->getMapGaussPts(),
482
483 OpPPMap::DataMapVec{{"U", u_ptr}},
484
485 OpPPMap::DataMapMat{{"FLUX", flux_ptr}},
486
488
490
491 )
492
493 );
494
495 pipeline_mng->getDomainRhsFE() = post_proc_fe;
496 CHKERR pipeline_mng->loopFiniteElements();
497
498 std::ostringstream strm;
499 strm << "out_" << iter_num << ".h5m";
500 CHKERR post_proc_fe->writeFile(strm.str().c_str());
502}
503//! [Output results]
504
505//! [OpError]
507 EntData &data) {
509 const int nb_integration_pts = getGaussPts().size2();
510 const double area = getMeasure();
511 auto t_w = getFTensor0IntegrationWeight();
512 auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
513 auto t_val_grad = getFTensor1FromMat<2>(*(commonDataPtr->approxValsGrad));
514 auto t_flux = getFTensor1FromMat<3>(*(commonDataPtr->approxFlux));
515 auto t_coords = getFTensor1CoordsAtGaussPts();
518
519 double error_l2 = 0;
520 double error_h1 = 0;
521 double error_ind = 0;
522 for (int gg = 0; gg != nb_integration_pts; ++gg) {
523 const double alpha = t_w * area;
524
525 double diff = t_val - MixedPoisson::analyticalFunction(
526 t_coords(0), t_coords(1), t_coords(2));
527 error_l2 += alpha * sqr(diff);
528
530 t_coords(0), t_coords(1), t_coords(2));
531 auto t_fun_grad = getFTensor1FromArray<2, 2>(vec);
532 t_diff(i) = t_val_grad(i) - t_fun_grad(i);
533 error_h1 += alpha * t_diff(i) * t_diff(i);
534
535 t_diff(i) = t_val_grad(i) - t_flux(i);
536 error_ind += alpha * t_diff(i) * t_diff(i);
537
538 ++t_w;
539 ++t_val;
540 ++t_val_grad;
541 ++t_flux;
542 ++t_coords;
543 }
544
545 const EntityHandle ent = getFEEntityHandle();
546 Tag th_error_l2, th_error_h1, th_error_ind;
547 CHKERR MixedPoisson::getTagHandle(mField, "ERROR_L2_NORM", MB_TYPE_DOUBLE,
548 th_error_l2);
549 CHKERR MixedPoisson::getTagHandle(mField, "ERROR_H1_SEMINORM", MB_TYPE_DOUBLE,
550 th_error_h1);
551 CHKERR MixedPoisson::getTagHandle(mField, "ERROR_INDICATOR", MB_TYPE_DOUBLE,
552 th_error_ind);
553 CHKERR mField.get_moab().tag_set_data(th_error_l2, &ent, 1, &error_l2);
554 CHKERR mField.get_moab().tag_set_data(th_error_h1, &ent, 1, &error_h1);
555 CHKERR mField.get_moab().tag_set_data(th_error_ind, &ent, 1, &error_ind);
556
557 int index = CommonData::ERROR_L2_NORM;
558 constexpr std::array<int, 4> indices = {
561 std::array<double, 4> values;
562 values[0] = error_l2;
563 values[1] = error_h1;
564 values[2] = error_ind;
565 values[3] = 1.;
566 CHKERR VecSetValues(commonDataPtr->petscVec, 4, indices.data(), values.data(),
567 ADD_VALUES);
569}
570//! [OpError]
571
572int main(int argc, char *argv[]) {
573 // Initialisation of MoFEM/PETSc and MOAB data structures
574 const char param_file[] = "param_file.petsc";
576
577 // Add logging channel for example problem
578 auto core_log = logging::core::get();
579 core_log->add_sink(
581 LogManager::setLog("EXAMPLE");
582 MOFEM_LOG_TAG("EXAMPLE", "MixedPoisson");
583
584 try {
585 //! [Register MoFEM discrete manager in PETSc]
586 DMType dm_name = "DMMOFEM";
587 CHKERR DMRegister_MoFEM(dm_name);
588 //! [Register MoFEM discrete manager in PETSc
589
590 //! [Create MoAB]
591 moab::Core mb_instance; ///< mesh database
592 moab::Interface &moab = mb_instance; ///< mesh database interface
593 //! [Create MoAB]
594
595 //! [Create MoFEM]
596 MoFEM::Core core(moab); ///< finite element database
597 MoFEM::Interface &m_field = core; ///< finite element database interface
598 //! [Create MoFEM]
599
600 //! [MixedPoisson]
601 MixedPoisson ex(m_field);
602 CHKERR ex.runProblem();
603 //! [MixedPoisson]
604 }
606
608}
static Index< 'p', 3 > p
std::string param_file
int main()
Definition: adol-c_atom.cpp:46
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ HCURL
field with continuous tangents
Definition: definitions.h:86
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
@ 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
PetscErrorCode DMSetUp_MoFEM(DM dm)
Definition: DMMoFEM.cpp:1267
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.
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
FTensor::Index< 'i', SPACE_DIM > i
double D
double sqr(double x)
static char help[]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< 2 > OpHdivU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, 3 > OpHdivHdiv
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpDomainSource
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
auto createGhostVector(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[])
Create smart ghost vector.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
boost::shared_ptr< MatrixDouble > approxFlux
SmartPetscObj< Vec > petscVec
boost::shared_ptr< VectorDouble > approxVals
boost::shared_ptr< MatrixDouble > approxValsGrad
MoFEM::Interface & mField
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Output results]
OpError(boost::shared_ptr< CommonData > &common_data_ptr, MoFEM::Interface &m_field)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode solveRefineLoop()
[Refine]
MoFEMErrorCode assembleSystem()
[Create common data]
Simple * simpleInterface
Range domainEntities
MoFEMErrorCode outputResults(int iter_num=0)
[Check error]
MoFEMErrorCode refineOrder()
[Solve]
MoFEMErrorCode createCommonData()
[Set integration rule]
MoFEMErrorCode checkError(int iter_num=0)
[Solve and refine loop]
static double sourceFunction(const double x, const double y, const double z)
[Analytical function gradient]
MoFEMErrorCode solveSystem()
[Assemble system]
static VectorDouble analyticalFunctionGrad(const double x, const double y, const double z)
[Analytical function]
MoFEM::Interface & mField
boost::shared_ptr< CommonData > commonDataPtr
MixedPoisson(MoFEM::Interface &m_field)
MoFEMErrorCode setIntegrationRules()
[Set up problem]
static double analyticalFunction(const double x, const double y, const double z)
[Analytical function]
MoFEMErrorCode runProblem()
[Run programme]
MoFEMErrorCode setupProblem()
[Read mesh]
double errorIndicatorIntegral
static MoFEMErrorCode getTagHandle(MoFEM::Interface &m_field, const char *name, DataType type, Tag &tag_handle)
[Source function]
MoFEMErrorCode readMesh()
[Run programme]
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
@ OPROW
operator doWork function is executed on FE rows
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
Get vector field for H-div approximation.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Make Hdiv space from Hcurl space in 2d.
Post post-proc data at points from hash maps.
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
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
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:667
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:194
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
BitRefLevel & getBitRefLevel()
Get the BitRefLevel.
Definition: Simple.hpp:313
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, SPACE_DIM > OpHdivHdiv
[Linear elastic problem]