v0.14.0
Loading...
Searching...
No Matches
elastic.cpp
Go to the documentation of this file.
1/**
2 * @file elastic.cpp
3 * @brief elastic example
4 * @version 0.13.2
5 * @date 2022-09-19
6 *
7 * @copyright Copyright (c) 2022
8 *
9 */
10
11#include <MoFEM.hpp>
12
13using namespace MoFEM;
14
15constexpr int BASE_DIM = 1; //< Dimension of the base functions
16
17//! [Define dimension]
18constexpr int SPACE_DIM =
19 EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
20//! [Define dimension]
22 ? AssemblyType::SCHUR
23 : AssemblyType::PETSC; //< selected assembly type
24
25constexpr IntegrationType I =
26 IntegrationType::GAUSS; //< selected integration type
27
28//! [Define entities]
35//! [Define entities]
36
37//! [OpK]
39 I>::OpGradSymTensorGrad<BASE_DIM, SPACE_DIM, SPACE_DIM, 0>;
40//! [OpK]
41//! [OpInternalForce]
43 I>::OpGradTimesSymTensor<BASE_DIM, SPACE_DIM, SPACE_DIM>;
44//! [OpInternalForce]
45struct DomainBCs {};
46struct BoundaryBCs {};
47
54
55template <int DIM> struct PostProcEleByDim;
56
57template <> struct PostProcEleByDim<2> {
61};
62
63template <> struct PostProcEleByDim<3> {
67};
68
72
73#include <ElasticSpring.hpp>
74#include <FluidLevel.hpp>
75#include <CalculateTraction.hpp>
76#include <NaturalDomainBC.hpp>
77#include <NaturalBoundaryBC.hpp>
78
79constexpr double young_modulus = 1;
80constexpr double poisson_ratio = 0.3;
81constexpr double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
82constexpr double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
83
84PetscBool is_plane_strain = PETSC_FALSE;
85
86struct Example {
87
88 Example(MoFEM::Interface &m_field) : mField(m_field) {}
89
91
92private:
94
95 PetscBool doEvalField;
96 std::array<double, SPACE_DIM> fieldEvalCoords;
97 boost::shared_ptr<FieldEvaluatorInterface::SetPtsData> fieldEvalData;
98 boost::shared_ptr<MatrixDouble> vectorFieldPtr;
99
107
109 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
110 std::string field_name, std::string block_name,
111 boost::shared_ptr<MatrixDouble> mat_D_Ptr, Sev sev);
112};
113
115 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
116 std::string field_name, std::string block_name,
117 boost::shared_ptr<MatrixDouble> mat_D_Ptr, Sev sev) {
119
120 struct OpMatBlocks : public DomainEleOp {
121 OpMatBlocks(std::string field_name, boost::shared_ptr<MatrixDouble> m,
122 double bulk_modulus_K, double shear_modulus_G,
123 MoFEM::Interface &m_field, Sev sev,
124 std::vector<const CubitMeshSets *> meshset_vec_ptr)
126 bulkModulusKDefault(bulk_modulus_K),
127 shearModulusGDefault(shear_modulus_G) {
128 std::fill(&(doEntities[MBEDGE]), &(doEntities[MBMAXTYPE]), false);
129 CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
130 "Can not get data from block");
131 }
132
133 MoFEMErrorCode doWork(int side, EntityType type,
136
137 for (auto &b : blockData) {
138
139 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
140 CHKERR getMatDPtr(matDPtr, b.bulkModulusK, b.shearModulusG);
142 }
143 }
144
145 CHKERR getMatDPtr(matDPtr, bulkModulusKDefault, shearModulusGDefault);
147 }
148
149 private:
150 boost::shared_ptr<MatrixDouble> matDPtr;
151
152 struct BlockData {
153 double bulkModulusK;
154 double shearModulusG;
155 Range blockEnts;
156 };
157
158 double bulkModulusKDefault;
159 double shearModulusGDefault;
160 std::vector<BlockData> blockData;
161
163 extractBlockData(MoFEM::Interface &m_field,
164 std::vector<const CubitMeshSets *> meshset_vec_ptr,
165 Sev sev) {
167
168 for (auto m : meshset_vec_ptr) {
169 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
170 std::vector<double> block_data;
171 CHKERR m->getAttributes(block_data);
172 if (block_data.size() < 2) {
173 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
174 "Expected that block has two attributes");
175 }
176 auto get_block_ents = [&]() {
177 Range ents;
178 CHKERR
179 m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
180 return ents;
181 };
182
183 double young_modulus = block_data[0];
184 double poisson_ratio = block_data[1];
185 double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
186 double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
187
188 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
189 << "E = " << young_modulus << " nu = " << poisson_ratio;
190
191 blockData.push_back(
192 {bulk_modulus_K, shear_modulus_G, get_block_ents()});
193 }
194 MOFEM_LOG_CHANNEL("WORLD");
196 }
197
198 MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
199 double bulk_modulus_K, double shear_modulus_G) {
201 //! [Calculate elasticity tensor]
202 auto set_material_stiffness = [&]() {
208 double A = 1.;
209 if (SPACE_DIM == 2 && !is_plane_strain) {
210 A = 2 * shear_modulus_G /
211 (bulk_modulus_K + (4. / 3.) * shear_modulus_G);
212 }
213 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
214 t_D(i, j, k, l) =
215 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
216 A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
217 t_kd(k, l);
218 };
219 //! [Calculate elasticity tensor]
220 constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
221 mat_D_ptr->resize(size_symm * size_symm, 1);
222 set_material_stiffness();
224 }
225 };
226
227 pipeline.push_back(new OpMatBlocks(
229
230 // Get blockset using regular expression
232
233 (boost::format("%s(.*)") % block_name).str()
234
235 ))
236
237 ));
238
240}
241
242//! [Run problem]
253}
254//! [Run problem]
255
256//! [Read mesh]
260 CHKERR simple->getOptions();
261 CHKERR simple->loadFile();
263}
264//! [Read mesh]
265
266//! [Set up problem]
270
271 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
272 const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
273 PetscInt choice_base_value = AINSWORTH;
274 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
275 LASBASETOPT, &choice_base_value, PETSC_NULL);
276
278 switch (choice_base_value) {
279 case AINSWORTH:
281 MOFEM_LOG("WORLD", Sev::inform)
282 << "Set AINSWORTH_LEGENDRE_BASE for displacements";
283 break;
284 case DEMKOWICZ:
286 MOFEM_LOG("WORLD", Sev::inform)
287 << "Set DEMKOWICZ_JACOBI_BASE for displacements";
288 break;
289 default:
290 base = LASTBASE;
291 break;
292 }
293
294 // Add field
295 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
296 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
297 int order = 3;
298 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
299
300 CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
301
302 CHKERR simple->setFieldOrder("U", order);
303 CHKERR simple->setFieldOrder("GEOMETRY", 2);
304 CHKERR simple->setUp();
305
306 auto project_ho_geometry = [&]() {
307 Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
308 return mField.loop_dofs("GEOMETRY", ent_method);
309 };
310 CHKERR project_ho_geometry();
311
312 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-plane_strain", &is_plane_strain,
313 PETSC_NULL);
314
315 int coords_dim = SPACE_DIM;
316 CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
317 fieldEvalCoords.data(), &coords_dim,
318 &doEvalField);
319
320 if (doEvalField) {
321 vectorFieldPtr = boost::make_shared<MatrixDouble>();
323 mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
324
325 if constexpr (SPACE_DIM == 3) {
327 fieldEvalData, simple->getDomainFEName());
328 } else {
330 fieldEvalData, simple->getDomainFEName());
331 }
332
333 fieldEvalData->setEvalPoints(fieldEvalCoords.data(), 1);
334 auto no_rule = [](int, int, int) { return -1; };
335 auto field_eval_fe_ptr = fieldEvalData->feMethodPtr.lock();
336 field_eval_fe_ptr->getRuleHook = no_rule;
337
338 field_eval_fe_ptr->getOpPtrVector().push_back(
340 }
341
343}
344//! [Set up problem]
345
346//! [Boundary condition]
349 auto pip = mField.getInterface<PipelineManager>();
351 auto bc_mng = mField.getInterface<BcManager>();
352
353 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
354 "U", 0, 0);
355 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
356 "U", 1, 1);
357 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
358 "U", 2, 2);
359 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
360 "REMOVE_ALL", "U", 0, 3);
361 CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
362 simple->getProblemName(), "U");
363
364 // adding MPCs
365 CHKERR bc_mng->addBlockDOFsToMPCs(simple->getProblemName(), "U");
366
368}
369//! [Boundary condition]
370
371//! [Push operators to pipeline]
374 auto pip = mField.getInterface<PipelineManager>();
376 auto bc_mng = mField.getInterface<BcManager>();
377
378 //! [Integration rule]
379 auto integration_rule = [](int, int, int approx_order) {
380 return 2 * approx_order + 1;
381 };
382
383 auto integration_rule_bc = [](int, int, int approx_order) {
384 return 2 * approx_order + 1;
385 };
386
387 CHKERR pip->setDomainRhsIntegrationRule(integration_rule);
388 CHKERR pip->setDomainLhsIntegrationRule(integration_rule);
389 CHKERR pip->setBoundaryRhsIntegrationRule(integration_rule_bc);
390 CHKERR pip->setBoundaryLhsIntegrationRule(integration_rule_bc);
391 //! [Integration rule]
392
394 pip->getOpDomainLhsPipeline(), {H1}, "GEOMETRY");
396 pip->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
398 pip->getOpBoundaryRhsPipeline(), {NOSPACE}, "GEOMETRY");
400 pip->getOpBoundaryLhsPipeline(), {NOSPACE}, "GEOMETRY");
401
402
403 //! [Push domain stiffness matrix to pipeline]
404 auto mat_D_ptr = boost::make_shared<MatrixDouble>();
405
406 // Assemble domain stiffness matrix
407 CHKERR addMatBlockOps(pip->getOpDomainLhsPipeline(), "U", "MAT_ELASTIC",
408 mat_D_ptr, Sev::verbose);
409 pip->getOpDomainLhsPipeline().push_back(new OpK("U", "U", mat_D_ptr));
410 //! [Push domain stiffness matrix to pipeline]
411
412 //! [Push Internal forces]
413 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
414 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
415 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
416 pip->getOpDomainRhsPipeline().push_back(
418 mat_grad_ptr));
419 CHKERR addMatBlockOps(pip->getOpDomainRhsPipeline(), "U", "MAT_ELASTIC",
420 mat_D_ptr, Sev::inform);
421 pip->getOpDomainRhsPipeline().push_back(
422 new OpSymmetrizeTensor<SPACE_DIM>(mat_grad_ptr, mat_strain_ptr));
423 pip->getOpDomainRhsPipeline().push_back(
425 mat_strain_ptr, mat_stress_ptr, mat_D_ptr));
426
427 pip->getOpDomainRhsPipeline().push_back(
428 new OpInternalForce("U", mat_stress_ptr,
429 [](double, double, double) constexpr { return -1; }));
430 //! [Push Internal forces]
431
432 //! [Push Body forces]
434 pip->getOpDomainRhsPipeline(), mField, "U", Sev::inform);
435 //! [Push Body forces]
436
437 //! [Push natural boundary conditions]
438 // Add force boundary condition
440 pip->getOpBoundaryRhsPipeline(), mField, "U", -1, Sev::inform);
441 // Add case for mix type of BCs
443 pip->getOpBoundaryLhsPipeline(), mField, "U", Sev::verbose);
444 //! [Push natural boundary conditions]
446}
447//! [Push operators to pipeline]
448
449
450struct SetUpSchur {
451 static boost::shared_ptr<SetUpSchur>
454
455protected:
456 SetUpSchur() = default;
457};
458//! [Solve]
462 auto pip = mField.getInterface<PipelineManager>();
463 auto solver = pip->createKSP();
464 CHKERR KSPSetFromOptions(solver);
465
466 auto dm = simple->getDM();
467 auto D = createDMVector(dm);
468 auto F = vectorDuplicate(D);
469
470 auto set_essential_bc = [&]() {
472 // This is low level pushing finite elements (pipelines) to solver
473 auto ksp_ctx_ptr = getDMKspCtx(simple->getDM());
474
475 auto pre_proc_rhs = boost::make_shared<FEMethod>();
476 auto post_proc_rhs = boost::make_shared<FEMethod>();
477 auto post_proc_lhs = boost::make_shared<FEMethod>();
478
479 auto get_pre_proc_hook = [&]() {
481 {});
482 };
483 pre_proc_rhs->preProcessHook = get_pre_proc_hook();
484
485 auto get_post_proc_hook_rhs = [this, post_proc_rhs]() {
487
489 1.)();
490 CHKERR EssentialPostProcRhs<MPCsType>(mField, post_proc_rhs, 1.)();
492 };
493
494 auto get_post_proc_hook_lhs = [this, post_proc_lhs]() {
496
498 1.)();
501 };
502
503 post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
504 post_proc_lhs->postProcessHook = get_post_proc_hook_lhs;
505
506 ksp_ctx_ptr->getPreProcComputeRhs().push_front(pre_proc_rhs);
507 ksp_ctx_ptr->getPostProcComputeRhs().push_back(post_proc_rhs);
508 ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_lhs);
510 };
511
512 auto setup_and_solve = [&]() {
514 BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
515 MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp";
516 CHKERR KSPSetUp(solver);
517 MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp <= Done";
518 MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve";
519 CHKERR KSPSolve(solver, F, D);
520 MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve <= Done";
522 };
523
524 MOFEM_LOG_CHANNEL("TIMER");
525 MOFEM_LOG_TAG("TIMER", "timer");
526
527 CHKERR set_essential_bc();
528
529 if (A == AssemblyType::SCHUR) {
530 auto schur_ptr = SetUpSchur::createSetUpSchur(mField);
531 CHKERR schur_ptr->setUp(solver);
532 CHKERR setup_and_solve();
533 } else {
534 CHKERR setup_and_solve();
535 }
536
537 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
538 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
539 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
540
541 if (doEvalField) {
542 if constexpr (SPACE_DIM == 3) {
543 CHKERR mField.getInterface<FieldEvaluatorInterface>()->evalFEAtThePoint3D(
544 fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
545 simple->getDomainFEName(), fieldEvalData, mField.get_comm_rank(),
546 mField.get_comm_rank(), nullptr, MF_EXIST, QUIET);
547 } else {
548 CHKERR mField.getInterface<FieldEvaluatorInterface>()->evalFEAtThePoint2D(
549 fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
550 simple->getDomainFEName(), fieldEvalData, mField.get_comm_rank(),
551 mField.get_comm_rank(), nullptr, MF_EXIST, QUIET);
552 }
553
554 if (vectorFieldPtr->size1()) {
555 auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vectorFieldPtr);
556 if constexpr (SPACE_DIM == 2)
557 MOFEM_LOG("FieldEvaluator", Sev::inform)
558 << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1);
559 else
560 MOFEM_LOG("FieldEvaluator", Sev::inform)
561 << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1)
562 << " U_Z: " << t_disp(2);
563 }
564
566 }
568}
569//! [Solve]
570
571//! [Postprocess results]
575 auto pip = mField.getInterface<PipelineManager>();
576 auto det_ptr = boost::make_shared<VectorDouble>();
577 auto jac_ptr = boost::make_shared<MatrixDouble>();
578 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
579 //! [Postprocess clean]
580 pip->getDomainRhsFE().reset();
581 pip->getDomainLhsFE().reset();
582 pip->getBoundaryRhsFE().reset();
583 pip->getBoundaryLhsFE().reset();
584 //! [Postprocess clean]
585
586 //! [Postprocess initialise]
587 auto post_proc_mesh = boost::make_shared<moab::Core>();
588 auto post_proc_begin = boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(
589 mField, post_proc_mesh);
590 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
591 mField, post_proc_mesh);
592 //! [Postprocess initialise]
593
594 auto calculate_stress_ops = [&](auto &pip) {
596 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
597 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
598 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
600 "U", mat_grad_ptr));
601 auto mat_D_ptr = boost::make_shared<MatrixDouble>();
602 CHKERR addMatBlockOps(pip, "U", "MAT_ELASTIC", mat_D_ptr, Sev::verbose);
603 pip.push_back(
604 new OpSymmetrizeTensor<SPACE_DIM>(mat_grad_ptr, mat_strain_ptr));
606 mat_strain_ptr, mat_stress_ptr, mat_D_ptr));
607 auto u_ptr = boost::make_shared<MatrixDouble>();
608 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
609 auto x_ptr = boost::make_shared<MatrixDouble>();
610 pip.push_back(
611 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
612 return boost::make_tuple(u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr);
613 };
614
615 auto post_proc_domain = [&](auto post_proc_mesh) {
616 auto post_proc_fe =
617 boost::make_shared<PostProcEleDomain>(mField, post_proc_mesh);
619
620 auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr] =
621 calculate_stress_ops(post_proc_fe->getOpPtrVector());
622
623 post_proc_fe->getOpPtrVector().push_back(
624
625 new OpPPMap(
626
627 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
628
629 {},
630
631 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
632
633 {},
634
635 {{"STRAIN", mat_strain_ptr}, {"STRESS", mat_stress_ptr}}
636
637 )
638
639 );
640 return post_proc_fe;
641 };
642
643 auto post_proc_boundary = [&](auto post_proc_mesh) {
644 auto post_proc_fe =
645 boost::make_shared<PostProcEleBdy>(mField, post_proc_mesh);
647 post_proc_fe->getOpPtrVector(), {}, "GEOMETRY");
648 auto op_loop_side =
649 new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
650 // push ops to side element, through op_loop_side operator
651 auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr] =
652 calculate_stress_ops(op_loop_side->getOpPtrVector());
653 post_proc_fe->getOpPtrVector().push_back(op_loop_side);
654 auto mat_traction_ptr = boost::make_shared<MatrixDouble>();
655 post_proc_fe->getOpPtrVector().push_back(
656 new ElasticExample::OpCalculateTraction(mat_stress_ptr,
657 mat_traction_ptr));
658
660
661 post_proc_fe->getOpPtrVector().push_back(
662
663 new OpPPMap(
664
665 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
666
667 {},
668
669 {{"U", u_ptr}, {"GEOMETRY", x_ptr}, {"T", mat_traction_ptr}},
670
671 {},
672
673 {{"STRAIN", mat_strain_ptr}, {"STRESS", mat_stress_ptr}}
674
675 )
676
677 );
678 return post_proc_fe;
679 };
680
681 PetscBool post_proc_skin_only = PETSC_FALSE;
682 if (SPACE_DIM == 3) {
683 post_proc_skin_only = PETSC_TRUE;
684 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_skin_only",
685 &post_proc_skin_only, PETSC_NULL);
686 }
687 if (post_proc_skin_only == PETSC_FALSE) {
688 pip->getDomainRhsFE() = post_proc_domain(post_proc_mesh);
689 } else {
690 pip->getBoundaryRhsFE() = post_proc_boundary(post_proc_mesh);
691 }
692
694 post_proc_begin->getFEMethod());
695 CHKERR pip->loopFiniteElements();
697 post_proc_end->getFEMethod());
698
699 CHKERR post_proc_end->writeFile("out_elastic.h5m");
701}
702//! [Postprocess results]
703
704//! [Check]
706 MOFEM_LOG_CHANNEL("WORLD");
708 auto pip = mField.getInterface<PipelineManager>();
710 pip->getDomainRhsFE().reset();
711 pip->getDomainLhsFE().reset();
712 pip->getBoundaryRhsFE().reset();
713 pip->getBoundaryLhsFE().reset();
714
715 auto integration_rule = [](int, int, int p_data) { return 2 * p_data + 1; };
716 CHKERR pip->setDomainRhsIntegrationRule(integration_rule);
717 CHKERR pip->setBoundaryRhsIntegrationRule(integration_rule);
718
720 pip->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
722 pip->getOpBoundaryRhsPipeline(), {}, "GEOMETRY");
723
724 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
725 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
726 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
727
728 pip->getOpDomainRhsPipeline().push_back(
730 mat_grad_ptr));
731 pip->getOpDomainRhsPipeline().push_back(
732 new OpSymmetrizeTensor<SPACE_DIM>(mat_grad_ptr, mat_strain_ptr));
733
734 auto mat_D_ptr = boost::make_shared<MatrixDouble>();
735 CHKERR addMatBlockOps(pip->getOpDomainRhsPipeline(), "U", "MAT_ELASTIC",
736 mat_D_ptr, Sev::verbose);
737 pip->getOpDomainRhsPipeline().push_back(
739 mat_strain_ptr, mat_stress_ptr, mat_D_ptr));
740
741 pip->getOpDomainRhsPipeline().push_back(
742 new OpInternalForce("U", mat_stress_ptr));
743
744 pip->getOpBoundaryRhsPipeline().push_back(
746 mat_strain_ptr, mat_stress_ptr, mat_D_ptr));
748 pip->getOpDomainRhsPipeline(), mField, "U", Sev::verbose);
750 pip->getOpBoundaryRhsPipeline(), mField, "U", -1, Sev::verbose);
751
752 auto dm = simple->getDM();
753 auto res = createDMVector(dm);
754 pip->getDomainRhsFE()->f = res;
755 pip->getBoundaryRhsFE()->f = res;
756
757 CHKERR VecZeroEntries(res);
758
759 CHKERR mField.getInterface<FieldBlas>()->fieldScale(-1, "U");
760 CHKERR pip->loopFiniteElements();
761 CHKERR mField.getInterface<FieldBlas>()->fieldScale(-1, "U");
762
763 CHKERR VecGhostUpdateBegin(res, ADD_VALUES, SCATTER_REVERSE);
764 CHKERR VecGhostUpdateEnd(res, ADD_VALUES, SCATTER_REVERSE);
765 CHKERR VecAssemblyBegin(res);
766 CHKERR VecAssemblyEnd(res);
767
768 auto zero_residual_at_constrains = [&]() {
770 auto fe_post_proc_ptr = boost::make_shared<FEMethod>();
771 auto get_post_proc_hook_rhs = [this, fe_post_proc_ptr, res]() {
774 mField, fe_post_proc_ptr, res)();
776 mField, fe_post_proc_ptr, 0, res)();
778 mField, fe_post_proc_ptr, 0, res)();
780 };
781 fe_post_proc_ptr->postProcessHook = get_post_proc_hook_rhs;
782 CHKERR DMoFEMPostProcessFiniteElements(dm, fe_post_proc_ptr.get());
784 };
785
786 CHKERR zero_residual_at_constrains();
787
788 double nrm2;
789 CHKERR VecNorm(res, NORM_2, &nrm2);
790 MOFEM_LOG_CHANNEL("WORLD");
791 MOFEM_LOG_C("WORLD", Sev::inform, "residual = %3.4e\n", nrm2);
792
793 PetscBool test = PETSC_FALSE;
794 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test, PETSC_NULL);
795 if (test == PETSC_TRUE) {
796
797 auto post_proc_residual = [&](auto dm, auto f_res, auto out_name) {
799 auto post_proc_fe =
800 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
802 auto u_vec = boost::make_shared<MatrixDouble>();
803 post_proc_fe->getOpPtrVector().push_back(
804 new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_vec, f_res));
805 post_proc_fe->getOpPtrVector().push_back(
806
807 new OpPPMap(
808
809 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
810
811 {},
812
813 {{"RES", u_vec}},
814
815 {}, {})
816
817 );
818
819 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
820 post_proc_fe);
821 post_proc_fe->writeFile(out_name);
823 };
824
825 CHKERR post_proc_residual(simple->getDM(), res, "res.h5m");
826
827 constexpr double eps = 1e-8;
828 if (nrm2 > eps)
829 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
830 "Residual is not zero");
831 }
832
834}
835//! [Check]
836
837static char help[] = "...\n\n";
838
839int main(int argc, char *argv[]) {
840
841 // Initialisation of MoFEM/PETSc and MOAB data structures
842 const char param_file[] = "param_file.petsc";
844
845 auto core_log = logging::core::get();
846 core_log->add_sink(
848
849 core_log->add_sink(
851 LogManager::setLog("FieldEvaluator");
852 MOFEM_LOG_TAG("FieldEvaluator", "field_eval");
853
854 try {
855
856 //! [Register MoFEM discrete manager in PETSc]
857 DMType dm_name = "DMMOFEM";
858 CHKERR DMRegister_MoFEM(dm_name);
859 //! [Register MoFEM discrete manager in PETSc
860
861 //! [Create MoAB]
862 moab::Core mb_instance; ///< mesh database
863 moab::Interface &moab = mb_instance; ///< mesh database interface
864 //! [Create MoAB]
865
866 //! [Create MoFEM]
867 MoFEM::Core core(moab); ///< finite element database
868 MoFEM::Interface &m_field = core; ///< finite element database interface
869 //! [Create MoFEM]
870
871 //! [Example]
872 Example ex(m_field);
873 CHKERR ex.runProblem();
874 //! [Example]
875 }
877
879}
880
881struct SetUpSchurImpl : public SetUpSchur {
882
884 if (S) {
887 "Is expected that schur matrix is not allocated. This is "
888 "possible only is PC is set up twice");
889 }
890 }
891 virtual ~SetUpSchurImpl() { S.reset(); }
892
894
895private:
900
902
904
908};
909
912 auto pip = mField.getInterface<PipelineManager>();
913 PC pc;
914 CHKERR KSPGetPC(solver, &pc);
915 PetscBool is_pcfs = PETSC_FALSE;
916 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
917 if (is_pcfs) {
918 if (S) {
921 "Is expected that schur matrix is not allocated. This is "
922 "possible only is PC is set up twice");
923 }
927 CHKERR MatSetBlockSize(S, SPACE_DIM);
928 CHKERR MatSetOption(S, MAT_SYMMETRIC, PETSC_TRUE);
930 CHKERR setPC(pc);
931 } else {
932 pip->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
933 pip->getOpBoundaryLhsPipeline().push_back(
934 new OpSchurAssembleEnd<SCHUR_DSYSV>({}, {}, {}, {}, {}));
935 pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
936 pip->getOpDomainLhsPipeline().push_back(
937 new OpSchurAssembleEnd<SCHUR_DSYSV>({}, {}, {}, {}, {}));
938 }
940}
941
945 CHKERR mField.get_moab().get_entities_by_dimension(simple->getMeshset(),
947 CHKERR mField.get_moab().get_entities_by_handle(simple->getMeshset(),
948 subEnts);
949 subEnts = subtract(subEnts, volEnts);
951};
952
956 subDM = createDM(mField.get_comm(), "DMMOFEM");
957 CHKERR DMMoFEMCreateSubDM(subDM, simple->getDM(), "SUB");
959 CHKERR DMMoFEMAddElement(subDM, simple->getDomainFEName());
960 auto sub_ents_ptr = boost::make_shared<Range>(subEnts);
961 CHKERR DMMoFEMAddSubFieldRow(subDM, "U", sub_ents_ptr);
962 CHKERR DMMoFEMAddSubFieldCol(subDM, "U", sub_ents_ptr);
963 CHKERR DMSetUp(subDM);
965}
966
969 auto pip = mField.getInterface<PipelineManager>();
970
971 // Boundary
972 auto dm_is = getDMSubData(subDM)->getSmartRowIs();
973 auto ao_up = createAOMappingIS(dm_is, PETSC_NULL);
974 pip->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
975 pip->getOpBoundaryLhsPipeline().push_back(new OpSchurAssembleEnd<SCHUR_DSYSV>(
976 {"U"}, {boost::make_shared<Range>(volEnts)}, {ao_up}, {S}, {true}));
977 // Domain
978 pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
979 pip->getOpDomainLhsPipeline().push_back(new OpSchurAssembleEnd<SCHUR_DSYSV>(
980 {"U"}, {boost::make_shared<Range>(volEnts)}, {ao_up}, {S}, {true}));
981
982 auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
983 auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
984
985 pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
987 if (S) {
988 CHKERR MatZeroEntries(S);
989 }
990 MOFEM_LOG("TIMER", Sev::inform) << "Lhs Assemble Begin";
992 };
993
994 post_proc_schur_lhs_ptr->postProcessHook = [this, post_proc_schur_lhs_ptr,
995 ao_up]() {
997 CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
998 CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
1000 mField, post_proc_schur_lhs_ptr, 1, S, ao_up)();
1001 MOFEM_LOG("TIMER", Sev::inform) << "Lhs Assemble End";
1003 };
1004
1005 auto simple = mField.getInterface<Simple>();
1006 auto ksp_ctx_ptr = getDMKspCtx(simple->getDM());
1007
1008 ksp_ctx_ptr->getPreProcSetOperators().push_front(pre_proc_schur_lhs_ptr);
1009 ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_schur_lhs_ptr);
1010
1012}
1013
1016 auto simple = mField.getInterface<Simple>();
1017 SmartPetscObj<IS> vol_is;
1018 mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1019 simple->getProblemName(), ROW, "U", 0, SPACE_DIM, vol_is, &volEnts);
1020 CHKERR PCFieldSplitSetIS(pc, NULL, vol_is);
1021 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1023}
1024
1025boost::shared_ptr<SetUpSchur>
1027 return boost::shared_ptr<SetUpSchur>(new SetUpSchurImpl(m_field));
1028}
Calculate traction for linear problem.
std::string param_file
Implementation of elastic spring bc.
Natural boundary condition applying pressure from fluid.
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
Implementation of natural boundary conditions.
Boundary conditions in domain, i.e. body forces.
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
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
PostProcEleByDim< SPACE_DIM >::SideEle SideEle
PostProcEleByDim< SPACE_DIM >::PostProcEleBdy PostProcEleBdy
static const double eps
constexpr int SPACE_DIM
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
#define SCHUR_ASSEMBLE
Definition: contact.cpp:18
@ 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
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
constexpr int order
double bulk_modulus_K
double shear_modulus_G
static char help[]
[Check]
Definition: elastic.cpp:837
PetscBool is_plane_strain
Definition: elastic.cpp:84
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesSymTensor< BASE_DIM, SPACE_DIM, SPACE_DIM > OpInternalForce
Definition: elastic.cpp:43
constexpr int SPACE_DIM
[Define dimension]
Definition: elastic.cpp:18
constexpr double poisson_ratio
Definition: elastic.cpp:80
constexpr int BASE_DIM
Definition: elastic.cpp:15
constexpr double shear_modulus_G
Definition: elastic.cpp:82
constexpr IntegrationType I
Definition: elastic.cpp:25
constexpr double bulk_modulus_K
Definition: elastic.cpp:81
constexpr AssemblyType A
[Define dimension]
Definition: elastic.cpp:21
constexpr double young_modulus
Definition: elastic.cpp:79
FTensor::Index< 'm', SPACE_DIM > m
@ F
auto integration_rule
constexpr auto t_kd
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 DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:542
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 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 DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:275
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:988
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:532
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
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
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.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
FTensor::Index< 'i', SPACE_DIM > i
double D
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
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 PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
auto getDMKspCtx(DM dm)
Get KSP context data structure used by DM.
Definition: DMMoFEM.hpp:1017
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.
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
constexpr AssemblyType A
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double young_modulus
Young modulus.
Definition: plastic.cpp:172
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
constexpr auto size_symm
Definition: plastic.cpp:42
constexpr auto field_name
static constexpr int approx_order
[OpInternalForce]
Definition: elastic.cpp:45
[Example]
Definition: plastic.cpp:228
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
PetscBool doEvalField
Definition: elastic.cpp:95
FieldApproximationBase base
Definition: plot_base.cpp:68
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
Example(MoFEM::Interface &m_field)
Definition: elastic.cpp:88
MoFEMErrorCode runProblem()
MoFEM::Interface & mField
Definition: plastic.cpp:235
std::array< double, SPACE_DIM > fieldEvalCoords
Definition: elastic.cpp:96
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
MoFEMErrorCode addMatBlockOps(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string field_name, std::string block_name, boost::shared_ptr< MatrixDouble > mat_D_Ptr, Sev sev)
Definition: elastic.cpp:114
boost::shared_ptr< MatrixDouble > vectorFieldPtr
Definition: elastic.cpp:98
boost::shared_ptr< FieldEvaluatorInterface::SetPtsData > fieldEvalData
Definition: elastic.cpp:97
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
virtual moab::Interface & get_moab()=0
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:76
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
Class (Function) to calculate residual side diagonal.
Definition: Essential.hpp:49
Basic algebra on fields.
Definition: FieldBlas.hpp:21
Field evaluator interface.
@ OPROW
operator doWork function is executed on FE rows
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
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
Definition: LogManager.cpp:348
Interface for managing meshsets containing materials and boundary conditions.
Assembly methods.
Definition: Natural.hpp:65
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
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
PipelineManager::ElementsAndOpsByDim< 2 >::FaceSideEle SideEle
PipelineManager::ElementsAndOpsByDim< 3 >::FaceSideEle SideEle
[Push operators to pipeline]
Definition: plastic.cpp:762
SetUpSchur()=default
virtual MoFEMErrorCode setUp(SmartPetscObj< KSP > solver)=0
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:1482
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
SmartPetscObj< DM > subDM
field split sub dm
Definition: plastic.cpp:1304
MoFEMErrorCode setUp(KSP solver)
SmartPetscObj< Mat > S
SetUpSchurImpl(MoFEM::Interface &m_field)
Definition: elastic.cpp:883
MoFEMErrorCode setEntities()
virtual ~SetUpSchurImpl()
Definition: elastic.cpp:891
MoFEMErrorCode setPC(PC pc)
MoFEMErrorCode setUpSubDM()
Definition: elastic.cpp:953
MoFEMErrorCode setOperator()
MoFEM::Interface & mField