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
15//! [Define dimension]
16constexpr int BASE_DIM = 1; //< Dimension of the base functions
17constexpr int SPACE_DIM =
18 EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
19constexpr AssemblyType A = (SCHUR_ASSEMBLE)
20 ? AssemblyType::SCHUR
21 : AssemblyType::PETSC; //< selected assembly type
22
23constexpr IntegrationType I =
24 IntegrationType::GAUSS; //< selected integration type
25
32
34 I>::OpGradSymTensorGrad<BASE_DIM, SPACE_DIM, SPACE_DIM, 0>;
36 I>::OpGradTimesSymTensor<BASE_DIM, SPACE_DIM, SPACE_DIM>;
37
38struct DomainBCs {};
39struct BoundaryBCs {};
40
47
48template <int DIM> struct PostProcEleByDim;
49
50template <> struct PostProcEleByDim<2> {
54};
55
56template <> struct PostProcEleByDim<3> {
60};
61
65
66#include <ElasticSpring.hpp>
67#include <CalculateTraction.hpp>
68#include <NaturalDomainBC.hpp>
69#include <NaturalBoundaryBC.hpp>
70
71constexpr double young_modulus = 1;
72constexpr double poisson_ratio = 0.3;
73constexpr double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
74constexpr double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
75
76PetscBool is_plain_strain = PETSC_FALSE;
77
78struct Example {
79
80 Example(MoFEM::Interface &m_field) : mField(m_field) {}
81
83
84private:
86
87 PetscBool doEvalField;
88 std::array<double, SPACE_DIM> fieldEvalCoords;
89 boost::shared_ptr<FieldEvaluatorInterface::SetPtsData> fieldEvalData;
90 boost::shared_ptr<MatrixDouble> vectorFieldPtr;
91
99
101 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
102 std::string field_name, std::string block_name,
103 boost::shared_ptr<MatrixDouble> mat_D_Ptr, Sev sev);
104};
105
107 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
108 std::string field_name, std::string block_name,
109 boost::shared_ptr<MatrixDouble> mat_D_Ptr, Sev sev) {
111
112 struct OpMatBlocks : public DomainEleOp {
113 OpMatBlocks(std::string field_name, boost::shared_ptr<MatrixDouble> m,
114 double bulk_modulus_K, double shear_modulus_G,
115 MoFEM::Interface &m_field, Sev sev,
116 std::vector<const CubitMeshSets *> meshset_vec_ptr)
118 bulkModulusKDefault(bulk_modulus_K),
119 shearModulusGDefault(shear_modulus_G) {
120 std::fill(&(doEntities[MBEDGE]), &(doEntities[MBMAXTYPE]), false);
121 CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
122 "Can not get data from block");
123 }
124
125 MoFEMErrorCode doWork(int side, EntityType type,
128
129 for (auto &b : blockData) {
130
131 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
132 CHKERR getMatDPtr(matDPtr, b.bulkModulusK, b.shearModulusG);
134 }
135 }
136
137 CHKERR getMatDPtr(matDPtr, bulkModulusKDefault, shearModulusGDefault);
139 }
140
141 private:
142 boost::shared_ptr<MatrixDouble> matDPtr;
143
144 struct BlockData {
145 double bulkModulusK;
146 double shearModulusG;
147 Range blockEnts;
148 };
149
150 double bulkModulusKDefault;
151 double shearModulusGDefault;
152 std::vector<BlockData> blockData;
153
155 extractBlockData(MoFEM::Interface &m_field,
156 std::vector<const CubitMeshSets *> meshset_vec_ptr,
157 Sev sev) {
159
160 for (auto m : meshset_vec_ptr) {
161 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
162 std::vector<double> block_data;
163 CHKERR m->getAttributes(block_data);
164 if (block_data.size() < 2) {
165 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
166 "Expected that block has two attributes");
167 }
168 auto get_block_ents = [&]() {
169 Range ents;
170 CHKERR
171 m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
172 return ents;
173 };
174
175 double young_modulus = block_data[0];
176 double poisson_ratio = block_data[1];
177 double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
178 double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
179
180 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
181 << "E = " << young_modulus << " nu = " << poisson_ratio;
182
183 blockData.push_back(
184 {bulk_modulus_K, shear_modulus_G, get_block_ents()});
185 }
186 MOFEM_LOG_CHANNEL("WORLD");
188 }
189
190 MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
191 double bulk_modulus_K, double shear_modulus_G) {
193 //! [Calculate elasticity tensor]
194 auto set_material_stiffness = [&]() {
200 double A = 1.;
201 if (SPACE_DIM == 2 && !is_plain_strain) {
202 A = 2 * shear_modulus_G /
203 (bulk_modulus_K + (4. / 3.) * shear_modulus_G);
204 }
205 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
206 t_D(i, j, k, l) =
207 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
208 A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
209 t_kd(k, l);
210 };
211 //! [Calculate elasticity tensor]
212 constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
213 mat_D_ptr->resize(size_symm * size_symm, 1);
214 set_material_stiffness();
216 }
217 };
218
219 pipeline.push_back(new OpMatBlocks(
221
222 // Get blockset using regular expression
224
225 (boost::format("%s(.*)") % block_name).str()
226
227 ))
228
229 ));
230
232}
233
234//! [Run problem]
245}
246//! [Run problem]
247
248//! [Read mesh]
252 CHKERR simple->getOptions();
253 CHKERR simple->loadFile();
255}
256//! [Read mesh]
257
258//! [Set up problem]
262
263 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
264 const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
265 PetscInt choice_base_value = AINSWORTH;
266 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
267 LASBASETOPT, &choice_base_value, PETSC_NULL);
268
270 switch (choice_base_value) {
271 case AINSWORTH:
273 MOFEM_LOG("WORLD", Sev::inform)
274 << "Set AINSWORTH_LEGENDRE_BASE for displacements";
275 break;
276 case DEMKOWICZ:
278 MOFEM_LOG("WORLD", Sev::inform)
279 << "Set DEMKOWICZ_JACOBI_BASE for displacements";
280 break;
281 default:
282 base = LASTBASE;
283 break;
284 }
285
286 // Add field
287 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
288 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
289 int order = 3;
290 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
291
292 CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
293
294 CHKERR simple->setFieldOrder("U", order);
295 CHKERR simple->setFieldOrder("GEOMETRY", 2);
296 CHKERR simple->setUp();
297
298 auto project_ho_geometry = [&]() {
299 Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
300 return mField.loop_dofs("GEOMETRY", ent_method);
301 };
302 CHKERR project_ho_geometry();
303
304 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-plane_strain", &is_plain_strain,
305 PETSC_NULL);
306
307 int coords_dim = SPACE_DIM;
308 CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
309 fieldEvalCoords.data(), &coords_dim,
310 &doEvalField);
311
312 if (doEvalField) {
313 vectorFieldPtr = boost::make_shared<MatrixDouble>();
315 mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
316
317 if constexpr (SPACE_DIM == 3) {
319 fieldEvalData, simple->getDomainFEName());
320 } else {
322 fieldEvalData, simple->getDomainFEName());
323 }
324
325 fieldEvalData->setEvalPoints(fieldEvalCoords.data(), 1);
326 auto no_rule = [](int, int, int) { return -1; };
327 auto field_eval_fe_ptr = fieldEvalData->feMethodPtr.lock();
328 field_eval_fe_ptr->getRuleHook = no_rule;
329
330 field_eval_fe_ptr->getOpPtrVector().push_back(
332 }
333
335}
336//! [Set up problem]
337
338//! [Boundary condition]
341 auto pip = mField.getInterface<PipelineManager>();
343 auto bc_mng = mField.getInterface<BcManager>();
344
345 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
346 "U", 0, 0);
347 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
348 "U", 1, 1);
349 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
350 "U", 2, 2);
351 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
352 "REMOVE_ALL", "U", 0, 3);
353 CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
354 simple->getProblemName(), "U");
355
356 auto integration_rule = [](int, int, int approx_order) {
357 return 2 * approx_order + 1;
358 };
359
360 auto integration_rule_bc = [](int, int, int approx_order) {
361 return 2 * approx_order + 1;
362 };
363
364 CHKERR pip->setDomainRhsIntegrationRule(integration_rule);
365 CHKERR pip->setDomainLhsIntegrationRule(integration_rule);
366 CHKERR pip->setBoundaryRhsIntegrationRule(integration_rule_bc);
367 CHKERR pip->setBoundaryLhsIntegrationRule(integration_rule_bc);
368
370}
371//! [Boundary condition]
372
373//! [Push operators to pipeline]
376 auto pip = mField.getInterface<PipelineManager>();
378 auto bc_mng = mField.getInterface<BcManager>();
379
381 pip->getOpDomainLhsPipeline(), {H1}, "GEOMETRY");
383 pip->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
385 pip->getOpBoundaryRhsPipeline(), {NOSPACE}, "GEOMETRY");
387 pip->getOpBoundaryLhsPipeline(), {NOSPACE}, "GEOMETRY");
388
389 auto mat_D_ptr = boost::make_shared<MatrixDouble>();
390
391 // Assemble domain stiffness matrix
392 CHKERR addMatBlockOps(pip->getOpDomainLhsPipeline(), "U", "MAT_ELASTIC",
393 mat_D_ptr, Sev::verbose);
394 pip->getOpDomainLhsPipeline().push_back(new OpK("U", "U", mat_D_ptr));
395
396 // Infernal forces
397 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
398 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
399 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
400 pip->getOpDomainRhsPipeline().push_back(
402 mat_grad_ptr));
403 pip->getOpDomainRhsPipeline().push_back(
404 new OpSymmetrizeTensor<SPACE_DIM>("U", mat_grad_ptr, mat_strain_ptr));
405 CHKERR addMatBlockOps(pip->getOpDomainRhsPipeline(), "U", "MAT_ELASTIC",
406 mat_D_ptr, Sev::inform);
407 pip->getOpDomainRhsPipeline().push_back(
409 "U", mat_strain_ptr, mat_stress_ptr, mat_D_ptr));
410 // Internal forces
411 pip->getOpDomainRhsPipeline().push_back(
412 new OpInternalForce("U", mat_stress_ptr,
413 [](double, double, double) constexpr { return -1; }));
414
415 // Body forces
417 pip->getOpDomainRhsPipeline(), mField, "U", Sev::inform);
418
419 // Add force boundary condition
421 pip->getOpBoundaryRhsPipeline(), mField, "U", -1, Sev::inform);
422 // Add case for mix type of BCs
424 pip->getOpBoundaryLhsPipeline(), mField, "U", Sev::verbose);
425
427}
428//! [Push operators to pipeline]
429
430//! [Solve]
431struct SetUpSchur {
432 static boost::shared_ptr<SetUpSchur>
435
436protected:
437 SetUpSchur() = default;
438};
439
443 auto pip = mField.getInterface<PipelineManager>();
444 auto solver = pip->createKSP();
445 CHKERR KSPSetFromOptions(solver);
446
447 auto dm = simple->getDM();
448 auto D = createDMVector(dm);
449 auto F = vectorDuplicate(D);
450
451 auto set_essential_bc = [&]() {
453 // This is low level pushing finite elements (pipelines) to solver
454 auto ksp_ctx_ptr = getDMKspCtx(simple->getDM());
455
456 auto pre_proc_rhs = boost::make_shared<FEMethod>();
457 auto post_proc_rhs = boost::make_shared<FEMethod>();
458 auto post_proc_lhs = boost::make_shared<FEMethod>();
459
460 auto get_pre_proc_hook = [&]() {
462 {});
463 };
464 pre_proc_rhs->preProcessHook = get_pre_proc_hook();
465
466 auto get_post_proc_hook_rhs = [&]() {
468 post_proc_rhs, 1.);
469 };
470 auto get_post_proc_hook_lhs = [&]() {
472 post_proc_lhs, 1.);
473 };
474 post_proc_rhs->postProcessHook = get_post_proc_hook_rhs();
475 post_proc_lhs->postProcessHook = get_post_proc_hook_lhs();
476
477 ksp_ctx_ptr->getPreProcComputeRhs().push_front(pre_proc_rhs);
478 ksp_ctx_ptr->getPostProcComputeRhs().push_back(post_proc_rhs);
479 ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_lhs);
481 };
482
483 auto setup_and_solve = [&]() {
485 BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
486 MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp";
487 CHKERR KSPSetUp(solver);
488 MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp <= Done";
489 MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve";
490 CHKERR KSPSolve(solver, F, D);
491 MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve <= Done";
493 };
494
495 MOFEM_LOG_CHANNEL("TIMER");
496 MOFEM_LOG_TAG("TIMER", "timer");
497
498 CHKERR set_essential_bc();
499
500 if (A == AssemblyType::SCHUR) {
501 auto schur_ptr = SetUpSchur::createSetUpSchur(mField);
502 CHKERR schur_ptr->setUp(solver);
503 CHKERR setup_and_solve();
504 } else {
505 CHKERR setup_and_solve();
506 }
507
508 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
509 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
510 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
511
512 if (doEvalField) {
513 if constexpr (SPACE_DIM == 3) {
514 CHKERR mField.getInterface<FieldEvaluatorInterface>()->evalFEAtThePoint3D(
515 fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
516 simple->getDomainFEName(), fieldEvalData, mField.get_comm_rank(),
517 mField.get_comm_rank(), nullptr, MF_EXIST, QUIET);
518 } else {
519 CHKERR mField.getInterface<FieldEvaluatorInterface>()->evalFEAtThePoint2D(
520 fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
521 simple->getDomainFEName(), fieldEvalData, mField.get_comm_rank(),
522 mField.get_comm_rank(), nullptr, MF_EXIST, QUIET);
523 }
524
525 if (vectorFieldPtr->size1()) {
526 auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vectorFieldPtr);
527 if constexpr (SPACE_DIM == 2)
528 MOFEM_LOG("FieldEvaluator", Sev::inform)
529 << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1);
530 else
531 MOFEM_LOG("FieldEvaluator", Sev::inform)
532 << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1)
533 << " U_Z: " << t_disp(2);
534 }
535
537 }
539}
540//! [Solve]
541
542//! [Postprocess results]
546 auto pip = mField.getInterface<PipelineManager>();
547 auto det_ptr = boost::make_shared<VectorDouble>();
548 auto jac_ptr = boost::make_shared<MatrixDouble>();
549 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
550 pip->getDomainRhsFE().reset();
551 pip->getDomainLhsFE().reset();
552 pip->getBoundaryRhsFE().reset();
553 pip->getBoundaryLhsFE().reset();
554
555 auto post_proc_mesh = boost::make_shared<moab::Core>();
556 auto post_proc_begin = boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(
557 mField, post_proc_mesh);
558 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
559 mField, post_proc_mesh);
560
561 auto calculate_stress_ops = [&](auto &pip) {
563 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
564 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
565 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
567 "U", mat_grad_ptr));
568 pip.push_back(
569 new OpSymmetrizeTensor<SPACE_DIM>("U", mat_grad_ptr, mat_strain_ptr));
570
571 auto mat_D_ptr = boost::make_shared<MatrixDouble>();
572 CHKERR addMatBlockOps(pip, "U", "MAT_ELASTIC", mat_D_ptr, Sev::verbose);
574 "U", mat_strain_ptr, mat_stress_ptr, mat_D_ptr));
575 auto u_ptr = boost::make_shared<MatrixDouble>();
576 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
577 auto x_ptr = boost::make_shared<MatrixDouble>();
578 pip.push_back(
579 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
580 return boost::make_tuple(u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr);
581 };
582
583 auto post_proc_domain = [&](auto post_proc_mesh) {
584 auto post_proc_fe =
585 boost::make_shared<PostProcEleDomain>(mField, post_proc_mesh);
587
588 auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr] =
589 calculate_stress_ops(post_proc_fe->getOpPtrVector());
590
591 post_proc_fe->getOpPtrVector().push_back(
592
593 new OpPPMap(
594
595 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
596
597 {},
598
599 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
600
601 {},
602
603 {{"STRAIN", mat_strain_ptr}, {"STRESS", mat_stress_ptr}}
604
605 )
606
607 );
608 return post_proc_fe;
609 };
610
611 auto post_proc_boundary = [&](auto post_proc_mesh) {
612 auto post_proc_fe =
613 boost::make_shared<PostProcEleBdy>(mField, post_proc_mesh);
615 post_proc_fe->getOpPtrVector(), {}, "GEOMETRY");
616 auto op_loop_side =
617 new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
618 // push ops to side element, through op_loop_side operator
619 auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr] =
620 calculate_stress_ops(op_loop_side->getOpPtrVector());
621 post_proc_fe->getOpPtrVector().push_back(op_loop_side);
622 auto mat_traction_ptr = boost::make_shared<MatrixDouble>();
623 post_proc_fe->getOpPtrVector().push_back(
624 new ElasticExample::OpCalculateTraction(mat_stress_ptr,
625 mat_traction_ptr));
626
628
629 post_proc_fe->getOpPtrVector().push_back(
630
631 new OpPPMap(
632
633 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
634
635 {},
636
637 {{"U", u_ptr}, {"GEOMETRY", x_ptr}, {"T", mat_traction_ptr}},
638
639 {},
640
641 {{"STRAIN", mat_strain_ptr}, {"STRESS", mat_stress_ptr}}
642
643 )
644
645 );
646 return post_proc_fe;
647 };
648
649 PetscBool post_proc_skin_only = PETSC_FALSE;
650 if (SPACE_DIM == 3) {
651 post_proc_skin_only = PETSC_TRUE;
652 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_skin_only",
653 &post_proc_skin_only, PETSC_NULL);
654 }
655 if (post_proc_skin_only == PETSC_FALSE) {
656 pip->getDomainRhsFE() = post_proc_domain(post_proc_mesh);
657 }
658 pip->getBoundaryRhsFE() = post_proc_boundary(post_proc_mesh);
659
661 post_proc_begin->getFEMethod());
662 CHKERR pip->loopFiniteElements();
664 post_proc_end->getFEMethod());
665
666 CHKERR post_proc_end->writeFile("out_elastic.h5m");
668}
669//! [Postprocessing results]
670
671//! [Check]
673 MOFEM_LOG_CHANNEL("WORLD");
675 auto pip = mField.getInterface<PipelineManager>();
677 pip->getDomainRhsFE().reset();
678 pip->getDomainLhsFE().reset();
679 pip->getBoundaryRhsFE().reset();
680 pip->getBoundaryLhsFE().reset();
681
682 auto integration_rule = [](int, int, int p_data) { return 2 * p_data + 1; };
683 CHKERR pip->setDomainRhsIntegrationRule(integration_rule);
684 CHKERR pip->setBoundaryRhsIntegrationRule(integration_rule);
685
687 pip->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
689 pip->getOpBoundaryRhsPipeline(), {}, "GEOMETRY");
690
691 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
692 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
693 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
694
695 pip->getOpDomainRhsPipeline().push_back(
697 mat_grad_ptr));
698 pip->getOpDomainRhsPipeline().push_back(
699 new OpSymmetrizeTensor<SPACE_DIM>("U", mat_grad_ptr, mat_strain_ptr));
700
701 auto mat_D_ptr = boost::make_shared<MatrixDouble>();
702 CHKERR addMatBlockOps(pip->getOpDomainRhsPipeline(), "U", "MAT_ELASTIC",
703 mat_D_ptr, Sev::verbose);
704 pip->getOpDomainRhsPipeline().push_back(
706 "U", mat_strain_ptr, mat_stress_ptr, mat_D_ptr));
707
708 pip->getOpDomainRhsPipeline().push_back(
709 new OpInternalForce("U", mat_stress_ptr));
710
711 pip->getOpBoundaryRhsPipeline().push_back(
713 "U", mat_strain_ptr, mat_stress_ptr, mat_D_ptr));
715 pip->getOpDomainRhsPipeline(), mField, "U", Sev::verbose);
717 pip->getOpBoundaryRhsPipeline(), mField, "U", 1, Sev::verbose);
718
719 auto dm = simple->getDM();
720 auto res = createDMVector(dm);
721 pip->getDomainRhsFE()->f = res;
722 pip->getBoundaryRhsFE()->f = res;
723
724 CHKERR VecZeroEntries(res);
725
726 CHKERR mField.getInterface<FieldBlas>()->fieldScale(-1, "U");
727 CHKERR pip->loopFiniteElements();
728 CHKERR mField.getInterface<FieldBlas>()->fieldScale(-1, "U");
729
730 CHKERR VecGhostUpdateBegin(res, ADD_VALUES, SCATTER_REVERSE);
731 CHKERR VecGhostUpdateEnd(res, ADD_VALUES, SCATTER_REVERSE);
732 CHKERR VecAssemblyBegin(res);
733 CHKERR VecAssemblyEnd(res);
734
735 auto zero_residual_at_constrains = [&]() {
737 auto fe_post_proc_ptr = boost::make_shared<FEMethod>();
738 auto get_post_proc_hook_rhs = [this, fe_post_proc_ptr, res]() {
741 mField, fe_post_proc_ptr, res)();
743 mField, fe_post_proc_ptr, 0, res)();
745 };
746 fe_post_proc_ptr->postProcessHook = get_post_proc_hook_rhs;
747 CHKERR DMoFEMPostProcessFiniteElements(dm, fe_post_proc_ptr.get());
749 };
750
751 CHKERR zero_residual_at_constrains();
752
753 double nrm2;
754 CHKERR VecNorm(res, NORM_2, &nrm2);
755 MOFEM_LOG_CHANNEL("WORLD");
756 MOFEM_LOG_C("WORLD", Sev::inform, "residual = %3.4e\n", nrm2);
757
758 PetscBool test = PETSC_FALSE;
759 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test, PETSC_NULL);
760 if (test == PETSC_TRUE) {
761
762 auto post_proc_residual = [&](auto dm, auto f_res, auto out_name) {
764 auto post_proc_fe =
765 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
767 auto u_vec = boost::make_shared<MatrixDouble>();
768 post_proc_fe->getOpPtrVector().push_back(
769 new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_vec, f_res));
770 post_proc_fe->getOpPtrVector().push_back(
771
772 new OpPPMap(
773
774 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
775
776 {},
777
778 {{"RES", u_vec}},
779
780 {}, {})
781
782 );
783
784 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
785 post_proc_fe);
786 post_proc_fe->writeFile(out_name);
788 };
789
790 CHKERR post_proc_residual(simple->getDM(), res, "res.h5m");
791
792 constexpr double eps = 1e-8;
793 if (nrm2 > eps)
794 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
795 "Residual is not zero");
796 }
797
799}
800//! [Check]
801
802static char help[] = "...\n\n";
803
804int main(int argc, char *argv[]) {
805
806 // Initialisation of MoFEM/PETSc and MOAB data structures
807 const char param_file[] = "param_file.petsc";
809
810 auto core_log = logging::core::get();
811 core_log->add_sink(
813
814 core_log->add_sink(
816 LogManager::setLog("FieldEvaluator");
817 MOFEM_LOG_TAG("FieldEvaluator", "field_eval");
818
819 try {
820
821 //! [Register MoFEM discrete manager in PETSc]
822 DMType dm_name = "DMMOFEM";
823 CHKERR DMRegister_MoFEM(dm_name);
824 //! [Register MoFEM discrete manager in PETSc
825
826 //! [Create MoAB]
827 moab::Core mb_instance; ///< mesh database
828 moab::Interface &moab = mb_instance; ///< mesh database interface
829 //! [Create MoAB]
830
831 //! [Create MoFEM]
832 MoFEM::Core core(moab); ///< finite element database
833 MoFEM::Interface &m_field = core; ///< finite element database interface
834 //! [Create MoFEM]
835
836 //! [Example]
837 Example ex(m_field);
838 CHKERR ex.runProblem();
839 //! [Example]
840 }
842
844}
845
846struct SetUpSchurImpl : public SetUpSchur {
847
849 if (S) {
852 "Is expected that schur matrix is not allocated. This is "
853 "possible only is PC is set up twice");
854 }
855 }
856 virtual ~SetUpSchurImpl() { S.reset(); }
857
859
860private:
865
867
869
873};
874
877 auto pip = mField.getInterface<PipelineManager>();
878 PC pc;
879 CHKERR KSPGetPC(solver, &pc);
880 PetscBool is_pcfs = PETSC_FALSE;
881 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
882 if (is_pcfs) {
883 if (S) {
886 "Is expected that schur matrix is not allocated. This is "
887 "possible only is PC is set up twice");
888 }
892 CHKERR MatSetBlockSize(S, SPACE_DIM);
893 CHKERR MatSetOption(S, MAT_SYMMETRIC, PETSC_TRUE);
895 CHKERR setPC(pc);
896 } else {
897 pip->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
898 pip->getOpBoundaryLhsPipeline().push_back(
899 new OpSchurAssembleEnd<SCHUR_DSYSV>({}, {}, {}, {}, {}));
900 pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
901 pip->getOpDomainLhsPipeline().push_back(
902 new OpSchurAssembleEnd<SCHUR_DSYSV>({}, {}, {}, {}, {}));
903 }
905}
906
910 CHKERR mField.get_moab().get_entities_by_dimension(simple->getMeshset(),
912 CHKERR mField.get_moab().get_entities_by_handle(simple->getMeshset(),
913 subEnts);
914 subEnts = subtract(subEnts, volEnts);
916};
917
921 subDM = createDM(mField.get_comm(), "DMMOFEM");
922 CHKERR DMMoFEMCreateSubDM(subDM, simple->getDM(), "SUB");
924 CHKERR DMMoFEMAddElement(subDM, simple->getDomainFEName());
925 auto sub_ents_ptr = boost::make_shared<Range>(subEnts);
926 CHKERR DMMoFEMAddSubFieldRow(subDM, "U", sub_ents_ptr);
927 CHKERR DMMoFEMAddSubFieldCol(subDM, "U", sub_ents_ptr);
928 CHKERR DMSetUp(subDM);
930}
931
934 auto pip = mField.getInterface<PipelineManager>();
935
936 // Boundary
937 auto dm_is = getDMSubData(subDM)->getSmartRowIs();
938 auto ao_up = createAOMappingIS(dm_is, PETSC_NULL);
939 pip->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
940 pip->getOpBoundaryLhsPipeline().push_back(new OpSchurAssembleEnd<SCHUR_DSYSV>(
941 {"U"}, {boost::make_shared<Range>(volEnts)}, {ao_up}, {S}, {true}));
942 // Domain
943 pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
944 pip->getOpDomainLhsPipeline().push_back(new OpSchurAssembleEnd<SCHUR_DSYSV>(
945 {"U"}, {boost::make_shared<Range>(volEnts)}, {ao_up}, {S}, {true}));
946
947 auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
948 auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
949
950 pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
952 if (S) {
953 CHKERR MatZeroEntries(S);
954 }
955 MOFEM_LOG("TIMER", Sev::inform) << "Lhs Assemble Begin";
957 };
958
959 post_proc_schur_lhs_ptr->postProcessHook = [this, post_proc_schur_lhs_ptr,
960 ao_up]() {
962 CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
963 CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
965 mField, post_proc_schur_lhs_ptr, 1, S, ao_up)();
966 MOFEM_LOG("TIMER", Sev::inform) << "Lhs Assemble End";
968 };
969
971 auto ksp_ctx_ptr = getDMKspCtx(simple->getDM());
972
973 ksp_ctx_ptr->getPreProcSetOperators().push_front(pre_proc_schur_lhs_ptr);
974 ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_schur_lhs_ptr);
975
977}
978
982 SmartPetscObj<IS> vol_is;
983 mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
984 simple->getProblemName(), ROW, "U", 0, SPACE_DIM, vol_is, &volEnts);
985 CHKERR PCFieldSplitSetIS(pc, NULL, vol_is);
986 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
988}
989
990boost::shared_ptr<SetUpSchur>
992 return boost::shared_ptr<SetUpSchur>(new SetUpSchurImpl(m_field));
993}
Calculate traction for linear problem.
std::string param_file
Implementation of elastic spring bc.
#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
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.
@ 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
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
PostProcEleByDim< SPACE_DIM >::SideEle SideEle
PostProcEleByDim< SPACE_DIM >::PostProcEleBdy PostProcEleBdy
double bulk_modulus_K
double shear_modulus_G
static char help[]
[Check]
Definition: elastic.cpp:802
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesSymTensor< BASE_DIM, SPACE_DIM, SPACE_DIM > OpInternalForce
Definition: elastic.cpp:36
constexpr int SPACE_DIM
Definition: elastic.cpp:17
constexpr double poisson_ratio
Definition: elastic.cpp:72
constexpr int BASE_DIM
[Define dimension]
Definition: elastic.cpp:16
constexpr double shear_modulus_G
Definition: elastic.cpp:74
constexpr IntegrationType I
Definition: elastic.cpp:23
constexpr double bulk_modulus_K
Definition: elastic.cpp:73
constexpr AssemblyType A
Definition: elastic.cpp:19
PetscBool is_plain_strain
Definition: elastic.cpp:76
constexpr double young_modulus
Definition: elastic.cpp:71
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:165
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
constexpr auto size_symm
Definition: plastic.cpp:40
constexpr auto field_name
static constexpr int approx_order
[Example]
Definition: plastic.cpp:217
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
PetscBool doEvalField
Definition: elastic.cpp:87
FieldApproximationBase base
Definition: plot_base.cpp:68
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
Example(MoFEM::Interface &m_field)
Definition: elastic.cpp:80
MoFEMErrorCode runProblem()
MoFEM::Interface & mField
Definition: plastic.cpp:224
std::array< double, SPACE_DIM > fieldEvalCoords
Definition: elastic.cpp:88
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:106
boost::shared_ptr< MatrixDouble > vectorFieldPtr
Definition: elastic.cpp:90
boost::shared_ptr< FieldEvaluatorInterface::SetPtsData > fieldEvalData
Definition: elastic.cpp:89
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: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
Class (Function) to calculate residual side diagonal.
Definition: Essential.hpp:63
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:692
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:1411
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
SmartPetscObj< DM > subDM
field split sub dm
Definition: plastic.cpp:1233
SmartPetscObj< Mat > S
Definition: plastic.cpp:1230
SetUpSchurImpl(MoFEM::Interface &m_field)
Definition: elastic.cpp:848
MoFEMErrorCode setUp(TS solver)
Definition: plastic.cpp:1238
MoFEMErrorCode setEntities()
virtual ~SetUpSchurImpl()
Definition: elastic.cpp:856
MoFEMErrorCode setPC(PC pc)
MoFEMErrorCode setUpSubDM()
Definition: elastic.cpp:918
MoFEMErrorCode setOperator()
MoFEM::Interface & mField
Definition: plastic.cpp:1232