v0.14.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
30
31// struct DomainBCs {};
32// struct BoundaryBCs {};
33
34// using DomainRhsBCs = NaturalBC<DomainEleOp>::Assembly<AT>::LinearForm<IT>;
35// using OpDomainRhsBCs = DomainRhsBCs::OpFlux<DomainBCs, 1, SPACE_DIM>;
36// using BoundaryRhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::LinearForm<IT>;
37// using OpBoundaryRhsBCs = BoundaryRhsBCs::OpFlux<BoundaryBCs, 1, SPACE_DIM>;
38// using BoundaryLhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::BiLinearForm<IT>;
39// using OpBoundaryLhsBCs = BoundaryLhsBCs::OpFlux<BoundaryBCs, 1, SPACE_DIM>;
40
41PetscBool doEvalField;
43
46 std::pair<boost::shared_ptr<PostProcEle>,
47 boost::shared_ptr<SkinPostProcEle>>
48 pair_post_proc_fe,
49 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> ux_scatter,
50 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uy_scatter,
51 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uz_scatter,
52 std::array<double, SPACE_DIM> pass_field_eval_coords,
53 boost::shared_ptr<SetPtsData> pass_field_eval_data,
54 boost::shared_ptr<MatrixDouble> vec_field_ptr)
55 : dM(dm), uXScatter(ux_scatter), uYScatter(uy_scatter),
56 uZScatter(uz_scatter), fieldEvalCoords(pass_field_eval_coords),
57 fieldEvalData(pass_field_eval_data), vecFieldPtr(vec_field_ptr) {
58 postProcFe = pair_post_proc_fe.first;
59 skinPostProcFe = pair_post_proc_fe.second;
60 };
61
62 MoFEMErrorCode preProcess() { return 0; }
63 MoFEMErrorCode operator()() { return 0; }
64
67
68 MoFEM::Interface *m_field_ptr;
69 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
70 auto *simple = m_field_ptr->getInterface<Simple>();
71
72 if (doEvalField) {
73 if (SPACE_DIM == 3) {
75 ->evalFEAtThePoint3D(
76 fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
77 simple->getDomainFEName(), fieldEvalData,
78 m_field_ptr->get_comm_rank(), m_field_ptr->get_comm_rank(),
79 getCacheWeakPtr().lock(), MF_EXIST, QUIET);
80 } else {
82 ->evalFEAtThePoint2D(
83 fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
84 simple->getDomainFEName(), fieldEvalData,
85 m_field_ptr->get_comm_rank(), m_field_ptr->get_comm_rank(),
86 getCacheWeakPtr().lock(), MF_EXIST, QUIET);
87 }
88
89 if (vecFieldPtr->size1()) {
90 auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vecFieldPtr);
91 if constexpr (SPACE_DIM == 2)
92 MOFEM_LOG("SYNC", Sev::inform)
93 << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1);
94 else
95 MOFEM_LOG("SYNC", Sev::inform)
96 << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1)
97 << " U_Z: " << t_disp(2);
98 }
99 }
100 MOFEM_LOG_SEVERITY_SYNC(m_field_ptr->get_comm(), Sev::inform);
101
102 auto make_vtk = [&]() {
104 if (postProcFe) {
107 CHKERR postProcFe->writeFile("out_incomp_elasticity" +
108 boost::lexical_cast<std::string>(ts_step) +
109 ".h5m");
110 }
111 if (skinPostProcFe) {
114 CHKERR skinPostProcFe->writeFile(
115 "out_skin_incomp_elasticity_" +
116 boost::lexical_cast<std::string>(ts_step) + ".h5m");
117 }
119 };
120
121 auto print_max_min = [&](auto &tuple, const std::string msg) {
123 CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
124 INSERT_VALUES, SCATTER_FORWARD);
125 CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
126 INSERT_VALUES, SCATTER_FORWARD);
127 double max, min;
128 CHKERR VecMax(std::get<0>(tuple), PETSC_NULL, &max);
129 CHKERR VecMin(std::get<0>(tuple), PETSC_NULL, &min);
130 MOFEM_LOG_C("INCOMP_ELASTICITY", Sev::inform,
131 "%s time %3.4e min %3.4e max %3.4e", msg.c_str(), ts_t, min,
132 max);
134 };
135
136 CHKERR make_vtk();
137 CHKERR print_max_min(uXScatter, "Ux");
138 CHKERR print_max_min(uYScatter, "Uy");
139 if constexpr (SPACE_DIM == 3)
140 CHKERR print_max_min(uZScatter, "Uz");
141
143 }
144
145private:
147 boost::shared_ptr<PostProcEle> postProcFe;
148 boost::shared_ptr<SkinPostProcEle> skinPostProcFe;
149 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
150 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
151 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
152 std::array<double, SPACE_DIM> fieldEvalCoords;
153 boost::shared_ptr<SetPtsData> fieldEvalData;
154 boost::shared_ptr<MatrixDouble> vecFieldPtr;
155};
156
157// Assemble to A matrix, by default, however, some terms are assembled only to
158// preconditioning.
159
160template <>
164 const EntitiesFieldData::EntData &row_data,
165 const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
166 return MatSetValues<AssemblyTypeSelector<AT>>(
167 op_ptr->getKSPA(), row_data, col_data, m, ADD_VALUES);
168 };
169
170/**
171 * @brief Element used to specialise assembly
172 *
173 */
175 using DomainEleOp::DomainEleOp;
176};
177
178/**
179 * @brief Specialise assembly for Stabilised matrix
180 *
181 * @tparam
182 */
183template <>
187 const EntitiesFieldData::EntData &row_data,
188 const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
189 return MatSetValues<AssemblyTypeSelector<AT>>(
190 op_ptr->getKSPB(), row_data, col_data, m, ADD_VALUES);
191 };
192//! [Specialisation for assembly]
193
194int order = 2;
196inline static double young_modulus = 100;
197inline static double poisson_ratio = 0.25;
198inline static double mu;
199inline static double lambda;
200
201PetscBool isDiscontinuousPressure = PETSC_FALSE;
202
204
205 Incompressible(MoFEM::Interface &m_field) : mField(m_field) {}
206
208
209private:
211
218
219 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
220 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
221 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
222};
223
224template <int DIM>
226 OpCalculateLameStress(double m_u, boost::shared_ptr<MatrixDouble> stress_ptr,
227 boost::shared_ptr<MatrixDouble> strain_ptr,
228 boost::shared_ptr<VectorDouble> pressure_ptr)
230 stressPtr(stress_ptr), strainPtr(strain_ptr),
231 pressurePtr(pressure_ptr) {}
232
236 // Define Indicies
239
240 // Define Kronecker Delta
242
243 // Number of Gauss points
244 const size_t nb_gauss_pts = getGaussPts().size2();
245
246 stressPtr->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts);
247 auto t_stress = getFTensor2SymmetricFromMat<DIM>(*(stressPtr));
248 auto t_strain = getFTensor2SymmetricFromMat<DIM>(*(strainPtr));
249 auto t_pressure = getFTensor0FromVec(*(pressurePtr));
250
251 const double l_mu = mU;
252 // Extract matrix from data matrix
253 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
254
255 t_stress(i, j) = t_pressure * t_kd(i, j) + 2. * l_mu * t_strain(i, j);
256
257 ++t_strain;
258 ++t_stress;
259 ++t_pressure;
260 }
261
263 }
264
265private:
266 double mU;
267 boost::shared_ptr<MatrixDouble> stressPtr;
268 boost::shared_ptr<MatrixDouble> strainPtr;
269 boost::shared_ptr<VectorDouble> pressurePtr;
270};
271
272//! [Run problem]
277 CHKERR bC();
278 CHKERR OPs();
279 CHKERR tsSolve();
282}
283//! [Run problem]
284
285//! [Set up problem]
289
290 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
291 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-geom_order", &geom_order,
292 PETSC_NULL);
293 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_discontinuous_pressure",
294 &isDiscontinuousPressure, PETSC_NULL);
295
296 MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform) << "Order " << order;
297 MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform) << "Geom order " << geom_order;
298
299 // Select base
300 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
301 const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
302 PetscInt choice_base_value = AINSWORTH;
303 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
304 LASBASETOPT, &choice_base_value, PETSC_NULL);
305
307 switch (choice_base_value) {
308 case AINSWORTH:
310 MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
311 << "Set AINSWORTH_LEGENDRE_BASE for displacements";
312 break;
313 case DEMKOWICZ:
315 MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
316 << "Set DEMKOWICZ_JACOBI_BASE for displacements";
317 break;
318 default:
319 base = LASTBASE;
320 break;
321 }
322
323 // Note: For tets we have only H1 Ainsworth base, for Hex we have only H1
324 // Demkowicz base. We need to implement Demkowicz H1 base on tet.
325 CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
326 CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
327
328 // Adding fields related to incompressible elasticity
329 // Add displacement domain and boundary fields
330 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
331 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
332 CHKERR simple->setFieldOrder("U", order);
333
334 // Add pressure domain and boundary fields
335 // Choose either Crouzeix-Raviart element:
337 CHKERR simple->addDomainField("P", L2, base, 1);
338 CHKERR simple->setFieldOrder("P", order - 2);
339 } else {
340 // ... or Taylor-Hood element:
341 CHKERR simple->addDomainField("P", H1, base, 1);
342 CHKERR simple->setFieldOrder("P", order - 1);
343 }
344
345 // Add geometry data field
346 CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
347 CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
348
349 CHKERR simple->setUp();
350
351 auto project_ho_geometry = [&]() {
352 Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
353 return mField.loop_dofs("GEOMETRY", ent_method);
354 };
355 CHKERR project_ho_geometry();
356
358} //! [Set up problem]
359
360//! [Create common data]
363
364 auto get_options = [&]() {
366 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
367 &young_modulus, PETSC_NULL);
368 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
369 &poisson_ratio, PETSC_NULL);
370
371 MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
372 << "Young modulus " << young_modulus;
373 MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform)
374 << "Poisson_ratio " << poisson_ratio;
375
376 mu = young_modulus / (2. + 2. * poisson_ratio);
377 const double lambda_denom =
378 (1. + poisson_ratio) * (1. - 2. * poisson_ratio);
379 lambda = young_modulus * poisson_ratio / lambda_denom;
380
382 };
383
384 CHKERR get_options();
385
387}
388//! [Create common data]
389
390//! [Boundary condition]
393
394 auto pipeline_mng = mField.getInterface<PipelineManager>();
396 auto bc_mng = mField.getInterface<BcManager>();
397 auto time_scale = boost::make_shared<TimeScale>();
398
399 auto integration_rule = [](int, int, int approx_order) {
400 return 2 * (approx_order - 1);
401 };
402
403 using DomainNaturalBC =
405 using OpBodyForce =
407 using BoundaryNaturalBC =
410
411 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
413 pipeline_mng->getOpBoundaryRhsPipeline(), {NOSPACE}, "GEOMETRY");
414 //! [Natural boundary condition]
416 pipeline_mng->getOpBoundaryRhsPipeline(), mField, "U", {time_scale},
417 "FORCE", Sev::inform);
418
419 //! [Define gravity vector]
421 pipeline_mng->getOpDomainRhsPipeline(), mField, "U", {time_scale},
422 "BODY_FORCE", Sev::inform);
423
424 // Essential BC
425 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
426 "U", 0, 0);
427 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
428 "U", 1, 1);
429 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
430 "U", 2, 2);
431 CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
432 simple->getProblemName(), "U");
433
435}
436//! [Boundary condition]
437
438//! [Push operators to pip]
442 auto pip_mng = mField.getInterface<PipelineManager>();
443 auto bc_mng = mField.getInterface<BcManager>();
444
445 auto integration_rule_vol = [](int, int, int approx_order) {
446 return 2 * approx_order + geom_order - 1;
447 };
448
449 auto add_domain_base_ops = [&](auto &pip) {
452 "GEOMETRY");
454 };
455
456 auto add_domain_ops_lhs = [&](auto &pip) {
458
459 // This assemble A-matrix
460 using OpMassPressure = FormsIntegrators<DomainEleOp>::Assembly<
462 // This assemble B-matrix (preconditioned)
463 using OpMassPressureStab = FormsIntegrators<DomainEleOpStab>::Assembly<
465 //! [Operators used for RHS incompressible elasticity]
466
467 //! [Operators used for incompressible elasticity]
468 using OpGradSymTensorGrad = FormsIntegrators<DomainEleOp>::Assembly<
472 //! [Operators used for incompressible elasticity]
473
474 auto mat_D_ptr = boost::make_shared<MatrixDouble>();
475 constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
476 mat_D_ptr->resize(size_symm * size_symm, 1);
477
482
484 auto t_mat = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
485 t_mat(i, j, k, l) = -2. * mu * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
486
487 pip.push_back(new OpMixScalarTimesDiv(
488 "P", "U",
489 [](const double, const double, const double) constexpr { return -1.; },
490 true, false));
491 pip.push_back(new OpGradSymTensorGrad("U", "U", mat_D_ptr));
492
493 auto get_lambda_reciprocal = [](const double, const double, const double) {
494 return 1. / lambda;
495 };
496 if (poisson_ratio < 0.5)
497 pip.push_back(new OpMassPressure("P", "P", get_lambda_reciprocal));
498
499 double eps_stab = 0;
500 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,
501 PETSC_NULL);
502 if (eps_stab > 0)
503 pip.push_back(new OpMassPressureStab(
504 "P", "P", [eps_stab](double, double, double) { return eps_stab; }));
505
507 };
508
509 auto add_domain_ops_rhs = [&](auto &pip) {
511
512 //! [Operators used for RHS incompressible elasticity]
513 using OpDomainGradTimesTensor = FormsIntegrators<DomainEleOp>::Assembly<
515
516 using OpDivDeltaUTimesP = FormsIntegrators<DomainEleOp>::Assembly<
518
519 using OpBaseTimesScalarValues = FormsIntegrators<DomainEleOp>::Assembly<
521
522 //! [Operators used for RHS incompressible elasticity]
523
524 auto pressure_ptr = boost::make_shared<VectorDouble>();
525 pip.push_back(new OpCalculateScalarFieldValues("P", pressure_ptr));
526
527 auto div_u_ptr = boost::make_shared<VectorDouble>();
528 pip.push_back(
530
531 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
533 "U", grad_u_ptr));
534
535 auto strain_ptr = boost::make_shared<MatrixDouble>();
536 pip.push_back(
537 new OpSymmetrizeTensor<SPACE_DIM>("U", grad_u_ptr, strain_ptr));
538
539 auto get_four_mu = [](const double, const double, const double) {
540 return -2. * mu;
541 };
542 auto minus_one = [](const double, const double, const double) constexpr {
543 return -1.;
544 };
545
546 pip.push_back(new OpDomainGradTimesTensor("U", strain_ptr, get_four_mu));
547
548 pip.push_back(new OpDivDeltaUTimesP("U", pressure_ptr, minus_one));
549
550 pip.push_back(new OpBaseTimesScalarValues("P", div_u_ptr, minus_one));
551
552 auto get_lambda_reciprocal = [](const double, const double, const double) {
553 return 1. / lambda;
554 };
555 if (poisson_ratio < 0.5) {
556 pip.push_back(new OpBaseTimesScalarValues("P", pressure_ptr,
557 get_lambda_reciprocal));
558 }
559
561 };
562
563 CHKERR add_domain_base_ops(pip_mng->getOpDomainLhsPipeline());
564 CHKERR add_domain_base_ops(pip_mng->getOpDomainRhsPipeline());
565 CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
566 CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
567
568 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule_vol);
569 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule_vol);
570
572}
573//! [Push operators to pip]
574
575//! [Solve]
576struct SetUpSchur {
577 static boost::shared_ptr<SetUpSchur>
580
581protected:
582 SetUpSchur() = default;
583};
584
587
589 auto pip_mng = mField.getInterface<PipelineManager>();
590 auto is_manager = mField.getInterface<ISManager>();
591
592 auto set_section_monitor = [&](auto solver) {
594 SNES snes;
595 CHKERR TSGetSNES(solver, &snes);
596 PetscViewerAndFormat *vf;
597 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
598 PETSC_VIEWER_DEFAULT, &vf);
599 CHKERR SNESMonitorSet(
600 snes,
601 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
602 vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
604 };
605
606 auto scatter_create = [&](auto D, auto coeff) {
608 CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
609 ROW, "U", coeff, coeff, is);
610 int loc_size;
611 CHKERR ISGetLocalSize(is, &loc_size);
612 Vec v;
613 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
614 VecScatter scatter;
615 CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
616 return std::make_tuple(SmartPetscObj<Vec>(v),
618 };
619
620 auto create_post_process_elements = [&]() {
621 auto pp_fe = boost::make_shared<PostProcEle>(mField);
622 auto &pip = pp_fe->getOpPtrVector();
623
624 auto push_vol_ops = [this](auto &pip) {
627 "GEOMETRY");
628
630 };
631
632 auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
634
635 auto &pip = pp_fe->getOpPtrVector();
636
638
639 auto x_ptr = boost::make_shared<MatrixDouble>();
640 pip.push_back(
641 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
642 auto u_ptr = boost::make_shared<MatrixDouble>();
643 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
644
645 auto pressure_ptr = boost::make_shared<VectorDouble>();
646 pip.push_back(new OpCalculateScalarFieldValues("P", pressure_ptr));
647
648 auto div_u_ptr = boost::make_shared<VectorDouble>();
650 "U", div_u_ptr));
651
652 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
654 "U", grad_u_ptr));
655
656 auto strain_ptr = boost::make_shared<MatrixDouble>();
657 pip.push_back(
658 new OpSymmetrizeTensor<SPACE_DIM>("U", grad_u_ptr, strain_ptr));
659
660 auto stress_ptr = boost::make_shared<MatrixDouble>();
661 pip.push_back(new OpCalculateLameStress<SPACE_DIM>(
662 mu, stress_ptr, strain_ptr, pressure_ptr));
663
664 pip.push_back(
665
666 new OpPPMap(
667
668 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
669
670 {{"P", pressure_ptr}},
671
672 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
673
674 {},
675
676 {{"STRAIN", strain_ptr}, {"STRESS", stress_ptr}}
677
678 )
679
680 );
681
683 };
684
685 auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
686 PetscBool post_proc_vol = PETSC_FALSE;
687 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_vol",
688 &post_proc_vol, PETSC_NULL);
689 if (post_proc_vol == PETSC_FALSE)
690 return boost::shared_ptr<PostProcEle>();
691 auto pp_fe = boost::make_shared<PostProcEle>(mField);
693 push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
694 "push_vol_post_proc_ops");
695 return pp_fe;
696 };
697
698 auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
699 PetscBool post_proc_skin = PETSC_TRUE;
700 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_skin",
701 &post_proc_skin, PETSC_NULL);
702 if (post_proc_skin == PETSC_FALSE)
703 return boost::shared_ptr<SkinPostProcEle>();
704
706 auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
707 auto op_side = new OpLoopSide<DomainEle>(
708 mField, simple->getDomainFEName(), SPACE_DIM, Sev::verbose);
709 pp_fe->getOpPtrVector().push_back(op_side);
710 CHK_MOAB_THROW(push_vol_post_proc_ops(
711 pp_fe, push_vol_ops(op_side->getOpPtrVector())),
712 "push_vol_post_proc_ops");
713 return pp_fe;
714 };
715
716 return std::make_pair(vol_post_proc(), skin_post_proc());
717 };
718
719 boost::shared_ptr<SetPtsData> field_eval_data;
720 boost::shared_ptr<MatrixDouble> vector_field_ptr;
721
722 std::array<double, SPACE_DIM> field_eval_coords;
723 int coords_dim = SPACE_DIM;
724 CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
725 field_eval_coords.data(), &coords_dim,
726 &doEvalField);
727
728 if (doEvalField) {
729 vector_field_ptr = boost::make_shared<MatrixDouble>();
730 field_eval_data =
731 mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
732 if constexpr (SPACE_DIM == 3) {
733 CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree3D(
734 field_eval_data, simple->getDomainFEName());
735 } else {
736 CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree2D(
737 field_eval_data, simple->getDomainFEName());
738 }
739
740 field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
741 auto no_rule = [](int, int, int) { return -1; };
742 auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
743 field_eval_fe_ptr->getRuleHook = no_rule;
744 field_eval_fe_ptr->getOpPtrVector().push_back(
745 new OpCalculateVectorFieldValues<SPACE_DIM>("U", vector_field_ptr));
746 }
747
748 auto set_time_monitor = [&](auto dm, auto solver) {
750 boost::shared_ptr<MonitorIncompressible> monitor_ptr(
752 create_post_process_elements(), uXScatter,
753 uYScatter, uZScatter, field_eval_coords,
754 field_eval_data, vector_field_ptr));
755 boost::shared_ptr<ForcesAndSourcesCore> null;
756 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
757 monitor_ptr, null, null);
759 };
760
761 auto set_essential_bc = [&]() {
763 // This is low level pushing finite elements (pipelines) to solver
764 auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
765 auto pre_proc_ptr = boost::make_shared<FEMethod>();
766 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
767
768 // Add boundary condition scaling
769 auto time_scale = boost::make_shared<TimeScale>();
770
771 pre_proc_ptr->preProcessHook = EssentialPreProc<DisplacementCubitBcData>(
772 mField, pre_proc_ptr, {time_scale}, false);
773 post_proc_rhs_ptr->postProcessHook =
774 EssentialPostProcRhs<DisplacementCubitBcData>(mField, post_proc_rhs_ptr,
775 1.);
776 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
777 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
778 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
780 };
781
782 auto set_schur_pc = [&](auto solver) {
783 SNES snes;
784 CHKERR TSGetSNES(solver, &snes);
785 KSP ksp;
786 CHKERR SNESGetKSP(snes, &ksp);
787 PC pc;
788 CHKERR KSPGetPC(ksp, &pc);
789 PetscBool is_pcfs = PETSC_FALSE;
790 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
791 boost::shared_ptr<SetUpSchur> schur_ptr;
792 auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
793
794 auto A = createDMMatrix(simple->getDM());
795 auto B = matDuplicate(A, MAT_DO_NOT_COPY_VALUES);
796
798 TSSetIJacobian(solver, A, B, TsSetIJacobian, ts_ctx_ptr.get()),
799 "set operators");
800 auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
801 auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
802 pre_proc_schur_lhs_ptr->preProcessHook = [pre_proc_schur_lhs_ptr]() {
804 MOFEM_LOG("INCOMP_ELASTICITY", Sev::verbose) << "Lhs Zero matrices";
805 CHKERR MatZeroEntries(pre_proc_schur_lhs_ptr->A);
806 CHKERR MatZeroEntries(pre_proc_schur_lhs_ptr->B);
808 };
809 post_proc_schur_lhs_ptr->postProcessHook = [this,
810 post_proc_schur_lhs_ptr]() {
812 MOFEM_LOG("INCOMP_ELASTICITY", Sev::verbose) << "Lhs Assemble Begin";
813 *(post_proc_schur_lhs_ptr->matAssembleSwitch) = false;
814 CHKERR MatAssemblyBegin(post_proc_schur_lhs_ptr->A, MAT_FINAL_ASSEMBLY);
815 CHKERR MatAssemblyEnd(post_proc_schur_lhs_ptr->A, MAT_FINAL_ASSEMBLY);
817 mField, post_proc_schur_lhs_ptr, 1.,
818 SmartPetscObj<Mat>(post_proc_schur_lhs_ptr->A, true))();
819
820 CHKERR MatAssemblyBegin(post_proc_schur_lhs_ptr->B, MAT_FINAL_ASSEMBLY);
821 CHKERR MatAssemblyEnd(post_proc_schur_lhs_ptr->B, MAT_FINAL_ASSEMBLY);
822 CHKERR MatAXPY(post_proc_schur_lhs_ptr->B, 1, post_proc_schur_lhs_ptr->A,
823 SAME_NONZERO_PATTERN);
824 MOFEM_LOG("INCOMP_ELASTICITY", Sev::verbose) << "Lhs Assemble End";
826 };
827 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
828 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
829
830 if (is_pcfs == PETSC_TRUE) {
831 if (AT == AssemblyType::SCHUR) {
832 schur_ptr = SetUpSchur::createSetUpSchur(mField);
833 CHK_MOAB_THROW(schur_ptr->setUp(solver), "setup schur preconditioner");
834 } else {
835 auto set_pcfieldsplit_preconditioned_ts = [&](auto solver) {
837 auto bc_mng = mField.getInterface<BcManager>();
838 auto name_prb = simple->getProblemName();
840 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
841 name_prb, ROW, "P", 0, 1, is_p);
842 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_p);
844 };
845 CHK_MOAB_THROW(set_pcfieldsplit_preconditioned_ts(solver),
846 "set pcfieldsplit preconditioned");
847 }
848 return boost::make_tuple(schur_ptr, A, B);
849 }
850
851 return boost::make_tuple(schur_ptr, A, B);
852 };
853
854 auto dm = simple->getDM();
855 auto D = createDMVector(dm);
856
857 uXScatter = scatter_create(D, 0);
858 uYScatter = scatter_create(D, 1);
859 if (SPACE_DIM == 3)
860 uZScatter = scatter_create(D, 2);
861
862 // Add extra finite elements to SNES solver pipelines to resolve essential
863 // boundary conditions
864 CHKERR set_essential_bc();
865
866 auto solver = pip_mng->createTSIM();
867 CHKERR TSSetFromOptions(solver);
868
869 CHKERR set_section_monitor(solver);
870 CHKERR set_time_monitor(dm, solver);
871 auto [schur_pc_ptr, A, B] = set_schur_pc(solver);
872
873 CHKERR TSSetSolution(solver, D);
874 CHKERR TSSetUp(solver);
875 CHKERR TSSolve(solver, NULL);
876
878}
879//! [Solve]
880
881//! [Check]
883//! [Check]
884
885static char help[] = "...\n\n";
886
887int main(int argc, char *argv[]) {
888
889 // Initialisation of MoFEM/PETSc and MOAB data structures
890 const char param_file[] = "param_file.petsc";
892
893 // Add logging channel for CONTACT
894 auto core_log = logging::core::get();
895 core_log->add_sink(
896 LogManager::createSink(LogManager::getStrmWorld(), "INCOMP_ELASTICITY"));
897 LogManager::setLog("INCOMP_ELASTICITY");
898 MOFEM_LOG_TAG("INCOMP_ELASTICITY", "Indent");
899
900 try {
901
902 //! [Register MoFEM discrete manager in PETSc]
903 DMType dm_name = "DMMOFEM";
904 CHKERR DMRegister_MoFEM(dm_name);
905 //! [Register MoFEM discrete manager in PETSc
906
907 //! [Create MoAB]
908 moab::Core mb_instance; ///< mesh database
909 moab::Interface &moab = mb_instance; ///< mesh database interface
910 //! [Create MoAB]
911
912 //! [Create MoFEM]
913 MoFEM::Core core(moab); ///< finite element database
914 MoFEM::Interface &m_field = core; ///< finite element database interface
915 //! [Create MoFEM]
916
917 //! [Load mesh]
918 Simple *simple = m_field.getInterface<Simple>();
919 CHKERR simple->getOptions();
920 CHKERR simple->loadFile("");
921 //! [Load mesh]
922
923 //! [CONTACT]
924 Incompressible ex(m_field);
925 CHKERR ex.runProblem();
926 //! [CONTACT]
927 }
929
931
932 return 0;
933}
934
935struct SetUpSchurImpl : public SetUpSchur {
936
938 virtual ~SetUpSchurImpl() { S.reset(); }
940
941private:
944
946
948
950
952};
953
957 auto pip = mField.getInterface<PipelineManager>();
958 auto dm = simple->getDM();
959
960 SNES snes;
961 CHKERR TSGetSNES(solver, &snes);
962 KSP ksp;
963 CHKERR SNESGetKSP(snes, &ksp);
964 PC pc;
965 CHKERR KSPGetPC(ksp, &pc);
966
967 MOFEM_LOG("INCOMP_ELASTICITY", Sev::inform) << "Setup Schur pc";
968
969 if (S) {
971 "Is expected that schur matrix is not allocated. This is "
972 "possible only is PC is set up twice");
973 }
974
975 auto create_sub_dm = [&]() {
976 auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
977 auto set_up = [&]() {
979 CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "SUB");
980 CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
981 CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
982 CHKERR DMMoFEMAddSubFieldRow(sub_dm, "U");
983 CHKERR DMSetUp(sub_dm);
985 };
986 CHK_THROW_MESSAGE(set_up(), "sey up dm");
987 return sub_dm;
988 };
989
990 subDM = create_sub_dm();
992 CHKERR MatSetBlockSize(S, SPACE_DIM);
993
994 auto dm_is = getDMSubData(subDM)->getSmartRowIs();
995 auto ao_up = createAOMappingIS(dm_is, PETSC_NULL);
996 // Domain
997 pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
998 pip->getOpDomainLhsPipeline().push_back(new OpSchurAssembleEnd<SCHUR_DSYSV>(
999 {"P"}, {nullptr}, {ao_up}, {S}, {false}, false));
1000
1001 auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1002 auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1003
1004 pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
1006 CHKERR MatZeroEntries(S);
1008 };
1009
1010 post_proc_schur_lhs_ptr->postProcessHook = [this, ao_up,
1011 post_proc_schur_lhs_ptr]() {
1013
1014 auto print_mat_norm = [this](auto a, std::string prefix) {
1016 double nrm;
1017 CHKERR MatNorm(a, NORM_FROBENIUS, &nrm);
1018 MOFEM_LOG("INCOMP_ELASTICITY", Sev::verbose)
1019 << prefix << " norm = " << nrm;
1021 };
1022
1023 CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
1024 CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
1026 mField, post_proc_schur_lhs_ptr, 1, S, ao_up)();
1027
1028#ifndef NDEBUG
1029 CHKERR print_mat_norm(S, "S");
1030#endif // NDEBUG
1032 };
1033
1034 auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
1035 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
1036 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
1037
1038 SmartPetscObj<IS> is_p;
1039 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1040 simple->getProblemName(), ROW, "P", 0, 1, is_p);
1041 CHKERR PCFieldSplitSetIS(pc, NULL, is_p);
1042 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1043
1045}
1046
1047boost::shared_ptr<SetUpSchur>
1049 return boost::shared_ptr<SetUpSchur>(new SetUpSchurImpl(m_field));
1050}
static Index< 'p', 3 > p
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
Definition: LogManager.hpp:352
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
int main()
Definition: adol-c_atom.cpp:46
constexpr double a
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
@ QUIET
Definition: definitions.h:208
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ MF_EXIST
Definition: definitions.h:100
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ LASTBASE
Definition: definitions.h:69
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
@ 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 ...
Definition: definitions.h:346
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
CoordinateTypes
Coodinate system.
Definition: definitions.h:114
@ CARTESIAN
Definition: definitions.h:115
#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 SPACE_DIM
[Define dimension]
Definition: elastic.cpp:18
FTensor::Index< 'm', SPACE_DIM > m
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:219
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:483
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:442
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:242
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
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:572
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition: DMMoFEM.cpp:400
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:988
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
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
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
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
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:1042
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobina in TS solver.
Definition: TsCtx.cpp:131
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1045
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 getDMSubData(DM dm)
Get sub problem data structure.
Definition: DMMoFEM.hpp:1061
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.
Definition: Templates.hpp:135
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.
constexpr AssemblyType A
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
constexpr auto size_symm
Definition: plastic.cpp:42
static constexpr int approx_order
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
MoFEM::Interface & mField
MoFEMErrorCode runProblem()
[Run problem]
MoFEMErrorCode OPs()
[Boundary condition]
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
MoFEMErrorCode bC()
[Create common data]
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:25
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
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:47
Class (Function) to enforce essential constrains on the right hand side diagonal.
Definition: Essential.hpp:55
Class (Function) to enforce essential constrains.
Definition: Essential.hpp:39
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.
Definition: LogManager.cpp:298
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
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 filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed 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.
Clear Schur complement internal data.
Definition: Schur.hpp:22
Assemble Schur complement.
Definition: Schur.hpp:104
PipelineManager interface.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
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 refernce to pointer of interface.
MoFEMErrorCode preProcess()
function is run at the beginning of loop
boost::shared_ptr< MatrixDouble > vecFieldPtr
std::array< double, SPACE_DIM > 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, SPACE_DIM > 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)
[Push operators to pipeline]
Definition: plastic.cpp:755
SetUpSchur()=default
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, SmartPetscObj< DM > sub_dm, SmartPetscObj< IS > field_split_it, SmartPetscObj< AO > ao_map)
Create data structure for handling Schur complement.
Definition: plastic.cpp:1474
virtual MoFEMErrorCode setUp(SmartPetscObj< TS > solver)=0
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
SmartPetscObj< DM > subDM
field split sub dm
Definition: plastic.cpp:1296
MoFEMErrorCode setUp(KSP solver)
SmartPetscObj< Mat > S
SetUpSchurImpl(MoFEM::Interface &m_field)
SmartPetscObj< DM > createSubDM()
MoFEMErrorCode setUp(SmartPetscObj< TS > solver)
MoFEMErrorCode setPC(PC pc)
MoFEMErrorCode setOperator()
MoFEM::Interface & mField