v0.14.0
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
43
44 double thetaParam;
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
80 const char *name, DataType type,
81 Tag &tag_handle) {
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 }
99 }
100
107 MoFEMErrorCode checkError(int iter_num = 0);
108 MoFEMErrorCode refineOrder(int iter_num = 0);
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;
117
119
125 };
126 };
127
128 boost::shared_ptr<CommonData> commonDataPtr;
129 struct OpError : public DomainEleOp {
130 boost::shared_ptr<CommonData> commonDataPtr;
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]
151}
152//! [Run programme]
153
154//! [Read mesh]
161}
162//! [Read mesh]
163
164//! [Set up problem]
167
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
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
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);
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 }
203}
204//! [Set up problem]
205
206//! [Set integration rule]
209
210 auto rule = [](int, int, int p) -> int { return 2 * p; };
211
213 CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
214 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
215
217}
218//! [Set integration rule]
219
220//! [Create common data]
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>();
235}
236//! [Create common data]
237
238//! [Assemble system]
242 pipeline_mng->getDomainLhsFE().reset();
243 pipeline_mng->getDomainRhsFE().reset();
244 pipeline_mng->getOpDomainRhsPipeline().clear();
245 pipeline_mng->getOpDomainLhsPipeline().clear();
246
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));
261}
262//! [Assemble system]
263
264//! [Solve]
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);
282}
283//! [Solve]
284
285//! [Refine]
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");
334 mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
337}
338//! [Refine]
339
340//! [Solve and refine loop]
348
349 int iter_num = 1;
350 while (fabs(indicTolerance) > DBL_EPSILON &&
352 MOFEM_LOG("EXAMPLE", Sev::inform) << "Refinement iteration " << iter_num;
353
354 CHKERR refineOrder(iter_num);
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 }
367}
368//! [Solve and refine loop]
369
370//! [Check error]
374 pipeline_mng->getDomainLhsFE().reset();
375 pipeline_mng->getDomainRhsFE().reset();
376 pipeline_mng->getOpDomainRhsPipeline().clear();
377
379 {HDIV, L2});
380
381 pipeline_mng->getOpDomainRhsPipeline().push_back(
382 new OpCalculateScalarFieldValues("U", commonDataPtr->approxVals));
383 pipeline_mng->getOpDomainRhsPipeline().push_back(
385 commonDataPtr->approxValsGrad));
386 pipeline_mng->getOpDomainRhsPipeline().push_back(
387 new OpCalculateHVecVectorField<3>("Q", commonDataPtr->approxFlux));
388 pipeline_mng->getOpDomainRhsPipeline().push_back(
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
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());
442}
443//! [Check error]
444
445//! [Output results]
449 pipeline_mng->getDomainLhsFE().reset();
450
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
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
476
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());
490}
491//! [Output results]
492
493//! [OpError]
495 EntData &data) {
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();
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
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 = {
558 std::array<double, CommonData::LAST_ELEMENT> values;
559 values[0] = error_l2;
560 values[1] = error_h1;
561 values[2] = error_ind;
563 indices.data(), values.data(), ADD_VALUES);
565}
566//! [OpError]
567
568int main(int argc, char *argv[]) {
569 // Initialisation of MoFEM/PETSc and MOAB data structures
570 const char param_file[] = "param_file.petsc";
572
573 // Add logging channel for example problem
574 auto core_log = logging::core::get();
575 core_log->add_sink(
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 }
602
604}
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_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
@ 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
constexpr int order
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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, 2 > OpHdivHdiv
double sqr(double x)
static char help[]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< 2 > OpHdivU
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
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
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]
double indicTolerance
Simple * simpleInterface
Range domainEntities
MoFEMErrorCode outputResults(int iter_num=0)
[Check error]
double totErrorIndicator
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]
double maxErrorIndicator
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]
static MoFEMErrorCode getTagHandle(MoFEM::Interface &m_field, const char *name, DataType type, Tag &tag_handle)
[Source function]
MoFEMErrorCode refineOrder(int iter_num=0)
[Solve]
MoFEMErrorCode readMesh()
[Run programme]
Add operators pushing bases from local to physical configuration.
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.
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
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]