v0.14.0
Loading...
Searching...
No Matches
EshelbianPlasticity.cpp
Go to the documentation of this file.
1/**
2 * \file EshelbianPlasticity.cpp
3 * \example EshelbianPlasticity.cpp
4 *
5 * \brief Eshelbian plasticity implementation
6 */
7
8#include <MoFEM.hpp>
9using namespace MoFEM;
10
12
14#include <boost/math/constants/constants.hpp>
15
16#include <cholesky.hpp>
17#ifdef PYTHON_SFD
18#include <boost/python.hpp>
19#include <boost/python/def.hpp>
20namespace bp = boost::python;
21
22#pragma message "With PYTHON_SFD"
23
24#else
25
26#pragma message "Without PYTHON_SFD"
27
28#endif
29
30#include <EshelbianContact.hpp>
31
32namespace EshelbianPlasticity {
33
37
38double EshelbianCore::exponentBase = exp(1);
39boost::function<double(const double)> EshelbianCore::f = EshelbianCore::f_log;
40boost::function<double(const double)> EshelbianCore::d_f =
42boost::function<double(const double)> EshelbianCore::dd_f =
44boost::function<double(const double)> EshelbianCore::inv_f =
46boost::function<double(const double)> EshelbianCore::inv_d_f =
48boost::function<double(const double)> EshelbianCore::inv_dd_f =
50
52EshelbianCore::query_interface(boost::typeindex::type_index type_index,
53 UnknownInterface **iface) const {
54 *iface = const_cast<EshelbianCore *>(this);
55 return 0;
56}
57
60
61 if (evalRhs)
62 CHKERR evaluateRhs(data);
63
64 if (evalLhs)
65 CHKERR evaluateLhs(data);
66
68}
69
71 : mField(m_field), piolaStress("P"), eshelbyStress("S"),
72 spatialL2Disp("wL2"), spatialH1Disp("wH1"), materialL2Disp("W"),
73 StretchTensor("u"), rotAxis("omega"), materialGradient("G"),
74 tauField("TAU"), lambdaField("LAMBDA"), bubbleField("BUBBLE"),
75 elementVolumeName("EP"), naturalBcElement("NATURAL_BC"),
76 essentialBcElement("ESSENTIAL_BC"), skinElement("SKIN_ELEMENT"),
77 contactElement("CONTACT_ELEMENT") {
78
79 ierr = getOptions();
80 CHKERRABORT(PETSC_COMM_WORLD, ierr);
81}
82
84
87 const char *list_rots[] = {"small", "moderate", "large"};
88 PetscInt choice_rot = EshelbianCore::rotSelector;
89 PetscInt choice_grad = EshelbianCore::gradApperoximator;
90
91 const char *list_stretches[] = {"linear", "log"};
92 PetscInt choice_stretch = StretchSelector::LOG;
93
94 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Eshelbian plasticity",
95 "none");
96
97 spaceOrder = 2;
98 CHKERR PetscOptionsInt("-space_order", "approximation oder for space", "",
99 spaceOrder, &spaceOrder, PETSC_NULL);
100 spaceH1Order = -1;
101 CHKERR PetscOptionsInt("-space_h1_order", "approximation oder for space", "",
102 spaceH1Order, &spaceH1Order, PETSC_NULL);
103
104 materialOrder = 1;
105 CHKERR PetscOptionsInt("-material_order", "approximation oder for material",
106 "", materialOrder, &materialOrder, PETSC_NULL);
107
108 alphaU = 0;
109 CHKERR PetscOptionsScalar("-viscosity_alpha_u", "viscosity", "", alphaU,
110 &alphaU, PETSC_NULL);
111
112 alphaW = 0;
113 CHKERR PetscOptionsScalar("-viscosity_alpha_w", "viscosity", "", alphaW,
114 &alphaW, PETSC_NULL);
115
116 alphaRho = 0;
117 CHKERR PetscOptionsScalar("-density_alpha_rho", "density", "", alphaRho,
118 &alphaRho, PETSC_NULL);
119
120 precEps = 0;
121 CHKERR PetscOptionsScalar("-preconditioner_eps", "preconditioner_eps", "",
122 precEps, &precEps, PETSC_NULL);
123
124 precEpsOmega = 0;
125 CHKERR PetscOptionsScalar("-preconditioner_eps_omega", "preconditioner_eps",
126 "", precEpsOmega, &precEpsOmega, PETSC_NULL);
127 precEpsW = 0;
128 CHKERR PetscOptionsScalar("-preconditioner_eps_w", "preconditioner_eps", "",
129 precEpsW, &precEpsW, PETSC_NULL);
130
131 CHKERR PetscOptionsEList("-rotations", "rotations", "", list_rots,
132 LARGE_ROT + 1, list_rots[choice_rot], &choice_rot,
133 PETSC_NULL);
134 CHKERR PetscOptionsEList("-grad", "gradient of defomation approximator", "",
135 list_rots, LARGE_ROT + 1, list_rots[choice_grad],
136 &choice_grad, PETSC_NULL);
137
138 CHKERR PetscOptionsScalar("-exponent_base", "exponent_base", "", exponentBase,
139 &EshelbianCore::exponentBase, PETSC_NULL);
140 CHKERR PetscOptionsEList(
141 "-stretches", "stretches", "", list_stretches, StretchSelector::LOG + 1,
142 list_stretches[choice_stretch], &choice_stretch, PETSC_NULL);
143
144 ierr = PetscOptionsEnd();
145 CHKERRG(ierr);
146
147 EshelbianCore::rotSelector = static_cast<RotSelector>(choice_rot);
148 EshelbianCore::gradApperoximator = static_cast<RotSelector>(choice_grad);
149 EshelbianCore::stretchSelector = static_cast<StretchSelector>(choice_stretch);
150
159 break;
167 break;
168 default:
169 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown stretch");
170 break;
171 };
172
174 precEpsW += precEps;
175
176 MOFEM_LOG("EP", Sev::inform) << "spaceOrder " << spaceOrder;
177 MOFEM_LOG("EP", Sev::inform) << "spaceH1Order " << spaceH1Order;
178 MOFEM_LOG("EP", Sev::inform) << "materialOrder " << materialOrder;
179 MOFEM_LOG("EP", Sev::inform) << "alphaU " << alphaU;
180 MOFEM_LOG("EP", Sev::inform) << "alphaW " << alphaW;
181 MOFEM_LOG("EP", Sev::inform) << "alphaRho " << alphaRho;
182 MOFEM_LOG("EP", Sev::inform) << "precEps " << precEps;
183 MOFEM_LOG("EP", Sev::inform) << "precEpsOmega " << precEpsOmega;
184 MOFEM_LOG("EP", Sev::inform) << "precEpsW " << precEpsW;
185 MOFEM_LOG("EP", Sev::inform)
186 << "Rotations " << list_rots[EshelbianCore::rotSelector];
187 MOFEM_LOG("EP", Sev::inform) << "Gradient of deformation "
189 if (exponentBase != exp(1))
190 MOFEM_LOG("EP", Sev::inform)
191 << "Base exponent " << EshelbianCore::exponentBase;
192 else
193 MOFEM_LOG("EP", Sev::inform) << "Base exponent e";
194 MOFEM_LOG("EP", Sev::inform) << "Stretch " << list_stretches[choice_stretch];
195
196 if (spaceH1Order == -1)
198
200}
201
204
205 Range tets;
206 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
207 Range tets_skin_part;
208 Skinner skin(&mField.get_moab());
209 CHKERR skin.find_skin(0, tets, false, tets_skin_part);
210 ParallelComm *pcomm =
211 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
212 Range tets_skin;
213 CHKERR pcomm->filter_pstatus(tets_skin_part,
214 PSTATUS_SHARED | PSTATUS_MULTISHARED,
215 PSTATUS_NOT, -1, &tets_skin);
216
218 for (auto &v : *bcSpatialDispVecPtr) {
219 tets_skin = subtract(tets_skin, v.faces);
220 }
222 for (auto &v : *bcSpatialRotationVecPtr) {
223 tets_skin = subtract(tets_skin, v.faces);
224 }
226 for (auto &v : *bcSpatialTraction) {
227 tets_skin = subtract(tets_skin, v.faces);
228 }
229
230 auto subtract_faces_where_displacements_are_applied =
231 [&](const std::string block_name) {
233 auto msh_mng = mField.getInterface<MeshsetsManager>();
234 auto reg_exp = std::regex((boost::format("%s(.*)") % block_name).str());
235 for (auto m : msh_mng->getCubitMeshsetPtr(reg_exp)) {
236 Range faces;
237 CHKERR m->getMeshsetIdEntitiesByDimension(mField.get_moab(),
238 SPACE_DIM - 1, faces, true);
239 tets_skin = subtract(tets_skin, faces);
240 MOFEM_LOG("EP", Sev::inform)
241 << "Subtracting " << m->getName() << " " << faces.size();
242 }
244 };
245 CHKERR subtract_faces_where_displacements_are_applied("CONTACT");
246
247 Range faces;
248 CHKERR mField.get_moab().get_adjacencies(tets, 2, true, faces,
249 moab::Interface::UNION);
250 Range faces_not_on_the_skin = subtract(faces, tets_skin);
251
252 auto add_hdiv_field = [&](const std::string field_name, const int order,
253 const int dim) {
256 MB_TAG_SPARSE, MF_ZERO);
258 CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
259 CHKERR mField.set_field_order(faces_not_on_the_skin, field_name, order);
262 };
263
264 auto add_l2_field = [this, meshset](const std::string field_name,
265 const int order, const int dim) {
268 MB_TAG_DENSE, MF_ZERO);
270 CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
272 };
273
274 auto add_h1_field = [this, meshset](const std::string field_name,
275 const int order, const int dim) {
278 MB_TAG_DENSE, MF_ZERO);
280 CHKERR mField.set_field_order(meshset, MBVERTEX, field_name, 1);
281 CHKERR mField.set_field_order(meshset, MBEDGE, field_name, order);
282 CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
283 CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
285 };
286
287 auto add_bubble_field = [this, meshset](const std::string field_name,
288 const int order, const int dim) {
291 MF_ZERO);
292 // Modify field
293 auto field_ptr = mField.get_field_structure(field_name);
294 auto field_order_table =
295 const_cast<Field *>(field_ptr)->getFieldOrderTable();
296 auto get_cgg_bubble_order_zero = [](int p) { return 0; };
297 auto get_cgg_bubble_order_tet = [](int p) {
299 };
300 field_order_table[MBVERTEX] = get_cgg_bubble_order_zero;
301 field_order_table[MBEDGE] = get_cgg_bubble_order_zero;
302 field_order_table[MBTRI] = get_cgg_bubble_order_zero;
303 field_order_table[MBTET] = get_cgg_bubble_order_tet;
305 CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
306 CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
308 };
309
310 // spatial fields
311 CHKERR add_hdiv_field(piolaStress, spaceOrder, 3);
312 CHKERR add_bubble_field(bubbleField, spaceOrder, 1);
313 CHKERR add_l2_field(spatialL2Disp, spaceOrder - 1, 3);
314 CHKERR add_l2_field(rotAxis, spaceOrder - 1, 3);
315 CHKERR add_l2_field(StretchTensor, spaceOrder, 6);
316
317 // material fields
318 // CHKERR add_hdiv_field(eshelbyStress, materialOrder, 3);
319 // CHKERR add_l2_field(materialGradient, materialOrder - 1, 9);
320 // CHKERR add_l2_field(materialL2Disp, materialOrder - 1, 3);
321 // CHKERR add_l2_field(tauField, materialOrder - 1, 1);
322 // CHKERR add_l2_field(lambdaField, materialOrder - 1, 1);
323
324 // Add history filedes
325 // CHKERR add_l2_field(materialGradient + "0", materialOrder - 1, 9);
326
327 // spatial displacement
328 CHKERR add_h1_field(spatialH1Disp, spaceH1Order, 3);
329
331
333}
334
338
339 // set finite element fields
340 auto add_field_to_fe = [this](const std::string fe,
341 const std::string field_name) {
347 };
348
353
354 CHKERR add_field_to_fe(elementVolumeName, piolaStress);
355 CHKERR add_field_to_fe(elementVolumeName, bubbleField);
356 // CHKERR add_field_to_fe(elementVolumeName, eshelbyStress);
357 CHKERR add_field_to_fe(elementVolumeName, StretchTensor);
358 CHKERR add_field_to_fe(elementVolumeName, rotAxis);
359 CHKERR add_field_to_fe(elementVolumeName, spatialL2Disp);
360 CHKERR add_field_to_fe(elementVolumeName, StretchTensor);
361 CHKERR add_field_to_fe(elementVolumeName, spatialH1Disp);
362 // CHKERR add_field_to_fe(elementVolumeName, materialGradient);
363 // CHKERR mField.modify_finite_element_add_field_data(elementVolumeName,
364 // materialGradient +
365 // "0");
366 }
367
368 // build finite elements data structures
370
372}
373
377
378 auto bc_elements_add_to_range = [&](const std::string block_name, Range &r) {
380
381 auto mesh_mng = mField.getInterface<MeshsetsManager>();
382 auto bcs = mesh_mng->getCubitMeshsetPtr(
383
384 std::regex((boost::format("%s(.*)") % block_name).str())
385
386 );
387
388 for (auto bc : bcs) {
389 Range faces;
390 CHKERR bc->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
391 true);
392 r.merge(faces);
393 }
394
396 };
397
398 // set finite element fields
399 auto add_field_to_fe = [this](const std::string fe,
400 const std::string field_name) {
406 };
407
408 Range natural_bc_elements;
410 for (auto &v : *bcSpatialDispVecPtr) {
411 natural_bc_elements.merge(v.faces);
412 }
413 }
415 for (auto &v : *bcSpatialRotationVecPtr) {
416 natural_bc_elements.merge(v.faces);
417 }
418 }
419 Range essentail_bc_elements;
420 if (bcSpatialTraction) {
421 for (auto &v : *bcSpatialTraction) {
422 essentail_bc_elements.merge(v.faces);
423 }
424 }
425
427 CHKERR mField.add_ents_to_finite_element_by_type(natural_bc_elements, MBTRI,
429 CHKERR add_field_to_fe(naturalBcElement, piolaStress);
430 // CHKERR add_field_to_fe(naturalBcElement, eshelbyStress);
432
434 CHKERR mField.add_ents_to_finite_element_by_type(essentail_bc_elements, MBTRI,
436 CHKERR add_field_to_fe(essentialBcElement, piolaStress);
437 // CHKERR add_field_to_fe(essentialBcElement, eshelbyStress);
439
440 auto get_skin = [&]() {
441 Range body_ents;
442 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
443 Skinner skin(&mField.get_moab());
444 Range skin_ents;
445 CHKERR skin.find_skin(0, body_ents, false, skin_ents);
446 return skin_ents;
447 };
448
449 auto filter_true_skin = [&](auto &&skin) {
450 Range boundary_ents;
451 ParallelComm *pcomm =
452 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
453 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
454 PSTATUS_NOT, -1, &boundary_ents);
455 return boundary_ents;
456 };
457
458 auto skin = filter_true_skin(get_skin());
459
462 CHKERR add_field_to_fe(skinElement, piolaStress);
463 CHKERR add_field_to_fe(skinElement, spatialH1Disp);
464 // CHKERR add_field_to_fe(skinElement, eshelbyStress);
466
467 Range contact_range;
468 CHKERR bc_elements_add_to_range("CONTACT", contact_range);
469 MOFEM_LOG("EP", Sev::inform) << "Contact elements " << contact_range.size();
470
474 CHKERR add_field_to_fe(contactElement, piolaStress);
476
478}
479
482
483 // find adjacencies between finite elements and dofs
485
486 // Create coupled problem
487 dM = createDM(mField.get_comm(), "DMMOFEM");
488 CHKERR DMMoFEMCreateMoFEM(dM, &mField, "ESHELBY_PLASTICITY", bit,
489 BitRefLevel().set());
491 CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
497 mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
498 CHKERR DMSetUp(dM);
499 mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
500
501 auto remove_dofs_on_essential_spatial_stress_boundary =
502 [&](const std::string prb_name) {
504 for (int d : {0, 1, 2})
506 prb_name, piolaStress, (*bcSpatialFreeTraction)[d], d, d, 0,
509 };
510 CHKERR remove_dofs_on_essential_spatial_stress_boundary("ESHELBY_PLASTICITY");
511
512 // Create elastic sub-problem
513 dmElastic = createDM(mField.get_comm(), "DMMOFEM");
514 CHKERR DMMoFEMCreateSubDM(dmElastic, dM, "ELASTIC_PROBLEM");
527 CHKERR DMSetUp(dmElastic);
528
529 {
530 PetscSection section;
531 CHKERR mField.getInterface<ISManager>()->sectionCreate("ELASTIC_PROBLEM",
532 &section);
533 CHKERR DMSetSection(dmElastic, section);
534 CHKERR DMSetGlobalSection(dmElastic, section);
535 CHKERR PetscSectionDestroy(&section);
536 }
537
538 dmPrjSpatial = createDM(mField.get_comm(), "DMMOFEM");
539 CHKERR DMMoFEMCreateSubDM(dmPrjSpatial, dM, "PROJECT_SPATIAL");
545 CHKERR DMSetUp(dmPrjSpatial);
546
548 ->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
549 "PROJECT_SPATIAL", spatialH1Disp, true, false);
550
552}
553
554BcDisp::BcDisp(std::string name, std::vector<double> &attr, Range &faces)
555 : blockName(name), faces(faces) {
556 vals.resize(3, false);
557 flags.resize(3, false);
558 for (int ii = 0; ii != 3; ++ii) {
559 vals[ii] = attr[ii];
560 flags[ii] = static_cast<int>(attr[ii + 3]);
561 }
562
563 MOFEM_LOG("EP", Sev::inform) << "Add BCDisp " << name;
564 MOFEM_LOG("EP", Sev::inform)
565 << "Add BCDisp vals " << vals[0] << " " << vals[1] << " " << vals[2];
566 MOFEM_LOG("EP", Sev::inform)
567 << "Add BCDisp flags " << flags[0] << " " << flags[1] << " " << flags[2];
568 MOFEM_LOG("EP", Sev::inform) << "Add BCDisp nb. of faces " << faces.size();
569}
570
571BcRot::BcRot(std::string name, std::vector<double> &attr, Range &faces)
572 : blockName(name), faces(faces) {
573 vals.resize(3, false);
574 for (int ii = 0; ii != 3; ++ii) {
575 vals[ii] = attr[ii];
576 }
577 theta = attr[3];
578}
579
580TractionBc::TractionBc(std::string name, std::vector<double> &attr,
581 Range &faces)
582 : blockName(name), faces(faces) {
583 vals.resize(3, false);
584 flags.resize(3, false);
585 for (int ii = 0; ii != 3; ++ii) {
586 vals[ii] = attr[ii];
587 flags[ii] = static_cast<int>(attr[ii + 3]);
588 }
589
590 MOFEM_LOG("EP", Sev::inform) << "Add BCForce " << name;
591 MOFEM_LOG("EP", Sev::inform)
592 << "Add BCForce vals " << vals[0] << " " << vals[1] << " " << vals[2];
593 MOFEM_LOG("EP", Sev::inform)
594 << "Add BCForce flags " << flags[0] << " " << flags[1] << " " << flags[2];
595 MOFEM_LOG("EP", Sev::inform) << "Add BCForce nb. of faces " << faces.size();
596}
597
600 boost::shared_ptr<TractionFreeBc> &bc_ptr,
601 const std::string contact_set_name) {
603
604 // get skin from all tets
605 Range tets;
606 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
607 Range tets_skin_part;
608 Skinner skin(&mField.get_moab());
609 CHKERR skin.find_skin(0, tets, false, tets_skin_part);
610 ParallelComm *pcomm =
611 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
612 Range tets_skin;
613 CHKERR pcomm->filter_pstatus(tets_skin_part,
614 PSTATUS_SHARED | PSTATUS_MULTISHARED,
615 PSTATUS_NOT, -1, &tets_skin);
616
617 bc_ptr->resize(3);
618 for (int dd = 0; dd != 3; ++dd)
619 (*bc_ptr)[dd] = tets_skin;
620
622 for (auto &v : *bcSpatialDispVecPtr) {
623 if (v.flags[0])
624 (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
625 if (v.flags[1])
626 (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
627 if (v.flags[2])
628 (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
629 }
630
632 for (auto &v : *bcSpatialRotationVecPtr) {
633 (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
634 (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
635 (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
636 }
637
638 // remove contact
640 std::regex((boost::format("%s(.*)") % contact_set_name).str()))) {
641 Range faces;
642 CHKERR m->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
643 true);
644 (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
645 (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
646 (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
647 }
648
649 // for (int dd = 0; dd != 3; ++dd) {
650 // EntityHandle meshset;
651 // CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
652 // CHKERR mField.get_moab().add_entities(meshset, (*bc_ptr)[dd]);
653 // std::string file_name = disp_block_set_name +
654 // "_traction_free_bc_" + boost::lexical_cast<std::string>(dd) + ".vtk";
655 // CHKERR mField.get_moab().write_file(file_name.c_str(), " VTK ", "",
656 // &meshset, 1);
657 // CHKERR mField.get_moab().delete_entities(&meshset, 1);
658 // }
659
661}
662
663/**
664 * @brief Set integration rule on element
665 * \param order on row
666 * \param order on column
667 * \param order on data
668 *
669 * Use maximal oder on data in order to determine integration rule
670 *
671 */
672struct VolRule {
673 int operator()(int p_row, int p_col, int p_data) const {
675 return p_data + p_data + (p_data - 1);
676 else
677 return 2 * p_data;
678 }
679};
680
681struct FaceRule {
682 int operator()(int p_row, int p_col, int p_data) const { return 2 * p_data; }
683};
684
686
689
690 MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
691 BaseFunctionUnknownInterface **iface) const {
692 *iface = const_cast<CGGUserPolynomialBase *>(this);
693 return 0;
694 }
695
697 boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
699
701
702 int nb_gauss_pts = pts.size2();
703 if (!nb_gauss_pts) {
705 }
706
707 if (pts.size1() < 3) {
708 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
709 "Wrong dimension of pts, should be at least 3 rows with "
710 "coordinates");
711 }
712
713 switch (cTx->sPace) {
714 case HDIV:
716 break;
717 default:
718 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
719 }
720
722 }
723
724private:
726
728
731
732 // This should be used only in case USER_BASE is selected
733 if (cTx->bAse != USER_BASE) {
734 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
735 "Wrong base, should be USER_BASE");
736 }
737 // get access to data structures on element
738 EntitiesFieldData &data = cTx->dAta;
739 // Get approximation order on element. Note that bubble functions are only
740 // on tetrahedron.
741 const int order = data.dataOnEntities[MBTET][0].getOrder();
742 /// number of integration points
743 const int nb_gauss_pts = pts.size2();
744 // number of base functions
745
746 // calculate shape functions, i.e. barycentric coordinates
747 shapeFun.resize(nb_gauss_pts, 4, false);
748 CHKERR ShapeMBTET(&*shapeFun.data().begin(), &pts(0, 0), &pts(1, 0),
749 &pts(2, 0), nb_gauss_pts);
750 // direvatives of shape functions
751 double diff_shape_fun[12];
752 CHKERR ShapeDiffMBTET(diff_shape_fun);
753
754 const int nb_base_functions = NBVOLUMETET_CCG_BUBBLE(order);
755 // get base functions and set size
756 MatrixDouble &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
757 phi.resize(nb_gauss_pts, 9 * nb_base_functions, false);
758 // finally calculate base functions
760 &phi(0, 0), &phi(0, 1), &phi(0, 2),
761
762 &phi(0, 3), &phi(0, 4), &phi(0, 5),
763
764 &phi(0, 6), &phi(0, 7), &phi(0, 8));
765 CHKERR CGG_BubbleBase_MBTET(order, &shapeFun(0, 0), diff_shape_fun, t_phi,
766 nb_gauss_pts);
767
769 }
770};
771
773 const int tag, const bool do_rhs, const bool do_lhs,
774 boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe) {
776 fe = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
777
778 fe->getUserPolynomialBase() =
779 boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
781 {HDIV, H1, L2});
782
783 // set integration rule
784 fe->getRuleHook = VolRule();
785
786 if (!dataAtPts) {
787 dataAtPts =
788 boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
789 dataAtPts->physicsPtr = physicalEquations;
790 }
791
792 // calculate fields values
793 fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
794 piolaStress, dataAtPts->getApproxPAtPts()));
795 fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
796 bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
797 fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
798 piolaStress, dataAtPts->getDivPAtPts()));
799 fe->getOpPtrVector().push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
800 StretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
801 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
802 rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
803 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
804 spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
805
806 // velocities
807 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
808 spatialL2Disp, dataAtPts->getSmallWL2DotAtPts(), MBTET));
809 fe->getOpPtrVector().push_back(
811 StretchTensor, dataAtPts->getLogStretchDotTensorAtPts(), MBTET));
812 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
813 rotAxis, dataAtPts->getRotAxisDotAtPts(), MBTET));
814
815 // acceleration
816 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
817 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDotDot<3>(
818 spatialL2Disp, dataAtPts->getSmallWL2DotDotAtPts(), MBTET));
819 }
820
821 // H1 displacements
822 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
823 spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
824 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
825 spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
826
827 // calculate other derived quantities
828 fe->getOpPtrVector().push_back(
830
831 // evaluate integration points
832 fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
833 tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
834
836}
837
839 const int tag, const bool add_elastic, const bool add_material,
840 boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_rhs,
841 boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_lhs) {
843
844 auto time_scale = boost::make_shared<TimeScale>();
845
846 // Right hand side
847 CHKERR setBaseVolumeElementOps(tag, true, false, fe_rhs);
848
849 // elastic
850 if (add_elastic) {
851 fe_rhs->getOpPtrVector().push_back(
853 fe_rhs->getOpPtrVector().push_back(
855 fe_rhs->getOpPtrVector().push_back(
857 fe_rhs->getOpPtrVector().push_back(
859 fe_rhs->getOpPtrVector().push_back(
861 fe_rhs->getOpPtrVector().push_back(
863
864 // Body forces
865 using BodyNaturalBC =
867 Assembly<PETSC>::LinearForm<GAUSS>;
868 using OpBodyForce =
869 BodyNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 3>;
870 CHKERR BodyNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
871 fe_rhs->getOpPtrVector(), mField, "w", {time_scale}, "BODY_FORCE",
872 Sev::inform);
873 }
874
875 // Left hand side
876 CHKERR setBaseVolumeElementOps(tag, true, true, fe_lhs);
877
878 // elastic
879 if (add_elastic) {
880
881 const bool symmetric_system = gradApperoximator <= MODERATE_ROT;
882
883 fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_du(
885 fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dP(
887 fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dBubble(
889 if (!symmetric_system) {
890 fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
892 }
893
894 fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dP(
896 fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dw(
898
899 fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dP_domega(
900 piolaStress, rotAxis, dataAtPts, symmetric_system));
901 fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
902 bubbleField, rotAxis, dataAtPts, symmetric_system));
903
904 if (!symmetric_system) {
905 fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dP(
906 rotAxis, piolaStress, dataAtPts, false));
907 fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dBubble(
908 rotAxis, bubbleField, dataAtPts, false));
909 fe_lhs->getOpPtrVector().push_back(
911 }
912
913 // Stabilisation
914 if constexpr (A == AssemblyType::SCHUR) {
915 // Note that we assemble to AMat, however Assembly<SCHUR> assemble by
916 // default to PMat
919 fe_lhs->getOpPtrVector().push_back(
920 new OpMassStab(rotAxis, rotAxis, [this](double, double, double) {
921 return precEpsOmega;
922 }));
923 if (std::abs(alphaRho + alphaW) <
924 std::numeric_limits<double>::epsilon()) {
925
926 fe_lhs->getOpPtrVector().push_back(new OpMassStab(
928 [this](double, double, double) { return precEpsW; }));
929 }
930 }
931
932 if (add_material) {
933 }
934 }
935
937}
938
940 const bool add_elastic, const bool add_material,
941 boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
942 boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs) {
944
945 fe_rhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
946 fe_lhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
947
948 // set integration rule
949 fe_rhs->getRuleHook = FaceRule();
950 fe_lhs->getRuleHook = FaceRule();
951
953 {HDIV});
955 {HDIV});
956
957 if (add_elastic) {
958 fe_rhs->getOpPtrVector().push_back(
960 fe_rhs->getOpPtrVector().push_back(
962 }
963
965}
966
968 boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
969 boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs) {
971
972 /* The above code is written in C++ and it appears to be defining and using
973 various operations on boundary elements and side elements. */
976
977 fe_rhs = boost::make_shared<BoundaryEle>(mField);
978 fe_lhs = boost::make_shared<BoundaryEle>(mField);
979 fe_rhs->getRuleHook = FaceRule();
980 fe_lhs->getRuleHook = FaceRule();
981
983 {HDIV});
985 {HDIV});
986
987 auto adj_cache =
988 boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>();
989
990 auto get_op_side = [&](auto disp_ptr, auto col_ent_data_ptr) {
992 auto op_loop_side = new OpLoopSide<SideEle>(
993 mField, elementVolumeName, SPACE_DIM, Sev::noisy, adj_cache);
994 op_loop_side->getSideFEPtr()->getUserPolynomialBase() =
995 boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
997 op_loop_side->getOpPtrVector(), {HDIV, H1, L2});
998 op_loop_side->getOpPtrVector().push_back(
1000 if (col_ent_data_ptr != nullptr) {
1001 op_loop_side->getOpPtrVector().push_back(
1003 spatialL2Disp, col_ent_data_ptr));
1004 }
1005 return op_loop_side;
1006 };
1007
1008 auto add_ops_lhs = [&](auto &pip) {
1010
1011 auto col_ent_data_ptr = boost::make_shared<EntitiesFieldData::EntData>();
1012 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1013 pip.push_back(
1014 get_op_side(common_data_ptr->contactDispPtr(), col_ent_data_ptr));
1016 piolaStress, common_data_ptr->contactTractionPtr()));
1018 piolaStress, common_data_ptr, col_ent_data_ptr));
1020 piolaStress, common_data_ptr));
1022 };
1023
1024 auto add_ops_rhs = [&](auto &pip) {
1026 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1028 piolaStress, common_data_ptr->contactTractionPtr()));
1029 pip.push_back(get_op_side(common_data_ptr->contactDispPtr(), nullptr));
1031 piolaStress, common_data_ptr));
1033 };
1034
1035 CHKERR add_ops_lhs(fe_lhs->getOpPtrVector());
1036 CHKERR add_ops_rhs(fe_rhs->getOpPtrVector());
1037
1039}
1040
1045
1046 auto adj_cache =
1047 boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>();
1048
1049 auto get_op_contact_bc = [&]() {
1051 auto op_loop_side = new OpLoopSide<SideEle>(
1052 mField, contactElement, SPACE_DIM - 1, Sev::noisy, adj_cache);
1053 return op_loop_side;
1054 };
1055
1056 auto op_contact_bc = get_op_contact_bc();
1057 elasticFeLhs->getOpPtrVector().push_back(op_contact_bc);
1058 CHKERR setContactElementOps(contactRhs, op_contact_bc->getSideFEPtr());
1060}
1061
1064 boost::shared_ptr<FEMethod> null;
1065 boost::shared_ptr<FeTractionBc> spatial_traction_bc(
1067
1068 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
1069 CHKERR DMMoFEMTSSetI2Function(dm, DM_NO_ELEMENT, null, spatial_traction_bc,
1070 null);
1072 null);
1074 null);
1077 spatial_traction_bc);
1079 null);
1081 null);
1082
1083 } else {
1084 CHKERR DMMoFEMTSSetIFunction(dm, DM_NO_ELEMENT, null, spatial_traction_bc,
1085 null);
1087 null);
1089 null);
1092 spatial_traction_bc);
1094 null);
1096 null);
1097 }
1099}
1100
1103#include "impl/SetUpSchurImpl.cpp"
1104
1107
1108#ifdef PYTHON_SFD
1109
1110 boost::shared_ptr<ContactOps::SDFPython> sdf_python_ptr;
1111
1112 auto file_exists = [](std::string myfile) {
1113 std::ifstream file(myfile.c_str());
1114 if (file) {
1115 return true;
1116 }
1117 return false;
1118 };
1119
1120 if (file_exists("sdf.py")) {
1121 MOFEM_LOG("EP", Sev::inform) << "sdf.py file found";
1122 sdf_python_ptr = boost::make_shared<ContactOps::SDFPython>();
1123 CHKERR sdf_python_ptr->sdfInit("sdf.py");
1124 ContactOps::sdfPythonWeakPtr = sdf_python_ptr;
1125 } else {
1126 MOFEM_LOG("EP", Sev::warning) << "sdf.py file NOT found";
1127 }
1128#else
1129#endif
1130
1131 boost::shared_ptr<TsCtx> ts_ctx;
1133
1134 CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx.get(), PETSC_NULL);
1135
1136 boost::shared_ptr<FEMethod> monitor_ptr(new EshelbianMonitor(*this));
1137 ts_ctx->getLoopsMonitor().push_back(
1139
1140 CHKERR TSAppendOptionsPrefix(ts, "elastic_");
1141 CHKERR TSSetFromOptions(ts);
1142
1143 CHKERR TSSetDM(ts, dmElastic);
1144
1145 SNES snes;
1146 CHKERR TSGetSNES(ts, &snes);
1147
1148 PetscViewerAndFormat *vf;
1149 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
1150 PETSC_VIEWER_DEFAULT, &vf);
1151 CHKERR SNESMonitorSet(
1152 snes,
1153 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
1154 vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
1155
1156 PetscSection section;
1157 CHKERR DMGetSection(dmElastic, &section);
1158 int num_fields;
1159 CHKERR PetscSectionGetNumFields(section, &num_fields);
1160 for (int ff = 0; ff != num_fields; ff++) {
1161 const char *field_name;
1162 CHKERR PetscSectionGetFieldName(section, ff, &field_name);
1163 MOFEM_LOG_C("EP", Sev::inform, "Field %d name %s", ff, field_name);
1164 }
1165
1166 CHKERR DMoFEMMeshToLocalVector(dmElastic, x, INSERT_VALUES, SCATTER_FORWARD);
1167
1168 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
1169 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
1170
1171 // Adding field split solver
1172 boost::shared_ptr<SetUpSchur> schur_ptr;
1173
1174 if constexpr (A == AssemblyType::SCHUR) {
1175
1176 auto p = matDuplicate(m, MAT_DO_NOT_COPY_VALUES);
1177
1178 // If density is larger than zero, use dynamic time solver
1179 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
1180 CHKERR TSSetI2Function(ts, f, PETSC_NULL, PETSC_NULL);
1181 CHKERR TSSetI2Jacobian(ts, m, p, PETSC_NULL, PETSC_NULL);
1182 } else {
1183 CHKERR TSSetIFunction(ts, f, PETSC_NULL, PETSC_NULL);
1184 CHKERR TSSetIJacobian(ts, m, p, PETSC_NULL, PETSC_NULL);
1185 }
1186
1187 KSP ksp;
1188 CHKERR SNESGetKSP(snes, &ksp);
1189 PC pc;
1190 CHKERR KSPGetPC(ksp, &pc);
1191 schur_ptr = SetUpSchur::createSetUpSchur(
1192 mField, SmartPetscObj<Mat>(m, true), p, &*this);
1193 CHKERR schur_ptr->setUp(ksp);
1194
1195 elasticFeLhs->preProcessHook = [&]() {
1197 *(elasticFeLhs->matAssembleSwitch) = false;
1198 CHKERR schur_ptr->preProc();
1200 };
1201 elasticFeLhs->postProcessHook = [&]() {
1204 };
1205 elasticBcLhs->preProcessHook = [&]() {
1208 };
1209 elasticBcLhs->preProcessHook = [&]() {
1211 *(elasticBcLhs->matAssembleSwitch) = false;
1212 CHKERR schur_ptr->postProc();
1214 };
1215 }
1216
1217 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
1218 Vec xx;
1219 CHKERR VecDuplicate(x, &xx);
1220 CHKERR VecZeroEntries(xx);
1221 CHKERR TS2SetSolution(ts, x, xx);
1222 CHKERR VecDestroy(&xx);
1223 } else {
1224 CHKERR TSSetSolution(ts, x);
1225 }
1226
1227 CHKERR TSSetUp(ts);
1228 CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
1229 CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
1231 CHKERR TSSolve(ts, PETSC_NULL);
1233
1234 // CHKERR TSGetSNES(ts, &snes);
1235 int lin_solver_iterations;
1236 CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
1237 MOFEM_LOG("EP", Sev::inform)
1238 << "Number of linear solver iterations " << lin_solver_iterations;
1239
1240 PetscBool test_cook_flg = PETSC_FALSE;
1241 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test_cook", &test_cook_flg,
1242 PETSC_NULL);
1243 if (test_cook_flg) {
1244 constexpr int expected_lin_solver_iterations = 11;
1245 // if (lin_solver_iterations != expected_lin_solver_iterations)
1246 // SETERRQ2(
1247 // PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1248 // "Expected number of iterations is different than expected %d !=
1249 // %d", lin_solver_iterations, expected_lin_solver_iterations);
1250 }
1251
1253
1255}
1256
1257
1259 const std::string file) {
1261
1262 if (!dataAtPts) {
1263 dataAtPts =
1264 boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
1265 }
1266
1268
1269 auto domain_ops = [&](auto &fe) {
1271 fe.getUserPolynomialBase() =
1272 boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
1274 {HDIV, H1, L2});
1275 fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
1276 piolaStress, dataAtPts->getApproxPAtPts()));
1277 fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
1278 bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
1279 fe.getOpPtrVector().push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
1280 StretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
1281 fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1282 rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
1283 fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1284 spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
1285 fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1286 spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
1287 fe.getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
1288 spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
1289 // evaluate derived quantities
1290 fe.getOpPtrVector().push_back(
1292
1293 // evaluate integration points
1294 fe.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
1295 tag, true, false, dataAtPts, physicalEquations));
1296
1298 };
1299
1302 CHKERR domain_ops(*(op_loop_side->getSideFEPtr()));
1303 post_proc.getOpPtrVector().push_back(op_loop_side);
1304 post_proc.getOpPtrVector().push_back(new OpPostProcDataStructure(
1305 post_proc.getPostProcMesh(), post_proc.getMapGaussPts(), spatialL2Disp,
1306 dataAtPts));
1307
1308 CHKERR DMoFEMLoopFiniteElements(dM, skinElement.c_str(), &post_proc);
1309 CHKERR post_proc.writeFile(file.c_str());
1311}
1312
1313//! [Getting norms]
1316
1317 auto post_proc_norm_fe =
1318 boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
1319
1320 auto post_proc_norm_rule_hook = [](int, int, int p) -> int {
1321 return 2 * (p);
1322 };
1323 post_proc_norm_fe->getRuleHook = post_proc_norm_rule_hook;
1324
1325 post_proc_norm_fe->getUserPolynomialBase() =
1326 boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
1327
1329 post_proc_norm_fe->getOpPtrVector(), {L2, H1, HDIV});
1330
1331 enum NORMS { U_NORM_L2 = 0, U_NORM_H1, PIOLA_NORM, U_ERROR_L2, LAST_NORM };
1332 auto norms_vec =
1333 createVectorMPI(mField.get_comm(), LAST_NORM, PETSC_DETERMINE);
1334 CHKERR VecZeroEntries(norms_vec);
1335
1336 auto u_l2_ptr = boost::make_shared<MatrixDouble>();
1337 auto u_h1_ptr = boost::make_shared<MatrixDouble>();
1338 post_proc_norm_fe->getOpPtrVector().push_back(
1340 post_proc_norm_fe->getOpPtrVector().push_back(
1342 post_proc_norm_fe->getOpPtrVector().push_back(
1343 new OpCalcNormL2Tensor1<SPACE_DIM>(u_l2_ptr, norms_vec, U_NORM_L2));
1344 post_proc_norm_fe->getOpPtrVector().push_back(
1345 new OpCalcNormL2Tensor1<SPACE_DIM>(u_h1_ptr, norms_vec, U_NORM_H1));
1346 post_proc_norm_fe->getOpPtrVector().push_back(
1347 new OpCalcNormL2Tensor1<SPACE_DIM>(u_l2_ptr, norms_vec, U_ERROR_L2,
1348 u_h1_ptr));
1349
1350 auto piola_ptr = boost::make_shared<MatrixDouble>();
1351 post_proc_norm_fe->getOpPtrVector().push_back(
1353 post_proc_norm_fe->getOpPtrVector().push_back(
1354 new OpCalcNormL2Tensor2<3, 3>(piola_ptr, norms_vec, PIOLA_NORM));
1355
1357 *post_proc_norm_fe);
1358
1359 CHKERR VecAssemblyBegin(norms_vec);
1360 CHKERR VecAssemblyEnd(norms_vec);
1361 const double *norms;
1362 CHKERR VecGetArrayRead(norms_vec, &norms);
1363 MOFEM_LOG("EP", Sev::inform) << "norm_u: " << std::sqrt(norms[U_NORM_L2]);
1364 MOFEM_LOG("EP", Sev::inform) << "norm_u_h1: " << std::sqrt(norms[U_NORM_H1]);
1365 MOFEM_LOG("EP", Sev::inform)
1366 << "norm_error_u_l2: " << std::sqrt(norms[U_ERROR_L2]);
1367 MOFEM_LOG("EP", Sev::inform)
1368 << "norm_piola: " << std::sqrt(norms[PIOLA_NORM]);
1369 CHKERR VecRestoreArrayRead(norms_vec, &norms);
1370
1372}
1373//! [Getting norms]
1374
1377
1378 auto bc_mng = mField.getInterface<BcManager>();
1379 CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
1380 "", piolaStress, false, false);
1381
1382 bcSpatialDispVecPtr = boost::make_shared<BcDispVec>();
1383
1384 for (auto bc : bc_mng->getBcMapByBlockName()) {
1385 if (auto disp_bc = bc.second->dispBcPtr) {
1386
1387 MOFEM_LOG("EP", Sev::noisy) << *disp_bc;
1388
1389 std::vector<double> block_attributes(6, 0.);
1390 if (disp_bc->data.flag1 == 1) {
1391 block_attributes[0] = disp_bc->data.value1;
1392 block_attributes[3] = 1;
1393 }
1394 if (disp_bc->data.flag2 == 1) {
1395 block_attributes[1] = disp_bc->data.value2;
1396 block_attributes[4] = 1;
1397 }
1398 if (disp_bc->data.flag3 == 1) {
1399 block_attributes[2] = disp_bc->data.value3;
1400 block_attributes[5] = 1;
1401 }
1402 auto faces = bc.second->bcEnts.subset_by_dimension(2);
1403 bcSpatialDispVecPtr->emplace_back(bc.first, block_attributes, faces);
1404 }
1405 }
1406
1407 // old way of naming blocksets for displacement BCs
1408 CHKERR getBc(bcSpatialDispVecPtr, "SPATIAL_DISP_BC", 6);
1409
1411}
1412
1415
1416 auto bc_mng = mField.getInterface<BcManager>();
1417 CHKERR bc_mng->pushMarkDOFsOnEntities<ForceCubitBcData>("", piolaStress,
1418 false, false);
1419
1420 bcSpatialTraction = boost::make_shared<TractionBcVec>();
1421
1422 for (auto bc : bc_mng->getBcMapByBlockName()) {
1423 if (auto force_bc = bc.second->forceBcPtr) {
1424 std::vector<double> block_attributes(6, 0.);
1425 block_attributes[0] = -force_bc->data.value3 * force_bc->data.value1;
1426 block_attributes[3] = 1;
1427 block_attributes[1] = -force_bc->data.value4 * force_bc->data.value1;
1428 block_attributes[4] = 1;
1429 block_attributes[2] = -force_bc->data.value5 * force_bc->data.value1;
1430 block_attributes[5] = 1;
1431 auto faces = bc.second->bcEnts.subset_by_dimension(2);
1432 bcSpatialTraction->emplace_back(bc.first, block_attributes, faces);
1433 }
1434 }
1435
1436 CHKERR getBc(bcSpatialTraction, "SPATIAL_TRACTION_BC", 6);
1437
1439}
1440
1441} // namespace EshelbianPlasticity
1442
static Index< 'p', 3 > p
Implementation of tonsorial bubble base div(v) = 0.
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:10
Contains definition of EshelbianMonitor class.
Eshelbian plasticity interface.
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
cholesky decomposition
@ QUIET
Definition: definitions.h:208
@ NOISY
Definition: definitions.h:211
@ MF_ZERO
Definition: definitions.h:98
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ USER_BASE
user implemented approximation base
Definition: definitions.h:68
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#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
const int dim
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:319
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition: fem_tools.c:306
FTensor::Index< 'm', SPACE_DIM > m
static double phi
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1109
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 DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
Definition: DMMoFEM.cpp:786
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMoFEM.cpp:118
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 DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMoFEM.cpp:1128
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
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, 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 TS Jacobian evaluation function
Definition: DMMoFEM.cpp:839
PetscErrorCode DMMoFEMTSSetI2Jacobian(DM dm, 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 TS Jacobian evaluation function
Definition: DMMoFEM.cpp:1003
PetscErrorCode DMMoFEMTSSetI2Function(DM dm, 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 TS implicit function evaluation function
Definition: DMMoFEM.cpp:961
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
MoFEMErrorCode removeDofsOnEntities(const std::string problem_name, const std::string field_name, const Range ents, const int lo_coeff=0, const int hi_coeff=MAX_DOFS_ON_ENTITY, const int lo_order=0, const int hi_order=100, int verb=VERBOSE, const bool debug=false)
Remove DOFs from problem.
auto bit
set bit
const double v
phase velocity of light in medium (cm/ns)
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1884
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
MoFEMErrorCode CGG_BubbleBase_MBTET(const int p, const double *N, const double *diffN, FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > &phi, const int gdim)
Calculate CGGT tonsorial bubble base.
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:221
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMoFEM.cpp:424
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
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
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition: plastic.cpp:61
constexpr auto field_name
BcDisp(std::string name, std::vector< double > &attr, Range &faces)
BcRot(std::string name, std::vector< double > &attr, Range &faces)
MoFEMErrorCode getValueHdivForCGGBubble(MatrixDouble &pts)
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, BaseFunctionUnknownInterface **iface) const
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, SmartPetscObj< Mat > m, SmartPetscObj< Mat > p, EshelbianCore *ep_core_ptr)
static double inv_d_f_linear(const double v)
MoFEMErrorCode getSpatialDispBc()
[Getting norms]
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Getting interface of core database.
static boost::function< double(const double)> dd_f
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcLhs
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeLhs
static boost::function< double(const double)> inv_f
static boost::function< double(const double)> inv_dd_f
MoFEMErrorCode postProcessResults(const int tag, const std::string file)
MoFEMErrorCode addBoundaryFiniteElement(const EntityHandle meshset=0)
static boost::function< double(const double)> inv_d_f
MoFEMErrorCode addDMs(const BitRefLevel bit=BitRefLevel().set(0))
static double inv_f_log(const double v)
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeRhs
MoFEMErrorCode setElasticElementOps(const int tag)
static double inv_d_f_log(const double v)
static double d_f_linear(const double v)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
static double dd_f_log(const double v)
MoFEMErrorCode solveElastic(TS ts, Mat m, Vec f, Vec x)
static enum StretchSelector stretchSelector
static double inv_dd_f_log(const double v)
static boost::function< double(const double)> d_f
MoFEMErrorCode getBc(boost::shared_ptr< BC > &bc_vec_ptr, const std::string block_name, const int nb_attributes)
MoFEMErrorCode setBaseVolumeElementOps(const int tag, const bool do_rhs, const bool do_lhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > &fe)
boost::shared_ptr< BcRotVec > bcSpatialRotationVecPtr
MoFEMErrorCode setContactElementOps(boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_rhs, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_lhs)
MoFEMErrorCode setFaceElementOps(const bool add_elastic, const bool add_material, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_rhs, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_lhs)
EshelbianCore(MoFEM::Interface &m_field)
SmartPetscObj< DM > dM
Coupled problem all fields.
static double d_f_log(const double v)
MoFEMErrorCode addVolumeFiniteElement(const EntityHandle meshset=0)
SmartPetscObj< DM > dmElastic
Elastic problem.
boost::shared_ptr< BcDispVec > bcSpatialDispVecPtr
boost::shared_ptr< PhysicalEquations > physicalEquations
static double inv_dd_f_linear(const double v)
static boost::function< double(const double)> f
MoFEMErrorCode setVolumeElementOps(const int tag, const bool add_elastic, const bool add_material, boost::shared_ptr< VolumeElementForcesAndSourcesCore > &fe_rhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > &fe_lhs)
boost::shared_ptr< FaceElementForcesAndSourcesCore > contactRhs
static double dd_f_linear(const double v)
MoFEMErrorCode gettingNorms()
[Getting norms]
boost::shared_ptr< TractionFreeBc > bcSpatialFreeTraction
static double f_linear(const double v)
boost::shared_ptr< TractionBcVec > bcSpatialTraction
MoFEMErrorCode addFields(const EntityHandle meshset=0)
static double f_log(const double v)
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcRhs
SmartPetscObj< DM > dmPrjSpatial
Projection spatial displacement.
MoFEMErrorCode getTractionFreeBc(const EntityHandle meshset, boost::shared_ptr< TractionFreeBc > &bc_ptr, const std::string contact_set_name)
Remove all, but entities where kinematic constrains are applied.
static double inv_f_linear(const double v)
int operator()(int p_row, int p_col, int p_data) const
virtual MoFEMErrorCode evaluateRhs(EntData &data)=0
virtual MoFEMErrorCode evaluateLhs(EntData &data)=0
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
TractionBc(std::string name, std::vector< double > &attr, Range &faces)
Set integration rule on element.
int operator()(int p_row, int p_col, int p_data) const
Add operators pushing bases from local to physical configuration.
Base class if inherited used to calculate base functions.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
virtual moab::Interface & get_moab()=0
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MPI_Comm & get_comm() const =0
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
Deprecated interface functions.
Definition of the displacement bc data structure.
Definition: BCData.hpp:72
Class used to pass element data to calculate base functions on tet,triangle,edge.
const FieldApproximationBase bAse
Data on single entity (This is passed as argument to DataOperator::doWork)
data structure for finite element entity
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Provide data structure for (tensor) field approximation.
Definition of the force bc data structure.
Definition: BCData.hpp:134
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
Interface for managing meshsets containing materials and boundary conditions.
Natural boundary conditions.
Definition: Natural.hpp:57
Get norm of input MatrixDouble for Tensor1.
Get norm of input MatrixDouble for Tensor2.
Calculate tenor field using tensor base, i.e. Hdiv/Hcurl.
Calculate divergence of tonsorial field using vectorial base.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Calculate trace of vector (Hdiv/Hcurl) space.
Calculate symmetric tensor field rates ant integratio pts.
Calculate symmetric tensor field values at integration pts.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Approximate field values for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Element used to execute operators on side of the element.
Problem manager is used to build and partition problems.
intrusive_ptr for managing petsc objects
FEMethodsSequence & getLoopsMonitor()
Get the loops to do Monitor object.
Definition: TsCtx.hpp:108
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Base volume element used to integrate on skeleton.
static MoFEMErrorCode postStepInitialise(EshelbianCore *ep_ptr)
static MoFEMErrorCode postStepDestroy()
static MoFEMErrorCode preStepFun(TS ts)
static MoFEMErrorCode postStepFun(TS ts)