v0.15.0
Loading...
Searching...
No Matches
incompressible_elasticity.cpp
Go to the documentation of this file.
1/**
2 * @file incompressible_elasticity.cpp
3 * @brief Incompressible elasticity problem
4 */
5
6#include <MoFEM.hpp>
7
8using namespace MoFEM;
9
10constexpr int SPACE_DIM =
11 EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
12
13constexpr AssemblyType AT =
14 (SCHUR_ASSEMBLE) ? AssemblyType::SCHUR
15 : AssemblyType::PETSC; //< selected assembly type
16
18 IntegrationType::GAUSS; //< selected integration type
20
23using DomainEleOp = DomainEle::UserDataOperator;
26using BoundaryEleOp = BoundaryEle::UserDataOperator;
30
31
32PetscBool doEvalField;
34
37 std::pair<boost::shared_ptr<PostProcEle>,
38 boost::shared_ptr<SkinPostProcEle>>
39 pair_post_proc_fe,
40 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> ux_scatter,
41 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uy_scatter,
42 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uz_scatter,
43 std::array<double, 3> pass_field_eval_coords,
44 boost::shared_ptr<SetPtsData> pass_field_eval_data,
45 boost::shared_ptr<MatrixDouble> vec_field_ptr)
46 : dM(dm), uXScatter(ux_scatter), uYScatter(uy_scatter),
47 uZScatter(uz_scatter), fieldEvalCoords(pass_field_eval_coords),
48 fieldEvalData(pass_field_eval_data), vecFieldPtr(vec_field_ptr) {
49 postProcFe = pair_post_proc_fe.first;
50 skinPostProcFe = pair_post_proc_fe.second;
51 };
52
53 MoFEMErrorCode preProcess() { return 0; }
54 MoFEMErrorCode operator()() { return 0; }
55
58
59 MoFEM::Interface *m_field_ptr;
60 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
61 auto *simple = m_field_ptr->getInterface<Simple>();
62
63 if (doEvalField) {
65 ->evalFEAtThePoint<SPACE_DIM>(
66 fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
67 simple->getDomainFEName(), fieldEvalData,
68 m_field_ptr->get_comm_rank(), m_field_ptr->get_comm_rank(),
69 getCacheWeakPtr().lock(), MF_EXIST, QUIET);
70
71 if (vecFieldPtr->size1()) {
73 if constexpr (SPACE_DIM == 2)
74 MOFEM_LOG("SYNC", Sev::inform)
75 << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1);
76 else
77 MOFEM_LOG("SYNC", Sev::inform)
78 << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1)
79 << " U_Z: " << t_disp(2);
80 }
81 }
82 MOFEM_LOG_SEVERITY_SYNC(m_field_ptr->get_comm(), Sev::inform);
83
84 auto make_vtk = [&]() {
86 if (postProcFe) {
89 CHKERR postProcFe->writeFile("out_incomp_elasticity_" +
90 boost::lexical_cast<std::string>(ts_step) +
91 ".h5m");
92 }
93 if (skinPostProcFe) {
96 CHKERR skinPostProcFe->writeFile(
97 "out_skin_incomp_elasticity_" +
98 boost::lexical_cast<std::string>(ts_step) + ".h5m");
99 }
101 };
102
103 auto print_max_min = [&](auto &tuple, const std::string msg) {
105 CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
106 INSERT_VALUES, SCATTER_FORWARD);
107 CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
108 INSERT_VALUES, SCATTER_FORWARD);
109 double max, min;
110 CHKERR VecMax(std::get<0>(tuple), PETSC_NULLPTR, &max);
111 CHKERR VecMin(std::get<0>(tuple), PETSC_NULLPTR, &min);
112 MOFEM_LOG_C("INCOMP_ELASTICITY", Sev::inform,
113 "%s time %3.4e min %3.4e max %3.4e", msg.c_str(), ts_t, min,
114 max);
116 };
117
118 CHKERR make_vtk();
119 CHKERR print_max_min(uXScatter, "Ux");
120 CHKERR print_max_min(uYScatter, "Uy");
121 if constexpr (SPACE_DIM == 3)
122 CHKERR print_max_min(uZScatter, "Uz");
123
125 }
126
127private:
129 boost::shared_ptr<PostProcEle> postProcFe;
130 boost::shared_ptr<SkinPostProcEle> skinPostProcFe;
131 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
132 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
133 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
134 std::array<double, 3> fieldEvalCoords;
135 boost::shared_ptr<SetPtsData> fieldEvalData;
136 boost::shared_ptr<MatrixDouble> vecFieldPtr;
137};
138
139// Assemble to A matrix, by default, however, some terms are assembled only to
140// preconditioning.
141
142template <>
145 [](ForcesAndSourcesCore::UserDataOperator *op_ptr,
146 const EntitiesFieldData::EntData &row_data,
147 const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
149 op_ptr->getKSPA(), row_data, col_data, m, ADD_VALUES);
150 };
151
152/**
153 * @brief Element used to specialise assembly
154 *
155 */
157 using DomainEleOp::DomainEleOp;
158};
159
160/**
161 * @brief Specialise assembly for Stabilised matrix
162 *
163 * @tparam
164 */
165template <>
168 [](ForcesAndSourcesCore::UserDataOperator *op_ptr,
169 const EntitiesFieldData::EntData &row_data,
170 const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
172 op_ptr->getKSPB(), row_data, col_data, m, ADD_VALUES);
173 };
174//! [Specialisation for assembly]
175
176int order = 2;
178inline static double young_modulus = 100;
179inline static double poisson_ratio = 0.25;
180inline static double mu;
181inline static double lambda;
182
183PetscBool isDiscontinuousPressure = PETSC_FALSE;
184
186
187 Incompressible(MoFEM::Interface &m_field) : mField(m_field) {}
188
190
191private:
193
200
201 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
202 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
203 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
204 std::array<double, 3> fieldEvalCoords;
205};
206
207template <int DIM>
208struct OpCalculateLameStress : public ForcesAndSourcesCore::UserDataOperator {
209 OpCalculateLameStress(double m_u, boost::shared_ptr<MatrixDouble> stress_ptr,
210 boost::shared_ptr<MatrixDouble> strain_ptr,
211 boost::shared_ptr<VectorDouble> pressure_ptr)
213 stressPtr(stress_ptr), strainPtr(strain_ptr),
214 pressurePtr(pressure_ptr) {}
215
216 MoFEMErrorCode doWork(int side, EntityType type,
219 // Define Indicies
222
223 // Define Kronecker Delta
225
226 // Number of Gauss points
227 const size_t nb_gauss_pts = getGaussPts().size2();
228
229 stressPtr->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts);
232 auto t_pressure = getFTensor0FromVec(*(pressurePtr));
233
234 const double l_mu = mU;
235 // Extract matrix from data matrix
236 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
237
238 t_stress(i, j) = t_pressure * t_kd(i, j) + 2. * l_mu * t_strain(i, j);
239
240 ++t_strain;
241 ++t_stress;
242 ++t_pressure;
243 }
244
246 }
247
248private:
249 double mU;
250 boost::shared_ptr<MatrixDouble> stressPtr;
251 boost::shared_ptr<MatrixDouble> strainPtr;
252 boost::shared_ptr<VectorDouble> pressurePtr;
253};
254
255//! [Run problem]
266//! [Run problem]
267
268//! [Set up problem]
272
273 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
274 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-geom_order", &geom_order,
275 PETSC_NULLPTR);
276 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-is_discontinuous_pressure",
277 &isDiscontinuousPressure, PETSC_NULLPTR);
278
279 MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform) << "Order " << order;
280 MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform) << "Geom order " << geom_order;
281
282 // Select base
283 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
284 const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
285 PetscInt choice_base_value = AINSWORTH;
286 CHKERR PetscOptionsGetEList(PETSC_NULLPTR, PETSC_NULLPTR, "-base", list_bases,
287 LASBASETOPT, &choice_base_value, PETSC_NULLPTR);
288 int coords_dim = 3;
289 CHKERR PetscOptionsGetRealArray(PETSC_NULLPTR, PETSC_NULLPTR,
290 "-field_eval_coords", fieldEvalCoords.data(),
291 &coords_dim, &doEvalField);
292
294 switch (choice_base_value) {
295 case AINSWORTH:
297 MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
298 << "Set AINSWORTH_LEGENDRE_BASE for displacements";
299 break;
300 case DEMKOWICZ:
302 MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
303 << "Set DEMKOWICZ_JACOBI_BASE for displacements";
304 break;
305 default:
306 base = LASTBASE;
307 break;
308 }
309
310 // Note: For tets we have only H1 Ainsworth base, for Hex we have only H1
311 // Demkowicz base. We need to implement Demkowicz H1 base on tet.
312 CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
313 CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
314
315 // Adding fields related to incompressible elasticity
316 // Add displacement domain and boundary fields
317 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
318 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
319 CHKERR simple->setFieldOrder("U", order);
320
321 // Add pressure domain and boundary fields
322 // Choose either Crouzeix-Raviart element:
324 CHKERR simple->addDomainField("P", L2, base, 1);
325 CHKERR simple->setFieldOrder("P", order - 2);
326 } else {
327 // ... or Taylor-Hood element:
328 CHKERR simple->addDomainField("P", H1, base, 1);
329 CHKERR simple->setFieldOrder("P", order - 1);
330 }
331
332 // Add geometry data field
333 CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
334 CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
335
336 CHKERR simple->setUp();
337
338 auto project_ho_geometry = [&]() {
339 Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
340 return mField.loop_dofs("GEOMETRY", ent_method);
341 };
342 CHKERR project_ho_geometry();
343
345} //! [Set up problem]
346
347//! [Create common data]
350
351 auto get_options = [&]() {
353 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-young_modulus",
354 &young_modulus, PETSC_NULLPTR);
355 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-poisson_ratio",
356 &poisson_ratio, PETSC_NULLPTR);
357
358 MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
359 << "Young modulus " << young_modulus;
360 MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
361 << "Poisson_ratio " << poisson_ratio;
362
363 mu = young_modulus / (2. + 2. * poisson_ratio);
364 const double lambda_denom =
365 (1. + poisson_ratio) * (1. - 2. * poisson_ratio);
366 lambda = young_modulus * poisson_ratio / lambda_denom;
367
369 };
370
371 CHKERR get_options();
372
374}
375//! [Create common data]
376
377//! [Boundary condition]
380
381 auto pipeline_mng = mField.getInterface<PipelineManager>();
383 auto bc_mng = mField.getInterface<BcManager>();
384 auto time_scale = boost::make_shared<TimeScale>();
385
386 auto integration_rule = [](int, int, int approx_order) {
387 return 2 * (approx_order - 1);
388 };
389
390 using DomainNaturalBC =
392 using OpBodyForce =
394 using BoundaryNaturalBC =
397
398 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
400 pipeline_mng->getOpBoundaryRhsPipeline(), {NOSPACE}, "GEOMETRY");
401 //! [Natural boundary condition]
403 pipeline_mng->getOpBoundaryRhsPipeline(), mField, "U", {time_scale},
404 "FORCE", Sev::inform);
405
406 //! [Define gravity vector]
408 pipeline_mng->getOpDomainRhsPipeline(), mField, "U", {time_scale},
409 "BODY_FORCE", Sev::inform);
410
411 // Essential BC
412 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
413 "U", 0, 0);
414 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
415 "U", 1, 1);
416 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
417 "U", 2, 2);
418 CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
419 simple->getProblemName(), "U");
420
422}
423//! [Boundary condition]
424
425//! [Push operators to pip]
428 auto pip_mng = mField.getInterface<PipelineManager>();
429
430 auto integration_rule_vol = [](int, int, int approx_order) {
431 return 2 * approx_order + geom_order - 1;
432 };
433
434 auto add_domain_base_ops = [&](auto &pip) {
437 "GEOMETRY");
439 };
440
441 auto add_domain_ops_lhs = [&](auto &pip) {
443
444 // This assemble A-matrix
445 using OpMassPressure = FormsIntegrators<DomainEleOp>::Assembly<
447 // This assemble B-matrix (preconditioned)
448 using OpMassPressureStab = FormsIntegrators<DomainEleOpStab>::Assembly<
450 //! [Operators used for RHS incompressible elasticity]
451
452 //! [Operators used for incompressible elasticity]
453 using OpGradSymTensorGrad = FormsIntegrators<DomainEleOp>::Assembly<
457 //! [Operators used for incompressible elasticity]
458
459 auto mat_D_ptr = boost::make_shared<MatrixDouble>();
460 constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
461 mat_D_ptr->resize(size_symm * size_symm, 1);
462
467
469 auto t_mat = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
470 t_mat(i, j, k, l) = -2. * mu * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
471
472 pip.push_back(new OpMixScalarTimesDiv(
473 "P", "U",
474 [](const double, const double, const double) constexpr { return -1.; },
475 true, false));
476 pip.push_back(new OpGradSymTensorGrad("U", "U", mat_D_ptr));
477
478 auto get_lambda_reciprocal = [](const double, const double, const double) {
479 return 1. / lambda;
480 };
481 if (poisson_ratio < 0.5)
482 pip.push_back(new OpMassPressure("P", "P", get_lambda_reciprocal));
483
484 double eps_stab = 0;
485 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-eps_stab", &eps_stab,
486 PETSC_NULLPTR);
487 if (eps_stab > 0)
488 pip.push_back(new OpMassPressureStab(
489 "P", "P", [eps_stab](double, double, double) { return eps_stab; }));
490
492 };
493
494 auto add_domain_ops_rhs = [&](auto &pip) {
496
497 //! [Operators used for RHS incompressible elasticity]
498 using OpDomainGradTimesTensor = FormsIntegrators<DomainEleOp>::Assembly<
500
501 using OpDivDeltaUTimesP = FormsIntegrators<DomainEleOp>::Assembly<
503
504 using OpBaseTimesScalarValues = FormsIntegrators<DomainEleOp>::Assembly<
506
507 //! [Operators used for RHS incompressible elasticity]
508
509 auto pressure_ptr = boost::make_shared<VectorDouble>();
510 pip.push_back(new OpCalculateScalarFieldValues("P", pressure_ptr));
511
512 auto div_u_ptr = boost::make_shared<VectorDouble>();
513 pip.push_back(
515
516 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
518 "U", grad_u_ptr));
519
520 auto strain_ptr = boost::make_shared<MatrixDouble>();
521 pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(grad_u_ptr, strain_ptr));
522
523 auto get_four_mu = [](const double, const double, const double) {
524 return -2. * mu;
525 };
526 auto minus_one = [](const double, const double, const double) constexpr {
527 return -1.;
528 };
529
530 pip.push_back(new OpDomainGradTimesTensor("U", strain_ptr, get_four_mu));
531
532 pip.push_back(new OpDivDeltaUTimesP("U", pressure_ptr, minus_one));
533
534 pip.push_back(new OpBaseTimesScalarValues("P", div_u_ptr, minus_one));
535
536 auto get_lambda_reciprocal = [](const double, const double, const double) {
537 return 1. / lambda;
538 };
539 if (poisson_ratio < 0.5) {
540 pip.push_back(new OpBaseTimesScalarValues("P", pressure_ptr,
541 get_lambda_reciprocal));
542 }
543
545 };
546
547 CHKERR add_domain_base_ops(pip_mng->getOpDomainLhsPipeline());
548 CHKERR add_domain_base_ops(pip_mng->getOpDomainRhsPipeline());
549 CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
550 CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
551
552 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule_vol);
553 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule_vol);
554
556}
557//! [Push operators to pip]
558
559//! [Solve]
560struct SetUpSchur {
561 static boost::shared_ptr<SetUpSchur>
564
565protected:
566 SetUpSchur() = default;
567};
568
571
573 auto pip_mng = mField.getInterface<PipelineManager>();
574 auto is_manager = mField.getInterface<ISManager>();
575
576 auto set_section_monitor = [&](auto solver) {
578 SNES snes;
579 CHKERR TSGetSNES(solver, &snes);
580 PetscViewerAndFormat *vf;
581 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
582 PETSC_VIEWER_DEFAULT, &vf);
583 CHKERR SNESMonitorSet(
584 snes,
585 (MoFEMErrorCode (*)(SNES, PetscInt, PetscReal,
586 void *))SNESMonitorFields,
587 vf, (MoFEMErrorCode (*)(void **))PetscViewerAndFormatDestroy);
589 };
590
591 auto scatter_create = [&](auto D, auto coeff) {
593 CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
594 ROW, "U", coeff, coeff, is);
595 int loc_size;
596 CHKERR ISGetLocalSize(is, &loc_size);
597 Vec v;
598 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
599 VecScatter scatter;
600 CHKERR VecScatterCreate(D, is, v, PETSC_NULLPTR, &scatter);
601 return std::make_tuple(SmartPetscObj<Vec>(v),
603 };
604
605 auto create_post_process_elements = [&]() {
606 auto pp_fe = boost::make_shared<PostProcEle>(mField);
607
608 auto push_vol_ops = [this](auto &pip) {
611 "GEOMETRY");
612
614 };
615
616 auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
618
619 auto &pip = pp_fe->getOpPtrVector();
620
622
623 auto x_ptr = boost::make_shared<MatrixDouble>();
624 pip.push_back(
625 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
626 auto u_ptr = boost::make_shared<MatrixDouble>();
627 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
628
629 auto pressure_ptr = boost::make_shared<VectorDouble>();
630 pip.push_back(new OpCalculateScalarFieldValues("P", pressure_ptr));
631
632 auto div_u_ptr = boost::make_shared<VectorDouble>();
634 "U", div_u_ptr));
635
636 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
638 "U", grad_u_ptr));
639
640 auto strain_ptr = boost::make_shared<MatrixDouble>();
641 pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(grad_u_ptr, strain_ptr));
642
643 auto stress_ptr = boost::make_shared<MatrixDouble>();
644 pip.push_back(new OpCalculateLameStress<SPACE_DIM>(
645 mu, stress_ptr, strain_ptr, pressure_ptr));
646
647 pip.push_back(
648
649 new OpPPMap(
650
651 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
652
653 {{"P", pressure_ptr}},
654
655 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
656
657 {},
658
659 {{"STRAIN", strain_ptr}, {"STRESS", stress_ptr}}
660
661 )
662
663 );
664
666 };
667
668 auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
669 PetscBool post_proc_vol = PETSC_FALSE;
670 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_vol",
671 &post_proc_vol, PETSC_NULLPTR);
672 if (post_proc_vol == PETSC_FALSE)
673 return boost::shared_ptr<PostProcEle>();
674 auto pp_fe = boost::make_shared<PostProcEle>(mField);
676 push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
677 "push_vol_post_proc_ops");
678 return pp_fe;
679 };
680
681 auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
682 PetscBool post_proc_skin = PETSC_TRUE;
683 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_skin",
684 &post_proc_skin, PETSC_NULLPTR);
685 if (post_proc_skin == PETSC_FALSE)
686 return boost::shared_ptr<SkinPostProcEle>();
687
689 auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
690 auto op_side = new OpLoopSide<DomainEle>(
691 mField, simple->getDomainFEName(), SPACE_DIM, Sev::verbose);
692 pp_fe->getOpPtrVector().push_back(op_side);
693 CHK_MOAB_THROW(push_vol_post_proc_ops(
694 pp_fe, push_vol_ops(op_side->getOpPtrVector())),
695 "push_vol_post_proc_ops");
696 return pp_fe;
697 };
698
699 return std::make_pair(vol_post_proc(), skin_post_proc());
700 };
701
702 boost::shared_ptr<SetPtsData> field_eval_data;
703 boost::shared_ptr<MatrixDouble> vector_field_ptr;
704
705 if (doEvalField) {
706 vector_field_ptr = boost::make_shared<MatrixDouble>();
707 field_eval_data =
708 mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
709 CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree<SPACE_DIM>(
710 field_eval_data, simple->getDomainFEName());
711
712 field_eval_data->setEvalPoints(fieldEvalCoords.data(), 1);
713 auto no_rule = [](int, int, int) { return -1; };
714 auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
715 field_eval_fe_ptr->getRuleHook = no_rule;
716 field_eval_fe_ptr->getOpPtrVector().push_back(
717 new OpCalculateVectorFieldValues<SPACE_DIM>("U", vector_field_ptr));
718 }
719
720 auto set_time_monitor = [&](auto dm, auto solver) {
722 boost::shared_ptr<MonitorIncompressible> monitor_ptr(
724 create_post_process_elements(), uXScatter,
725 uYScatter, uZScatter, fieldEvalCoords,
726 field_eval_data, vector_field_ptr));
727 boost::shared_ptr<ForcesAndSourcesCore> null;
728 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
729 monitor_ptr, null, null);
731 };
732
733 auto set_essential_bc = [&]() {
735 // This is low level pushing finite elements (pipelines) to solver
736 auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
737 auto pre_proc_ptr = boost::make_shared<FEMethod>();
738 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
739
740 // Add boundary condition scaling
741 auto time_scale = boost::make_shared<TimeScale>();
742
743 pre_proc_ptr->preProcessHook = EssentialPreProc<DisplacementCubitBcData>(
744 mField, pre_proc_ptr, {time_scale}, false);
745 post_proc_rhs_ptr->postProcessHook =
746 EssentialPostProcRhs<DisplacementCubitBcData>(mField, post_proc_rhs_ptr,
747 1.);
748 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
749 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
750 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
752 };
753
754 auto set_schur_pc = [&](auto solver) {
755 SNES snes;
756 CHKERR TSGetSNES(solver, &snes);
757 KSP ksp;
758 CHKERR SNESGetKSP(snes, &ksp);
759 PC pc;
760 CHKERR KSPGetPC(ksp, &pc);
761 PetscBool is_pcfs = PETSC_FALSE;
762 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
763 boost::shared_ptr<SetUpSchur> schur_ptr;
764 auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
765
766 auto A = createDMMatrix(simple->getDM());
767 auto B = matDuplicate(A, MAT_DO_NOT_COPY_VALUES);
768
770 TSSetIJacobian(solver, A, B, TsSetIJacobian, ts_ctx_ptr.get()),
771 "set operators");
772 auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
773 auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
774 pre_proc_schur_lhs_ptr->preProcessHook = [pre_proc_schur_lhs_ptr]() {
776 MOFEM_LOG("INCOMP_ELASTICITY", Sev::verbose) << "Lhs Zero matrices";
777 CHKERR MatZeroEntries(pre_proc_schur_lhs_ptr->A);
778 CHKERR MatZeroEntries(pre_proc_schur_lhs_ptr->B);
780 };
781 post_proc_schur_lhs_ptr->postProcessHook = [this,
782 post_proc_schur_lhs_ptr]() {
784 MOFEM_LOG("INCOMP_ELASTICITY", Sev::verbose) << "Lhs Assemble Begin";
785 *(post_proc_schur_lhs_ptr->matAssembleSwitch) = false;
786 CHKERR MatAssemblyBegin(post_proc_schur_lhs_ptr->A, MAT_FINAL_ASSEMBLY);
787 CHKERR MatAssemblyEnd(post_proc_schur_lhs_ptr->A, MAT_FINAL_ASSEMBLY);
789 mField, post_proc_schur_lhs_ptr, 1.,
790 SmartPetscObj<Mat>(post_proc_schur_lhs_ptr->A, true))();
791
792 CHKERR MatAssemblyBegin(post_proc_schur_lhs_ptr->B, MAT_FINAL_ASSEMBLY);
793 CHKERR MatAssemblyEnd(post_proc_schur_lhs_ptr->B, MAT_FINAL_ASSEMBLY);
794 CHKERR MatAXPY(post_proc_schur_lhs_ptr->B, 1, post_proc_schur_lhs_ptr->A,
795 SAME_NONZERO_PATTERN);
796 MOFEM_LOG("INCOMP_ELASTICITY", Sev::verbose) << "Lhs Assemble End";
798 };
799 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
800 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
801
802 if (is_pcfs == PETSC_TRUE) {
803 if (AT == AssemblyType::SCHUR) {
804 schur_ptr = SetUpSchur::createSetUpSchur(mField);
805 CHK_MOAB_THROW(schur_ptr->setUp(solver), "setup schur preconditioner");
806 } else {
807 auto set_pcfieldsplit_preconditioned_ts = [&](auto solver) {
809 auto name_prb = simple->getProblemName();
811 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
812 name_prb, ROW, "P", 0, 1, is_p);
813 CHKERR PCFieldSplitSetIS(pc, PETSC_NULLPTR, is_p);
815 };
816 CHK_MOAB_THROW(set_pcfieldsplit_preconditioned_ts(solver),
817 "set pcfieldsplit preconditioned");
818 }
819 return boost::make_tuple(schur_ptr, A, B);
820 }
821
822 return boost::make_tuple(schur_ptr, A, B);
823 };
824
825 auto dm = simple->getDM();
826 auto D = createDMVector(dm);
827
828 uXScatter = scatter_create(D, 0);
829 uYScatter = scatter_create(D, 1);
830 if (SPACE_DIM == 3)
831 uZScatter = scatter_create(D, 2);
832
833 // Add extra finite elements to SNES solver pipelines to resolve essential
834 // boundary conditions
835 CHKERR set_essential_bc();
836
837 auto solver = pip_mng->createTSIM();
838 CHKERR TSSetFromOptions(solver);
839
840 CHKERR set_section_monitor(solver);
841 CHKERR set_time_monitor(dm, solver);
842 auto [schur_pc_ptr, A, B] = set_schur_pc(solver);
843
844 CHKERR TSSetSolution(solver, D);
845 CHKERR TSSetUp(solver);
846 CHKERR TSSolve(solver, NULL);
847
849}
850//! [Solve]
851
852//! [Check]
855
856 int atom_test = 0;
857 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-atom_test", &atom_test,
858 PETSC_NULLPTR);
859
860 if (atom_test) {
861 double eps = 1e-6;
862 double Ux_ref, Uy_ref, Uz_ref;
863 string test_name;
864
865 auto vector_field_ptr = boost::make_shared<MatrixDouble>();
866
867 auto field_eval_data =
868 mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
869 CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree<SPACE_DIM>(
870 field_eval_data, mField.getInterface<Simple>()->getDomainFEName());
871
872 field_eval_data->setEvalPoints(fieldEvalCoords.data(), 1);
873 auto no_rule = [](int, int, int) { return -1; };
874 auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
875 field_eval_fe_ptr->getRuleHook = no_rule;
876 field_eval_fe_ptr->getOpPtrVector().push_back(
877 new OpCalculateVectorFieldValues<SPACE_DIM>("U", vector_field_ptr));
878
880 ->evalFEAtThePoint<SPACE_DIM>(
881 fieldEvalCoords.data(), 1e-12,
883 mField.getInterface<Simple>()->getDomainFEName(), field_eval_data,
885 QUIET);
886
887 auto eval_num_vec = createVectorMPI(mField.get_comm(), PETSC_DECIDE, 1);
888 CHKERR VecZeroEntries(eval_num_vec);
889 if (vector_field_ptr->size1()) {
890 CHKERR VecSetValue(eval_num_vec, 0, 1, ADD_VALUES);
891 }
892 CHKERR VecAssemblyBegin(eval_num_vec);
893 CHKERR VecAssemblyEnd(eval_num_vec);
894
895 double eval_num;
896 CHKERR VecSum(eval_num_vec, &eval_num);
897 if (!(int)eval_num) { /// no elements to evaluate
898 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
899 "atom test %d failed: did not find elements to evaluate "
900 "the field, check the coordinates",
901 atom_test);
902 }
903
904 if (vector_field_ptr->size1()) {
905 auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vector_field_ptr);
906 switch (atom_test) {
907 case 1: // Cooke's Membrane
908 Ux_ref = -5.89757;
909 Uy_ref = 8.15644;
910 Uz_ref = 0.0;
911 test_name = "Cooke's Membrane";
912 break;
913 default:
914 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
915 "atom test %d does not exist", atom_test);
916 }
917
918 auto calculate_error = [](double computed, double reference,
919 double tolerance) {
920 return fabs(computed - reference) > tolerance * fabs(reference);
921 };
922
923 if (calculate_error(t_disp(0), Ux_ref, eps) ||
924 calculate_error(t_disp(1), Uy_ref, eps) ||
925 (SPACE_DIM == 3 && calculate_error(t_disp(2), Uz_ref, eps))) {
926 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
927 "atom test %d:%s failed: Displacement exceeds tolerance: Ux = "
928 "%f, Uy = %f",
929 atom_test, test_name.c_str(), t_disp(0), t_disp(1));
930 }
931 }
932 }
933
935}
936//! [Check]
937
938static char help[] = "...\n\n";
939
940int main(int argc, char *argv[]) {
941
942 // Initialisation of MoFEM/PETSc and MOAB data structures
943 const char param_file[] = "param_file.petsc";
944 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
945
946 // Add logging channel for CONTACT
947 auto core_log = logging::core::get();
948 core_log->add_sink(
949 LogManager::createSink(LogManager::getStrmWorld(), "INCOMP_ELASTICITY"));
950 LogManager::setLog("INCOMP_ELASTICITY");
951 MOFEM_LOG_TAG("INCOMP_ELASTICITY", "Indent");
952
953 try {
954
955 //! [Register MoFEM discrete manager in PETSc]
956 DMType dm_name = "DMMOFEM";
957 CHKERR DMRegister_MoFEM(dm_name);
958 //! [Register MoFEM discrete manager in PETSc
959
960 //! [Create MoAB]
961 moab::Core mb_instance; ///< mesh database
962 moab::Interface &moab = mb_instance; ///< mesh database interface
963 //! [Create MoAB]
964
965 //! [Create MoFEM]
966 MoFEM::Core core(moab); ///< finite element database
967 MoFEM::Interface &m_field = core; ///< finite element database interface
968 //! [Create MoFEM]
969
970 //! [Load mesh]
971 Simple *simple = m_field.getInterface<Simple>();
973 CHKERR simple->loadFile("");
974 //! [Load mesh]
975
976 //! [CONTACT]
977 Incompressible ex(m_field);
978 CHKERR ex.runProblem();
979 //! [CONTACT]
980 }
982
984
985 return 0;
986}
987
988struct SetUpSchurImpl : public SetUpSchur {
989
991 virtual ~SetUpSchurImpl() { S.reset(); }
993
994private:
997
999
1001
1003
1005};
1006
1009 auto simple = mField.getInterface<Simple>();
1010 auto pip = mField.getInterface<PipelineManager>();
1011 auto dm = simple->getDM();
1012
1013 SNES snes;
1014 CHKERR TSGetSNES(solver, &snes);
1015 KSP ksp;
1016 CHKERR SNESGetKSP(snes, &ksp);
1017 PC pc;
1018 CHKERR KSPGetPC(ksp, &pc);
1019
1020 MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform) << "Setup Schur pc";
1021
1022 if (S) {
1024 "Is expected that schur matrix is not allocated. This is "
1025 "possible only is PC is set up twice");
1026 }
1027
1028 auto create_sub_dm = [&]() {
1029 auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
1030 auto set_up = [&]() {
1032 CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "SUB");
1033 CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
1034 CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
1035 CHKERR DMMoFEMAddSubFieldRow(sub_dm, "U");
1036 CHKERR DMSetUp(sub_dm);
1038 };
1039 CHK_THROW_MESSAGE(set_up(), "sey up dm");
1040 return sub_dm;
1041 };
1042
1043 subDM = create_sub_dm();
1045 CHKERR MatSetBlockSize(S, SPACE_DIM);
1046
1047 auto dm_is = getDMSubData(subDM)->getSmartRowIs();
1048 auto ao_up = createAOMappingIS(dm_is, PETSC_NULLPTR);
1049 // Domain
1050 pip->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
1051 pip->getOpDomainLhsPipeline().push_back(
1052 createOpSchurAssembleEnd({"P"}, {nullptr}, ao_up, S, false, false));
1053
1054 auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1055 auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1056
1057 pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
1059 CHKERR MatZeroEntries(S);
1061 };
1062
1063 post_proc_schur_lhs_ptr->postProcessHook = [this, ao_up,
1064 post_proc_schur_lhs_ptr]() {
1066
1067 CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
1068 CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
1070 mField, post_proc_schur_lhs_ptr, 1, S, ao_up)();
1071
1072#ifndef NDEBUG
1073 auto print_mat_norm = [this](auto a, std::string prefix) {
1075 double nrm;
1076 CHKERR MatNorm(a, NORM_FROBENIUS, &nrm);
1077 MOFEM_LOG("INCOMP_ELASTICITY", Sev::verbose)
1078 << prefix << " norm = " << nrm;
1080 };
1081 CHKERR print_mat_norm(S, "S");
1082#endif // NDEBUG
1084 };
1085
1086 auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
1087 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
1088 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
1089
1090 SmartPetscObj<IS> is_p;
1091 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1092 simple->getProblemName(), ROW, "P", 0, 1, is_p);
1093 CHKERR PCFieldSplitSetIS(pc, NULL, is_p);
1094 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1095
1097}
1098
1099boost::shared_ptr<SetUpSchur>
1101 return boost::shared_ptr<SetUpSchur>(new SetUpSchurImpl(m_field));
1102}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
int main()
constexpr double a
static const double eps
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
[Define dimension]
Kronecker Delta class symmetric.
#define SCHUR_ASSEMBLE
Definition contact.cpp:18
@ QUIET
@ ROW
#define CATCH_ERRORS
Catch errors.
@ MF_EXIST
FieldApproximationBase
approximation base
Definition definitions.h:58
@ LASTBASE
Definition definitions.h:69
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
@ L2
field with C-1 continuity
Definition definitions.h:88
@ H1
continuous field
Definition definitions.h:85
@ NOSPACE
Definition definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_INVALID_DATA
Definition definitions.h:36
CoordinateTypes
Coodinate system.
@ CARTESIAN
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
auto integration_rule
constexpr auto t_kd
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, COORD_TYPE > OpMixScalarTimesDiv
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition DMMoFEM.cpp:215
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:488
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition DMMoFEM.cpp:450
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition DMMoFEM.cpp:238
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:576
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1102
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition DMMoFEM.cpp:410
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:1059
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
FTensor::Index< 'i', SPACE_DIM > i
static double young_modulus
constexpr AssemblyType AT
constexpr IntegrationType IT
static char help[]
[Check]
constexpr CoordinateTypes coord_type
static double lambda
constexpr int SPACE_DIM
static double poisson_ratio
PetscBool isDiscontinuousPressure
static double mu
int order
[Specialisation for assembly]
PetscBool doEvalField
double D
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition DMMoFEM.cpp:1046
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
Definition TsCtx.cpp:169
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition DMMoFEM.hpp:1144
OpSchurAssembleBase * createOpSchurAssembleEnd(std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range > > field_ents, SmartPetscObj< AO > ao, SmartPetscObj< Mat > schur, bool sym_schur, bool symm_op)
Construct a new Op Schur Assemble End object.
Definition Schur.cpp:2585
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
static FTensor::Ddg< FTensor::PackPtr< T *, 1 >, Tensor_Dim01, Tensor_Dim23 > getFTensor4DdgFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 4 on first two and last indices from form data matrix.
auto getDMSubData(DM dm)
Get sub problem data structure.
Definition DMMoFEM.hpp:1160
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
static auto getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
OpSchurAssembleBase * createOpSchurAssembleBegin()
Definition Schur.cpp:2580
constexpr AssemblyType A
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
int atom_test
Atom test.
Definition plastic.cpp:121
#define EXECUTABLE_DIMENSION
Definition plastic.cpp:13
constexpr auto size_symm
Definition plastic.cpp:42
static constexpr int approx_order
FTensor::Index< 'm', 3 > m
Element used to specialise assembly.
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
MoFEMErrorCode checkResults()
[Solve]
MoFEMErrorCode createCommonData()
[Set up problem]
MoFEMErrorCode setupProblem()
[Run problem]
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
MoFEMErrorCode runProblem()
[Run problem]
MoFEMErrorCode OPs()
[Boundary condition]
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
MoFEMErrorCode bC()
[Create common data]
std::array< double, 3 > fieldEvalCoords
Incompressible(MoFEM::Interface &m_field)
Add operators pushing bases from local to physical configuration.
boost::weak_ptr< CacheTuple > getCacheWeakPtr() const
Get the cache weak ptr object.
Simple interface for fast problem set-up.
Definition BcManager.hpp:29
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:118
Deprecated interface functions.
Definition of the displacement bc data structure.
Definition BCData.hpp:72
Data on single entity (This is passed as argument to DataOperator::doWork)
Class (Function) to enforce essential constrains on the left hand side diagonal.
Definition Essential.hpp:33
Class (Function) to enforce essential constrains on the right hand side diagonal.
Definition Essential.hpp:41
Class (Function) to enforce essential constrains.
Definition Essential.hpp:25
structure for User Loop Methods on finite elements
Field evaluator interface.
structure to get information form mofem into EntitiesFieldData
Section manager is used to create indexes and sections.
Definition ISManager.hpp:23
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Assembly methods.
Definition Natural.hpp:65
boost::function< MoFEMErrorCode( ForcesAndSourcesCore::UserDataOperator *op_ptr, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, MatrixDouble &m)> MatSetValuesHook
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Get values at integration pts for tensor field rank 1, i.e. vector field.
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
PipelineManager interface.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
const std::string getProblemName() const
Get the Problem Name.
Definition Simple.hpp:410
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition Simple.hpp:389
intrusive_ptr for managing petsc objects
PetscReal ts_t
time
PetscInt ts_step
time step number
Vec & ts_u
state vector
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
MoFEMErrorCode preProcess()
function is run at the beginning of loop
boost::shared_ptr< MatrixDouble > vecFieldPtr
std::array< double, 3 > fieldEvalCoords
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
boost::shared_ptr< SkinPostProcEle > skinPostProcFe
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
MonitorIncompressible(SmartPetscObj< DM > dm, std::pair< boost::shared_ptr< PostProcEle >, boost::shared_ptr< SkinPostProcEle > > pair_post_proc_fe, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > ux_scatter, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uy_scatter, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uz_scatter, std::array< double, 3 > pass_field_eval_coords, boost::shared_ptr< SetPtsData > pass_field_eval_data, boost::shared_ptr< MatrixDouble > vec_field_ptr)
MoFEMErrorCode operator()()
function is run for every finite element
boost::shared_ptr< SetPtsData > fieldEvalData
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
boost::shared_ptr< PostProcEle > postProcFe
MoFEMErrorCode postProcess()
function is run at the end of loop
OpCalculateLameStress(double m_u, boost::shared_ptr< MatrixDouble > stress_ptr, boost::shared_ptr< MatrixDouble > strain_ptr, boost::shared_ptr< VectorDouble > pressure_ptr)
boost::shared_ptr< MatrixDouble > stressPtr
boost::shared_ptr< VectorDouble > pressurePtr
boost::shared_ptr< MatrixDouble > strainPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
SmartPetscObj< DM > subDM
field split sub dm
Definition plastic.cpp:1680
SmartPetscObj< Mat > S
SetUpSchurImpl(MoFEM::Interface &m_field)
SmartPetscObj< DM > createSubDM()
MoFEMErrorCode setUp(SmartPetscObj< TS > solver)
MoFEMErrorCode setUp(SmartPetscObj< KSP >)
MoFEMErrorCode setPC(PC pc)
MoFEMErrorCode setOperator()
MoFEM::Interface & mField
[Push operators to pipeline]
SetUpSchur()=default
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
virtual MoFEMErrorCode setUp(SmartPetscObj< TS > solver)=0
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)