v0.14.0
Loading...
Searching...
No Matches
plastic.cpp
Go to the documentation of this file.
1/**
2 * \file plastic.cpp
3 * \example plastic.cpp
4 *
5 * Plasticity in 2d and 3d
6 *
7 */
8
9/* The above code is a preprocessor directive in C++ that checks if the macro
10"EXECUTABLE_DIMENSION" has been defined. If it has not been defined, it replaces
11the " */
12#ifndef EXECUTABLE_DIMENSION
13#define EXECUTABLE_DIMENSION 3
14#endif
15
16// #define ADD_CONTACT
17
18#include <MoFEM.hpp>
19#include <MatrixFunction.hpp>
20#include <IntegrationRules.hpp>
21
22using namespace MoFEM;
23
24template <int DIM> struct ElementsAndOps;
25
26template <> struct ElementsAndOps<2> {
30 static constexpr FieldSpace CONTACT_SPACE = HCURL;
31};
32
33template <> struct ElementsAndOps<3> {
37 static constexpr FieldSpace CONTACT_SPACE = HDIV;
38};
39
40constexpr int SPACE_DIM =
41 EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
42constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
43
44constexpr AssemblyType AT =
45 (SCHUR_ASSEMBLE) ? AssemblyType::SCHUR
46 : AssemblyType::PETSC; //< selected assembly type
48 IntegrationType::GAUSS; //< selected integration type
49
53
62
63#ifdef ADD_CONTACT
64//! [Specialisation for assembly]
65
66// Assemble to A matrix, by default, however, some terms are assembled only to
67// preconditioning.
68
69template <>
73 const EntitiesFieldData::EntData &row_data,
74 const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
75 return MatSetValues<AssemblyTypeSelector<AT>>(
76 op_ptr->getKSPA(), row_data, col_data, m, ADD_VALUES);
77 };
78
79template <>
83 const EntitiesFieldData::EntData &row_data,
84 const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
85 return MatSetValues<AssemblyTypeSelector<AT>>(
86 op_ptr->getKSPA(), row_data, col_data, m, ADD_VALUES);
87 };
88
89/**
90 * @brief Element used to specialise assembly
91 *
92 */
93struct BoundaryEleOpStab : public BoundaryEleOp {
94 using BoundaryEleOp::BoundaryEleOp;
95};
96
97/**
98 * @brief Specialise assembly for Stabilised matrix
99 *
100 * @tparam
101 */
102template <>
106 const EntitiesFieldData::EntData &row_data,
107 const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
108 return MatSetValues<AssemblyTypeSelector<AT>>(
109 op_ptr->getKSPB(), row_data, col_data, m, ADD_VALUES);
110 };
111//! [Specialisation for assembly]
112#endif // ADD_CONTACT
113
114inline double iso_hardening_exp(double tau, double b_iso) {
115 return std::exp(
116 std::max(static_cast<double>(std::numeric_limits<float>::min_exponent10),
117 -b_iso * tau));
118}
119
120/**
121 * Isotropic hardening
122 */
123inline double iso_hardening(double tau, double H, double Qinf, double b_iso,
124 double sigmaY) {
125 return H * tau + Qinf * (1. - iso_hardening_exp(tau, b_iso)) + sigmaY;
126}
127
128inline double iso_hardening_dtau(double tau, double H, double Qinf,
129 double b_iso) {
130 auto r = [&](auto tau) {
131 return H + Qinf * b_iso * iso_hardening_exp(tau, b_iso);
132 };
133 constexpr double eps = 1e-12;
134 return std::max(r(tau), eps * r(0));
135}
136
137/**
138 * Kinematic hardening
139 */
140template <typename T, int DIM>
141inline auto
143 double C1_k) {
147 if (C1_k < std::numeric_limits<double>::epsilon()) {
148 t_alpha(i, j) = 0;
149 return t_alpha;
150 }
151 t_alpha(i, j) = C1_k * t_plastic_strain(i, j);
152 return t_alpha;
153}
154
155template <int DIM>
163 t_diff(i, j, k, l) = C1_k * (t_kd(i, k) ^ t_kd(j, l)) / 4.;
164 return t_diff;
165}
166
167PetscBool is_large_strains = PETSC_TRUE; ///< Large strains
168PetscBool set_timer = PETSC_FALSE; ///< Set timer
169
170double scale = 1.;
171
172double young_modulus = 206913; ///< Young modulus
173double poisson_ratio = 0.29; ///< Poisson ratio
174double sigmaY = 450; ///< Yield stress
175double H = 129; ///< Hardening
176double visH = 0; ///< Viscous hardening
177double zeta = 5e-2; ///< Viscous hardening
178double Qinf = 265; ///< Saturation yield stress
179double b_iso = 16.93; ///< Saturation exponent
180double C1_k = 0; ///< Kinematic hardening
181
182double cn0 = 1;
183double cn1 = 1;
184
185int order = 2; ///< Order displacement
186int tau_order = order - 2; ///< Order of tau files
187int ep_order = order - 1; ///< Order of ep files
188int geom_order = 2; ///< Order if fixed.
189
190PetscBool is_quasi_static = PETSC_TRUE;
191double rho = 0.0;
192double alpha_damping = 0;
193
194#include <HenckyOps.hpp>
195#include <PlasticOps.hpp>
196#include <PlasticNaturalBCs.hpp>
197
198#ifdef ADD_CONTACT
199#ifdef PYTHON_SDF
200#include <boost/python.hpp>
201#include <boost/python/def.hpp>
202#include <boost/python/numpy.hpp>
203namespace bp = boost::python;
204namespace np = boost::python::numpy;
205#endif
206
207namespace ContactOps {
208
209double cn_contact = 0.1;
210
211}; // namespace ContactOps
212
213#include <ContactOps.hpp>
214#endif // ADD_CONTACT
215
225
226using namespace PlasticOps;
227using namespace HenckyOps;
228struct Example {
229
230 Example(MoFEM::Interface &m_field) : mField(m_field) {}
231
233
234private:
236
242
243 boost::shared_ptr<DomainEle> reactionFe;
244
245 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
246 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
247 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
248
251 double getScale(const double time) {
252 return scale * MoFEM::TimeScale::getScale(time);
253 };
254 };
255
256#ifdef ADD_CONTACT
257#ifdef PYTHON_SDF
258 boost::shared_ptr<ContactOps::SDFPython> sdfPythonPtr;
259#endif
260#endif // ADD_CONTACT
261};
262
263//! [Run problem]
268 CHKERR bC();
269 CHKERR OPs();
270 CHKERR tsSolve();
272}
273//! [Run problem]
274
275//! [Set up problem]
279
280 Range domain_ents;
281 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, domain_ents,
282 true);
283 auto get_ents_by_dim = [&](const auto dim) {
284 if (dim == SPACE_DIM) {
285 return domain_ents;
286 } else {
287 Range ents;
288 if (dim == 0)
289 CHKERR mField.get_moab().get_connectivity(domain_ents, ents, true);
290 else
291 CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents, true);
292 return ents;
293 }
294 };
295
296 auto get_base = [&]() {
297 auto domain_ents = get_ents_by_dim(SPACE_DIM);
298 if (domain_ents.empty())
299 CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Empty mesh");
300 const auto type = type_from_handle(domain_ents[0]);
301 switch (type) {
302 case MBQUAD:
304 case MBHEX:
306 case MBTRI:
308 case MBTET:
310 default:
311 CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Element type not handled");
312 }
313 return NOBASE;
314 };
315
316 const auto base = get_base();
317 MOFEM_LOG("PLASTICITY", Sev::inform)
318 << "Base " << ApproximationBaseNames[base];
319
320 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
321 CHKERR simple->addDomainField("EP", L2, base, size_symm);
322 CHKERR simple->addDomainField("TAU", L2, base, 1);
323 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
324
325 CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
326
327 PetscBool order_edge = PETSC_FALSE;
328 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-order_edge", &order_edge,
329 PETSC_NULL);
330 PetscBool order_face = PETSC_FALSE;
331 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-order_face", &order_face,
332 PETSC_NULL);
333 PetscBool order_volume = PETSC_FALSE;
334 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-order_volume", &order_volume,
335 PETSC_NULL);
336
337 if (order_edge || order_face || order_volume) {
338
339 MOFEM_LOG("PLASTICITY", Sev::inform) << "Order edge " << order_edge
340 ? "true"
341 : "false";
342 MOFEM_LOG("PLASTICITY", Sev::inform) << "Order face " << order_face
343 ? "true"
344 : "false";
345 MOFEM_LOG("PLASTICITY", Sev::inform) << "Order volume " << order_volume
346 ? "true"
347 : "false";
348
349 auto ents = get_ents_by_dim(0);
350 if (order_edge)
351 ents.merge(get_ents_by_dim(1));
352 if (order_face)
353 ents.merge(get_ents_by_dim(2));
354 if (order_volume)
355 ents.merge(get_ents_by_dim(3));
356 CHKERR simple->setFieldOrder("U", order, &ents);
357 } else {
358 CHKERR simple->setFieldOrder("U", order);
359 }
360 CHKERR simple->setFieldOrder("EP", ep_order);
361 CHKERR simple->setFieldOrder("TAU", tau_order);
362
363 CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
364
365#ifdef ADD_CONTACT
366 CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
367 SPACE_DIM);
368 CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
369 SPACE_DIM);
370
371 auto get_skin = [&]() {
372 Range body_ents;
373 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
374 Skinner skin(&mField.get_moab());
375 Range skin_ents;
376 CHKERR skin.find_skin(0, body_ents, false, skin_ents);
377 return skin_ents;
378 };
379
380 auto filter_blocks = [&](auto skin) {
381 bool is_contact_block = false;
382 Range contact_range;
383 for (auto m :
385
386 (boost::format("%s(.*)") % "CONTACT").str()
387
388 ))
389
390 ) {
391 is_contact_block =
392 true; ///< bloks interation is collectibe, so that is set irrespective
393 ///< if there are enerities in given rank or not in the block
394 MOFEM_LOG("CONTACT", Sev::inform)
395 << "Find contact block set: " << m->getName();
396 auto meshset = m->getMeshset();
397 Range contact_meshset_range;
398 CHKERR mField.get_moab().get_entities_by_dimension(
399 meshset, SPACE_DIM - 1, contact_meshset_range, true);
400
401 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
402 contact_meshset_range);
403 contact_range.merge(contact_meshset_range);
404 }
405 if (is_contact_block) {
406 MOFEM_LOG("SYNC", Sev::inform)
407 << "Nb entities in contact surface: " << contact_range.size();
409 skin = intersect(skin, contact_range);
410 }
411 return skin;
412 };
413
414 auto filter_true_skin = [&](auto skin) {
415 Range boundary_ents;
416 ParallelComm *pcomm =
417 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
418 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
419 PSTATUS_NOT, -1, &boundary_ents);
420 return boundary_ents;
421 };
422
423 auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
424 CHKERR simple->setFieldOrder("SIGMA", 0);
425 CHKERR simple->setFieldOrder("SIGMA", order - 1, &boundary_ents);
426#endif
427
428 CHKERR simple->setUp();
429 CHKERR simple->addFieldToEmptyFieldBlocks("U", "TAU");
430
431 auto project_ho_geometry = [&]() {
432 Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
433 return mField.loop_dofs("GEOMETRY", ent_method);
434 };
435 PetscBool project_geometry = PETSC_TRUE;
436 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-project_geometry",
437 &project_geometry, PETSC_NULL);
438 if (project_geometry){
439 CHKERR project_ho_geometry();
440 }
441
443}
444//! [Set up problem]
445
446//! [Create common data]
449
450 auto get_command_line_parameters = [&]() {
452
453 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-scale", &scale, PETSC_NULL);
454 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
455 &young_modulus, PETSC_NULL);
456 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
457 &poisson_ratio, PETSC_NULL);
458 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening", &H, PETSC_NULL);
459 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening_viscous", &visH,
460 PETSC_NULL);
461 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-yield_stress", &sigmaY,
462 PETSC_NULL);
463 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn0", &cn0, PETSC_NULL);
464 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn1", &cn1, PETSC_NULL);
465 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-zeta", &zeta, PETSC_NULL);
466 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-Qinf", &Qinf, PETSC_NULL);
467 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-b_iso", &b_iso, PETSC_NULL);
468 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-C1_k", &C1_k, PETSC_NULL);
469 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-large_strains",
470 &is_large_strains, PETSC_NULL);
471 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-set_timer", &set_timer,
472 PETSC_NULL);
473
474 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
475 PetscBool tau_order_is_set; ///< true if tau order is set
476 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-tau_order", &tau_order,
477 &tau_order_is_set);
478 PetscBool ep_order_is_set; ///< true if tau order is set
479 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-ep_order", &ep_order,
480 &ep_order_is_set);
481 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-geom_order", &geom_order,
482 PETSC_NULL);
483
484 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho", &rho, PETSC_NULL);
485 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-alpha_damping",
486 &alpha_damping, PETSC_NULL);
487
488 MOFEM_LOG("PLASTICITY", Sev::inform) << "Young modulus " << young_modulus;
489 MOFEM_LOG("PLASTICITY", Sev::inform) << "Poisson ratio " << poisson_ratio;
490 MOFEM_LOG("PLASTICITY", Sev::inform) << "Yield stress " << sigmaY;
491 MOFEM_LOG("PLASTICITY", Sev::inform) << "Hardening " << H;
492 MOFEM_LOG("PLASTICITY", Sev::inform) << "Viscous hardening " << visH;
493 MOFEM_LOG("PLASTICITY", Sev::inform) << "Saturation yield stress " << Qinf;
494 MOFEM_LOG("PLASTICITY", Sev::inform) << "Saturation exponent " << b_iso;
495 MOFEM_LOG("PLASTICITY", Sev::inform) << "Kinematic hardening " << C1_k;
496 MOFEM_LOG("PLASTICITY", Sev::inform) << "cn0 " << cn0;
497 MOFEM_LOG("PLASTICITY", Sev::inform) << "cn1 " << cn1;
498 MOFEM_LOG("PLASTICITY", Sev::inform) << "zeta " << zeta;
499
500 if (tau_order_is_set == PETSC_FALSE)
501 tau_order = order - 2;
502 if (ep_order_is_set == PETSC_FALSE)
503 ep_order = order - 1;
504
505 MOFEM_LOG("PLASTICITY", Sev::inform) << "Approximation order " << order;
506 MOFEM_LOG("PLASTICITY", Sev::inform)
507 << "Ep approximation order " << ep_order;
508 MOFEM_LOG("PLASTICITY", Sev::inform)
509 << "Tau approximation order " << tau_order;
510 MOFEM_LOG("PLASTICITY", Sev::inform)
511 << "Geometry approximation order " << geom_order;
512
513 MOFEM_LOG("PLASTICITY", Sev::inform) << "Density " << rho;
514 MOFEM_LOG("PLASTICITY", Sev::inform) << "alpha_damping " << alpha_damping;
515
516 PetscBool is_scale = PETSC_TRUE;
517 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_scale", &is_scale,
518 PETSC_NULL);
519 if (is_scale) {
521 }
522
523 MOFEM_LOG("PLASTICITY", Sev::inform) << "Scale " << scale;
524
525#ifdef ADD_CONTACT
526 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn_contact",
527 &ContactOps::cn_contact, PETSC_NULL);
528 MOFEM_LOG("CONTACT", Sev::inform)
529 << "cn_contact " << ContactOps::cn_contact;
530#endif // ADD_CONTACT
531
532 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-quasi_static",
533 &is_quasi_static, PETSC_NULL);
534 MOFEM_LOG("PLASTICITY", Sev::inform)
535 << "Is quasi static: " << (is_quasi_static ? "true" : "false");
536
538 };
539
540 CHKERR get_command_line_parameters();
541
542#ifdef ADD_CONTACT
543#ifdef PYTHON_SDF
544 sdfPythonPtr = boost::make_shared<ContactOps::SDFPython>();
545 CHKERR sdfPythonPtr->sdfInit("sdf.py");
546 ContactOps::sdfPythonWeakPtr = sdfPythonPtr;
547#endif
548#endif // ADD_CONTACT
549
551}
552//! [Create common data]
553
554//! [Boundary condition]
557
559 auto bc_mng = mField.getInterface<BcManager>();
560 auto prb_mng = mField.getInterface<ProblemsManager>();
561
562 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
563 "U", 0, 0);
564 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
565 "U", 1, 1);
566 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
567 "U", 2, 2);
568 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
569 "REMOVE_ALL", "U", 0, 3);
570
571#ifdef ADD_CONTACT
572 for (auto b : {"FIX_X", "REMOVE_X"})
573 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
574 "SIGMA", 0, 0, false, true);
575 for (auto b : {"FIX_Y", "REMOVE_Y"})
576 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
577 "SIGMA", 1, 1, false, true);
578 for (auto b : {"FIX_Z", "REMOVE_Z"})
579 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
580 "SIGMA", 2, 2, false, true);
581 for (auto b : {"FIX_ALL", "REMOVE_ALL"})
582 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
583 "SIGMA", 0, 3, false, true);
584 CHKERR bc_mng->removeBlockDOFsOnEntities(
585 simple->getProblemName(), "NO_CONTACT", "SIGMA", 0, 3, false, true);
586#endif
587
588 CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
589 simple->getProblemName(), "U");
590
591 auto &bc_map = bc_mng->getBcMapByBlockName();
592 for (auto bc : bc_map)
593 MOFEM_LOG("PLASTICITY", Sev::verbose) << "Marker " << bc.first;
594
596}
597//! [Boundary condition]
598
599//! [Push operators to pipeline]
602 auto pip_mng = mField.getInterface<PipelineManager>();
604 auto bc_mng = mField.getInterface<BcManager>();
605
606 auto integration_rule_bc = [](int, int, int ao) { return 2 * ao; };
607
608 auto vol_rule = [](int, int, int ao) { return 2 * ao + geom_order - 1; };
609
610 auto add_boundary_ops_lhs_mechanical = [&](auto &pip) {
612
614 "GEOMETRY");
615 pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
616
617 // Add Natural BCs to LHS
619 pip, mField, "U", Sev::inform);
620
621#ifdef ADD_CONTACT
622 CHKERR
623 ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
624 pip, "SIGMA", "U");
625 CHKERR
626 ContactOps::opFactoryBoundaryToDomainLhs<SPACE_DIM, AT, IT, DomainEle>(
627 mField, pip, simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",
628 vol_rule);
629#endif // ADD_CONTACT
630
632 };
633
634 auto add_boundary_ops_rhs_mechanical = [&](auto &pip) {
636
638 "GEOMETRY");
639 pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
640
641 // Add Natural BCs to RHS
643 pip, mField, "U", {boost::make_shared<ScaledTimeScale>()}, Sev::inform);
644
645#ifdef ADD_CONTACT
646 CHKERR ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, IT, BoundaryEleOp>(
647 pip, "SIGMA", "U");
648#endif // ADD_CONTACT
649
651 };
652
653 auto add_domain_ops_lhs = [this](auto &pip) {
656 "GEOMETRY");
657
658 if (is_quasi_static == PETSC_FALSE) {
659
660 //! [Only used for dynamics]
663 //! [Only used for dynamics]
664
665 auto get_inertia_and_mass_damping = [this](const double, const double,
666 const double) {
667 auto *pip = mField.getInterface<PipelineManager>();
668 auto &fe_domain_lhs = pip->getDomainLhsFE();
669 return (rho / scale) * fe_domain_lhs->ts_aa +
670 (alpha_damping / scale) * fe_domain_lhs->ts_a;
671 };
672 pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
673 }
674
675 CHKERR PlasticOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
676 mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
677
679 };
680
681 auto add_domain_ops_rhs = [this](auto &pip) {
683
685 "GEOMETRY");
686
688 pip, mField, "U",
689 {boost::make_shared<ScaledTimeScale>("body_force_hist.txt")},
690 Sev::inform);
691
692 // only in case of dynamics
693 if (is_quasi_static == PETSC_FALSE) {
694
695 //! [Only used for dynamics]
698 //! [Only used for dynamics]
699
700 auto mat_acceleration = boost::make_shared<MatrixDouble>();
702 "U", mat_acceleration));
703 pip.push_back(
704 new OpInertiaForce("U", mat_acceleration, [](double, double, double) {
705 return rho / scale;
706 }));
707 if (alpha_damping > 0) {
708 auto mat_velocity = boost::make_shared<MatrixDouble>();
709 pip.push_back(
710 new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", mat_velocity));
711 pip.push_back(
712 new OpInertiaForce("U", mat_velocity, [](double, double, double) {
713 return alpha_damping / scale;
714 }));
715 }
716 }
717
718 CHKERR PlasticOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
719 mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
720
721#ifdef ADD_CONTACT
722 CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
723 pip, "SIGMA", "U");
724#endif // ADD_CONTACT
725
727 };
728
729 CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
730 CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
731
732 // Boundary
733 CHKERR add_boundary_ops_lhs_mechanical(pip_mng->getOpBoundaryLhsPipeline());
734 CHKERR add_boundary_ops_rhs_mechanical(pip_mng->getOpBoundaryRhsPipeline());
735
736 CHKERR pip_mng->setDomainRhsIntegrationRule(vol_rule);
737 CHKERR pip_mng->setDomainLhsIntegrationRule(vol_rule);
738
739 CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_bc);
740 CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_bc);
741
742 auto create_reaction_pipeline = [&](auto &pip) {
745 "GEOMETRY");
746 CHKERR PlasticOps::opFactoryDomainReactions<SPACE_DIM, AT, IT, DomainEleOp>(
747 mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
749 };
750
751 reactionFe = boost::make_shared<DomainEle>(mField);
752 reactionFe->getRuleHook = vol_rule;
753 CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
754 reactionFe->postProcessHook =
756
758}
759//! [Push operators to pipeline]
760
761//! [Solve]
763
764 /**
765 * @brief Create data structure for handling Schur complement
766 *
767 * @param m_field
768 * @param sub_dm Schur complement sub dm
769 * @param field_split_it IS of Schur block
770 * @param ao_map AO map from sub dm to main problem
771 * @return boost::shared_ptr<SetUpSchur>
772 */
773 static boost::shared_ptr<SetUpSchur> createSetUpSchur(
774
775 MoFEM::Interface &m_field, SmartPetscObj<DM> sub_dm,
776 SmartPetscObj<IS> field_split_it, SmartPetscObj<AO> ao_map
777
778 );
779 virtual MoFEMErrorCode setUp(TS solver) = 0;
780
781protected:
782 SetUpSchur() = default;
783};
784
787
790 ISManager *is_manager = mField.getInterface<ISManager>();
791
792 auto snes_ctx_ptr = getDMSnesCtx(simple->getDM());
793
794 auto set_section_monitor = [&](auto solver) {
796 SNES snes;
797 CHKERR TSGetSNES(solver, &snes);
798 CHKERR SNESMonitorSet(snes,
799 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
801 (void *)(snes_ctx_ptr.get()), nullptr);
803 };
804
805 auto create_post_process_elements = [&]() {
806 auto pp_fe = boost::make_shared<PostProcEle>(mField);
807 auto &pip = pp_fe->getOpPtrVector();
808
809 auto push_vol_ops = [this](auto &pip) {
811 "GEOMETRY");
812
813 auto [common_plastic_ptr, common_henky_ptr] =
814 PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
815 mField, "MAT_PLASTIC", pip, "U", "EP", "TAU", 1., Sev::inform);
816
817 if (common_henky_ptr) {
818 if (common_plastic_ptr->mGradPtr != common_henky_ptr->matGradPtr)
819 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong pointer for grad");
820 }
821
822 return std::make_pair(common_plastic_ptr, common_henky_ptr);
823 };
824
825 auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
827
828 auto &pip = pp_fe->getOpPtrVector();
829
830 auto [common_plastic_ptr, common_henky_ptr] = p;
831
833
834 auto x_ptr = boost::make_shared<MatrixDouble>();
835 pip.push_back(
836 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
837 auto u_ptr = boost::make_shared<MatrixDouble>();
838 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
839
840 if (is_large_strains) {
841
842 pip.push_back(
843
844 new OpPPMap(
845
846 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
847
848 {{"PLASTIC_SURFACE",
849 common_plastic_ptr->getPlasticSurfacePtr()},
850 {"PLASTIC_MULTIPLIER",
851 common_plastic_ptr->getPlasticTauPtr()}},
852
853 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
854
855 {{"GRAD", common_henky_ptr->matGradPtr},
856 {"FIRST_PIOLA", common_henky_ptr->getMatFirstPiolaStress()}},
857
858 {{"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
859 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
860
861 )
862
863 );
864
865 } else {
866
867 pip.push_back(
868
869 new OpPPMap(
870
871 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
872
873 {{"PLASTIC_SURFACE",
874 common_plastic_ptr->getPlasticSurfacePtr()},
875 {"PLASTIC_MULTIPLIER",
876 common_plastic_ptr->getPlasticTauPtr()}},
877
878 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
879
880 {},
881
882 {{"STRAIN", common_plastic_ptr->mStrainPtr},
883 {"STRESS", common_plastic_ptr->mStressPtr},
884 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
885 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
886
887 )
888
889 );
890 }
891
893 };
894
895 auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
896 PetscBool post_proc_vol = PETSC_FALSE;
897 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_vol",
898 &post_proc_vol, PETSC_NULL);
899 if (post_proc_vol == PETSC_FALSE)
900 return boost::shared_ptr<PostProcEle>();
901 auto pp_fe = boost::make_shared<PostProcEle>(mField);
903 push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
904 "push_vol_post_proc_ops");
905 return pp_fe;
906 };
907
908 auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
909 PetscBool post_proc_skin = PETSC_TRUE;
910 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_skin",
911 &post_proc_skin, PETSC_NULL);
912 if (post_proc_skin == PETSC_FALSE)
913 return boost::shared_ptr<SkinPostProcEle>();
914
915 auto simple = mField.getInterface<Simple>();
916 auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
917 auto op_side = new OpLoopSide<SideEle>(
918 mField, simple->getDomainFEName(), SPACE_DIM, Sev::verbose);
919 pp_fe->getOpPtrVector().push_back(op_side);
920 CHK_MOAB_THROW(push_vol_post_proc_ops(
921 pp_fe, push_vol_ops(op_side->getOpPtrVector())),
922 "push_vol_post_proc_ops");
923 return pp_fe;
924 };
925
926 return std::make_pair(vol_post_proc(), skin_post_proc());
927 };
928
929 auto scatter_create = [&](auto D, auto coeff) {
931 CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
932 ROW, "U", coeff, coeff, is);
933 int loc_size;
934 CHKERR ISGetLocalSize(is, &loc_size);
935 Vec v;
936 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
937 VecScatter scatter;
938 CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
939 return std::make_tuple(SmartPetscObj<Vec>(v),
941 };
942
943 auto set_time_monitor = [&](auto dm, auto solver) {
945 boost::shared_ptr<Monitor> monitor_ptr(
946 new Monitor(dm, create_post_process_elements(), reactionFe, uXScatter,
947 uYScatter, uZScatter));
948 boost::shared_ptr<ForcesAndSourcesCore> null;
949 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
950 monitor_ptr, null, null);
952 };
953
954 auto set_schur_pc = [&](auto solver,
955 boost::shared_ptr<SetUpSchur> &schur_ptr) {
957
958 auto bc_mng = mField.getInterface<BcManager>();
959 auto name_prb = simple->getProblemName();
960
961 // create sub dm for Schur complement
962 auto create_sub_u_dm = [&](SmartPetscObj<DM> base_dm,
963 SmartPetscObj<DM> &dm_sub) {
965 dm_sub = createDM(mField.get_comm(), "DMMOFEM");
966 CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "SUB_U");
967 CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
968 CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
969 CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
970 for (auto f : {"U"}) {
973 }
974 CHKERR DMSetUp(dm_sub);
975
977 };
978
979 // Create nested (sub BC) Schur DM
980 if constexpr (AT == AssemblyType::SCHUR) {
981 SmartPetscObj<IS> is_epp;
982 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
983 simple->getProblemName(), ROW, "EP", 0, MAX_DOFS_ON_ENTITY, is_epp);
984 SmartPetscObj<IS> is_tau;
985 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
986 simple->getProblemName(), ROW, "TAU", 0, MAX_DOFS_ON_ENTITY, is_tau);
987
988 IS is_union_raw;
989 CHKERR ISExpand(is_epp, is_tau, &is_union_raw);
990 SmartPetscObj<IS> is_union(is_union_raw);
991
992#ifdef ADD_CONTACT
993 auto add_sigma_to_is = [&](auto is_union) {
994 SmartPetscObj<IS> is_union_sigma;
995 auto add_sigma_to_is_impl = [&]() {
997 SmartPetscObj<IS> is_sigma;
998 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
999 simple->getProblemName(), ROW, "SIGMA", 0, MAX_DOFS_ON_ENTITY,
1000 is_sigma);
1001 IS is_union_raw_sigma;
1002 CHKERR ISExpand(is_union, is_sigma, &is_union_raw_sigma);
1003 is_union_sigma = SmartPetscObj<IS>(is_union_raw_sigma);
1005 };
1006 CHK_THROW_MESSAGE(add_sigma_to_is_impl(), "Can not add sigma to IS");
1007 return is_union_sigma;
1008 };
1009 is_union = add_sigma_to_is(is_union);
1010#endif // ADD_CONTACT
1011
1012 SmartPetscObj<DM> dm_u_sub;
1013 CHKERR create_sub_u_dm(simple->getDM(), dm_u_sub);
1014
1015 // Indices has to be map fro very to level, while assembling Schur
1016 // complement.
1017 auto is_up = getDMSubData(dm_u_sub)->getSmartRowIs();
1018 auto ao_up = createAOMappingIS(is_up, PETSC_NULL);
1019 schur_ptr =
1020 SetUpSchur::createSetUpSchur(mField, dm_u_sub, is_union, ao_up);
1021 CHKERR schur_ptr->setUp(solver);
1022 }
1023
1025 };
1026
1027 auto dm = simple->getDM();
1028 auto D = createDMVector(dm);
1029 auto DD = vectorDuplicate(D);
1030 uXScatter = scatter_create(D, 0);
1031 uYScatter = scatter_create(D, 1);
1032 if constexpr (SPACE_DIM == 3)
1033 uZScatter = scatter_create(D, 2);
1034
1035 auto create_solver = [pip_mng]() {
1036 if (is_quasi_static == PETSC_TRUE)
1037 return pip_mng->createTSIM();
1038 else
1039 return pip_mng->createTSIM2();
1040 };
1041
1042 auto solver = create_solver();
1043
1044 auto active_pre_lhs = []() {
1046 std::fill(PlasticOps::CommonData::activityData.begin(),
1049 };
1050
1051 auto active_post_lhs = [&]() {
1053 auto get_iter = [&]() {
1054 SNES snes;
1055 CHK_THROW_MESSAGE(TSGetSNES(solver, &snes), "Can not get SNES");
1056 int iter;
1057 CHK_THROW_MESSAGE(SNESGetIterationNumber(snes, &iter),
1058 "Can not get iter");
1059 return iter;
1060 };
1061
1062 auto iter = get_iter();
1063 if (iter >= 0) {
1064
1065 std::array<int, 5> activity_data;
1066 std::fill(activity_data.begin(), activity_data.end(), 0);
1067 MPI_Allreduce(PlasticOps::CommonData::activityData.data(),
1068 activity_data.data(), activity_data.size(), MPI_INT,
1069 MPI_SUM, mField.get_comm());
1070
1071 int &active_points = activity_data[0];
1072 int &avtive_full_elems = activity_data[1];
1073 int &avtive_elems = activity_data[2];
1074 int &nb_points = activity_data[3];
1075 int &nb_elements = activity_data[4];
1076
1077 if (nb_points) {
1078
1079 double proc_nb_points =
1080 100 * static_cast<double>(active_points) / nb_points;
1081 double proc_nb_active =
1082 100 * static_cast<double>(avtive_elems) / nb_elements;
1083 double proc_nb_full_active = 100;
1084 if (avtive_elems)
1085 proc_nb_full_active =
1086 100 * static_cast<double>(avtive_full_elems) / avtive_elems;
1087
1088 MOFEM_LOG_C("PLASTICITY", Sev::inform,
1089 "Iter %d nb pts %d nb avtive pts %d (%3.3f\%) nb active "
1090 "elements %d "
1091 "(%3.3f\%) nb full active elems %d (%3.3f\%)",
1092 iter, nb_points, active_points, proc_nb_points,
1093 avtive_elems, proc_nb_active, avtive_full_elems,
1094 proc_nb_full_active, iter);
1095 }
1096 }
1097
1099 };
1100
1101 auto add_active_dofs_elem = [&](auto dm) {
1103 auto fe_pre_proc = boost::make_shared<FEMethod>();
1104 fe_pre_proc->preProcessHook = active_pre_lhs;
1105 auto fe_post_proc = boost::make_shared<FEMethod>();
1106 fe_post_proc->postProcessHook = active_post_lhs;
1107 auto ts_ctx_ptr = getDMTsCtx(dm);
1108 ts_ctx_ptr->getPreProcessIJacobian().push_front(fe_pre_proc);
1109 ts_ctx_ptr->getPostProcessIJacobian().push_back(fe_post_proc);
1111 };
1112
1113 auto set_essential_bc = [&](auto dm, auto solver) {
1115 // This is low level pushing finite elements (pipelines) to solver
1116
1117 auto pre_proc_ptr = boost::make_shared<FEMethod>();
1118 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
1119 auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
1120
1121 // Add boundary condition scaling
1122 auto disp_time_scale = boost::make_shared<TimeScale>();
1123
1124 auto get_bc_hook_rhs = [this, pre_proc_ptr, disp_time_scale]() {
1125 EssentialPreProc<DisplacementCubitBcData> hook(mField, pre_proc_ptr,
1126 {disp_time_scale}, false);
1127 return hook;
1128 };
1129 pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
1130
1131 auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
1134 mField, post_proc_rhs_ptr, nullptr, Sev::verbose)();
1136 mField, post_proc_rhs_ptr, 1.)();
1138 };
1139 auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
1141 mField, post_proc_lhs_ptr, 1.);
1142 };
1143 post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
1144
1145 auto ts_ctx_ptr = getDMTsCtx(dm);
1146 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
1147 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
1148 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
1149
1150 SNES snes;
1151 CHKERR TSGetSNES(solver, &snes);
1152 KSP ksp;
1153 CHKERR SNESGetKSP(snes, &ksp);
1154 PC pc;
1155 CHKERR KSPGetPC(ksp, &pc);
1156 PetscBool is_pcfs = PETSC_FALSE;
1157 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1158
1159 if (is_pcfs == PETSC_FALSE) {
1160 post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
1161 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
1162 }
1164 };
1165
1166 if (is_quasi_static == PETSC_TRUE) {
1167 CHKERR TSSetSolution(solver, D);
1168 } else {
1169 CHKERR TS2SetSolution(solver, D, DD);
1170 }
1171
1172 CHKERR set_section_monitor(solver);
1173 CHKERR set_time_monitor(dm, solver);
1174 CHKERR TSSetFromOptions(solver);
1175 CHKERR TSSetUp(solver);
1176
1177 CHKERR add_active_dofs_elem(dm);
1178 boost::shared_ptr<SetUpSchur> schur_ptr;
1179 CHKERR set_schur_pc(solver, schur_ptr);
1180 CHKERR set_essential_bc(dm, solver);
1181
1182 MOFEM_LOG_CHANNEL("TIMER");
1183 MOFEM_LOG_TAG("TIMER", "timer");
1184 if (set_timer)
1185 BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
1186 MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp";
1187 CHKERR TSSetUp(solver);
1188 MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp <= done";
1189 MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve";
1190 CHKERR TSSolve(solver, NULL);
1191 MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve <= done";
1192
1194}
1195//! [Solve]
1196
1197static char help[] = "...\n\n";
1198
1199int main(int argc, char *argv[]) {
1200
1201#ifdef ADD_CONTACT
1202#ifdef PYTHON_SDF
1203 Py_Initialize();
1204 np::initialize();
1205#endif
1206#endif // ADD_CONTACT
1207
1208 // Initialisation of MoFEM/PETSc and MOAB data structures
1209 const char param_file[] = "param_file.petsc";
1210 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
1211
1212 // Add logging channel for example
1213 auto core_log = logging::core::get();
1214 core_log->add_sink(
1216 core_log->add_sink(
1218 LogManager::setLog("PLASTICITY");
1219 MOFEM_LOG_TAG("PLASTICITY", "Plasticity");
1220
1221#ifdef ADD_CONTACT
1222 core_log->add_sink(
1224 LogManager::setLog("CONTACT");
1225 MOFEM_LOG_TAG("CONTACT", "Contact");
1226#endif // ADD_CONTACT
1227
1228 try {
1229
1230 //! [Register MoFEM discrete manager in PETSc]
1231 DMType dm_name = "DMMOFEM";
1232 CHKERR DMRegister_MoFEM(dm_name);
1233 //! [Register MoFEM discrete manager in PETSc
1234
1235 //! [Create MoAB]
1236 moab::Core mb_instance; ///< mesh database
1237 moab::Interface &moab = mb_instance; ///< mesh database interface
1238 //! [Create MoAB]
1239
1240 //! [Create MoFEM]
1241 MoFEM::Core core(moab); ///< finite element database
1242 MoFEM::Interface &m_field = core; ///< finite element database interface
1243 //! [Create MoFEM]
1244
1245 //! [Load mesh]
1246 Simple *simple = m_field.getInterface<Simple>();
1247 CHKERR simple->getOptions();
1248 CHKERR simple->loadFile();
1249 //! [Load mesh]
1250
1251 //! [Example]
1252 Example ex(m_field);
1253 CHKERR ex.runProblem();
1254 //! [Example]
1255 }
1257
1259
1260#ifdef ADD_CONTACT
1261#ifdef PYTHON_SDF
1262 if (Py_FinalizeEx() < 0) {
1263 exit(120);
1264 }
1265#endif
1266#endif // ADD_CONTACT
1267
1268 return 0;
1269}
1270
1271struct SetUpSchurImpl : public SetUpSchur {
1272
1274 SmartPetscObj<IS> field_split_is, SmartPetscObj<AO> ao_up)
1275 : SetUpSchur(), mField(m_field), subDM(sub_dm),
1276 fieldSplitIS(field_split_is), aoUp(ao_up) {
1277 if (S) {
1280 "Is expected that schur matrix is not allocated. This is "
1281 "possible only is PC is set up twice");
1282 }
1283 }
1285#ifdef ADD_CONTACT
1286 A.reset();
1287 P.reset();
1288#endif // ADD_CONTACT
1289 S.reset();
1290 }
1291
1292 MoFEMErrorCode setUp(TS solver);
1295
1296private:
1297#ifdef ADD_CONTACT
1300#endif // ADD_CONTACT
1302
1304 SmartPetscObj<DM> subDM; ///< field split sub dm
1305 SmartPetscObj<IS> fieldSplitIS; ///< IS for split Schur block
1306 SmartPetscObj<AO> aoUp; ///> main DM to subDM
1307};
1308
1311 auto simple = mField.getInterface<Simple>();
1312 auto pip_mng = mField.getInterface<PipelineManager>();
1313
1314 SNES snes;
1315 CHKERR TSGetSNES(solver, &snes);
1316 KSP ksp;
1317 CHKERR SNESGetKSP(snes, &ksp);
1318 CHKERR KSPSetFromOptions(ksp);
1319
1320 PC pc;
1321 CHKERR KSPSetFromOptions(ksp);
1322 CHKERR KSPGetPC(ksp, &pc);
1323 PetscBool is_pcfs = PETSC_FALSE;
1324 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1325 if (is_pcfs) {
1326 if (S) {
1329 "Is expected that schur matrix is not allocated. This is "
1330 "possible only is PC is set up twice");
1331 }
1332
1333#ifdef ADD_CONTACT
1334 auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
1335 A = createDMMatrix(simple->getDM());
1336 P = matDuplicate(A, MAT_DO_NOT_COPY_VALUES);
1337 CHKERR TSSetIJacobian(solver, A, P, TsSetIJacobian, ts_ctx_ptr.get());
1338#endif // ADD_CONTACT
1340 CHKERR MatSetBlockSize(S, SPACE_DIM);
1341
1342 auto set_ops = [&]() {
1344 auto pip_mng = mField.getInterface<PipelineManager>();
1345
1346#ifndef ADD_CONTACT
1347 // Boundary
1348 pip_mng->getOpBoundaryLhsPipeline().push_front(
1349 new OpSchurAssembleBegin());
1350 pip_mng->getOpBoundaryLhsPipeline().push_back(
1352
1353 {"EP", "TAU"}, {nullptr, nullptr}, {SmartPetscObj<AO>(), aoUp},
1354 {SmartPetscObj<Mat>(), S}, {false, false}
1355
1356 ));
1357 // Domain
1358 pip_mng->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
1359 pip_mng->getOpDomainLhsPipeline().push_back(
1361
1362 {"EP", "TAU"}, {nullptr, nullptr}, {SmartPetscObj<AO>(), aoUp},
1363 {SmartPetscObj<Mat>(), S}, {false, false}
1364
1365 ));
1366#else
1367
1368 double eps_stab = 1e-4;
1369 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,
1370 PETSC_NULL);
1371
1374 using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
1375
1376 // Boundary
1377 pip_mng->getOpBoundaryLhsPipeline().push_front(
1378 new OpSchurAssembleBegin());
1379 pip_mng->getOpBoundaryLhsPipeline().push_back(
1380 new OpMassStab("SIGMA", "SIGMA", [eps_stab](double, double, double) {
1381 return eps_stab;
1382 }));
1383 pip_mng->getOpBoundaryLhsPipeline().push_back(
1385
1386 {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr},
1389 {false, false, false}
1390
1391 ));
1392 // Domain
1393 pip_mng->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
1394 pip_mng->getOpDomainLhsPipeline().push_back(
1396
1397 {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr},
1400 {false, false, false}
1401
1402 ));
1403#endif // ADD_CONTACT
1405 };
1406
1407 auto set_assemble_elems = [&]() {
1409 auto schur_asmb_pre_proc = boost::make_shared<FEMethod>();
1410 schur_asmb_pre_proc->preProcessHook = [this]() {
1412#ifdef ADD_CONTACT
1413 CHKERR MatZeroEntries(A);
1414 CHKERR MatZeroEntries(P);
1415#endif // ADD_CONTACT
1416 CHKERR MatZeroEntries(S);
1417 MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble Begin";
1419 };
1420 auto schur_asmb_post_proc = boost::make_shared<FEMethod>();
1421
1422 schur_asmb_post_proc->postProcessHook = [this, schur_asmb_post_proc]() {
1424 MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble End";
1425
1426#ifndef ADD_CONTACT
1428 mField, schur_asmb_post_proc, 1)();
1429#else // ADD_CONTACT
1430 CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1431 CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1432 // Apply essential constrains to A matrix
1434 mField, schur_asmb_post_proc, 1, A)();
1435 CHKERR MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
1436 CHKERR MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
1437 CHKERR MatAXPY(P, 1, A, SAME_NONZERO_PATTERN);
1438#endif // ADD_CONTACT
1439
1440 // Apply essential constrains to Schur complement
1441 CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
1442 CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
1444 mField, schur_asmb_post_proc, 1, S, aoUp)();
1445
1447 };
1448 auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
1449 ts_ctx_ptr->getPreProcessIJacobian().push_front(schur_asmb_pre_proc);
1450 ts_ctx_ptr->getPostProcessIJacobian().push_front(schur_asmb_post_proc);
1452 };
1453
1454 auto set_pc = [&]() {
1456 CHKERR PCFieldSplitSetIS(pc, NULL, fieldSplitIS);
1457 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1459 };
1460
1461 CHKERR set_ops();
1462 CHKERR set_pc();
1463 CHKERR set_assemble_elems();
1464
1465 } else {
1466 pip_mng->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
1467 pip_mng->getOpBoundaryLhsPipeline().push_back(
1468 new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
1469 pip_mng->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
1470 pip_mng->getOpDomainLhsPipeline().push_back(
1471 new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
1472 }
1473
1474 // we do not those anymore
1475 subDM.reset();
1476 fieldSplitIS.reset();
1477 // aoUp.reset();
1479}
1480
1481boost::shared_ptr<SetUpSchur>
1483 SmartPetscObj<DM> sub_dm, SmartPetscObj<IS> is_sub,
1484 SmartPetscObj<AO> ao_up) {
1485 return boost::shared_ptr<SetUpSchur>(
1486 new SetUpSchurImpl(m_field, sub_dm, is_sub, ao_up));
1487}
static Index< 'p', 3 > p
std::string param_file
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
int main()
Definition: adol-c_atom.cpp:46
PostProcEleByDim< SPACE_DIM >::SideEle SideEle
static const double eps
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
#define SCHUR_ASSEMBLE
Definition: contact.cpp:18
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
#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
@ NOBASE
Definition: definitions.h:59
@ 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
FieldSpace
approximation spaces
Definition: definitions.h:82
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
@ 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 CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
static const char *const ApproximationBaseNames[]
Definition: definitions.h:72
#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
constexpr int SPACE_DIM
[Define dimension]
Definition: elastic.cpp:18
const int dim
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FTensor::Index< 'm', SPACE_DIM > m
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 DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:242
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
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
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
#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
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1869
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMoFEM.cpp:1042
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
Definition: TsCtx.cpp:165
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1045
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition: SnesCtx.cpp:232
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
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.
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:1031
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double young_modulus
Young modulus.
Definition: plastic.cpp:172
constexpr AssemblyType AT
Definition: plastic.cpp:44
double C1_k
Kinematic hardening.
Definition: plastic.cpp:180
double Qinf
Saturation yield stress.
Definition: plastic.cpp:178
constexpr IntegrationType IT
Definition: plastic.cpp:47
static char help[]
[Solve]
Definition: plastic.cpp:1197
double rho
Definition: plastic.cpp:191
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
PetscBool is_quasi_static
Definition: plastic.cpp:190
double alpha_damping
Definition: plastic.cpp:192
constexpr int SPACE_DIM
Definition: plastic.cpp:40
double visH
Viscous hardening.
Definition: plastic.cpp:176
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:173
auto kinematic_hardening(FTensor::Tensor2_symmetric< T, DIM > &t_plastic_strain, double C1_k)
Definition: plastic.cpp:142
PetscBool set_timer
Set timer.
Definition: plastic.cpp:168
double iso_hardening_dtau(double tau, double H, double Qinf, double b_iso)
Definition: plastic.cpp:128
double scale
Definition: plastic.cpp:170
constexpr auto size_symm
Definition: plastic.cpp:42
double zeta
Viscous hardening.
Definition: plastic.cpp:177
double H
Hardening.
Definition: plastic.cpp:175
int tau_order
Order of tau files.
Definition: plastic.cpp:186
double iso_hardening_exp(double tau, double b_iso)
Definition: plastic.cpp:114
double cn0
Definition: plastic.cpp:182
int order
Order displacement.
Definition: plastic.cpp:185
double b_iso
Saturation exponent.
Definition: plastic.cpp:179
PetscBool is_large_strains
Large strains.
Definition: plastic.cpp:167
int geom_order
Order if fixed.
Definition: plastic.cpp:188
double sigmaY
Yield stress.
Definition: plastic.cpp:174
double iso_hardening(double tau, double H, double Qinf, double b_iso, double sigmaY)
Definition: plastic.cpp:123
auto kinematic_hardening_dplastic_strain(double C1_k)
Definition: plastic.cpp:156
int ep_order
Order of ep files.
Definition: plastic.cpp:187
double cn1
Definition: plastic.cpp:183
constexpr FieldSpace CONTACT_SPACE
Definition: plastic.cpp:52
FTensor::Index< 'm', 3 > m
Element used to specialise assembly.
Definition: contact.cpp:94
PipelineManager::ElementsAndOpsByDim< 2 >::FaceSideEle SideEle
Definition: plastic.cpp:29
PipelineManager::ElementsAndOpsByDim< 3 >::FaceSideEle SideEle
Definition: plastic.cpp:36
double getScale(const double time)
Get scaling at a given time.
Definition: plastic.cpp:251
[Example]
Definition: plastic.cpp:228
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:247
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:245
MoFEMErrorCode tsSolve()
Definition: plastic.cpp:785
FieldApproximationBase base
Definition: plot_base.cpp:68
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:447
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: plastic.cpp:246
Example(MoFEM::Interface &m_field)
Definition: plastic.cpp:230
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:600
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:264
MoFEM::Interface & mField
Definition: plastic.cpp:235
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:276
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:555
boost::shared_ptr< DomainEle > reactionFe
Definition: plastic.cpp:243
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
Data on single entity (This is passed as argument to DataOperator::doWork)
Class (Function) to enforce essential constrains on the left hand side diagonal.
Definition: Essential.hpp:33
Class (Function) to enforce essential constrains on the right hand side diagonal.
Definition: Essential.hpp:41
Class (Function) to enforce essential constrains.
Definition: Essential.hpp:25
Class (Function) to calculate residual side diagonal.
Definition: Essential.hpp:49
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
Interface for managing meshsets containing materials and boundary conditions.
Assembly methods.
Definition: Natural.hpp:65
boost::function< MoFEMErrorCode(ForcesAndSourcesCore::UserDataOperator *op_ptr, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, MatrixDouble &m)> MatSetValuesHook
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.
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.
MoFEM::VolumeElementForcesAndSourcesCore VolEle
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Problem manager is used to build and partition problems.
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
Force scale operator for reading two columns.
double getScale(const double time)
Get scaling at a given time.
TimeScale(std::string file_name="", bool error_if_file_not_given=false)
TimeScale constructor.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
[Push operators to pipeline]
static std::array< int, 5 > activityData
Definition: PlasticOps.hpp:107
[Push operators to pipeline]
Definition: plastic.cpp:762
SetUpSchur()=default
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, SmartPetscObj< DM > sub_dm, SmartPetscObj< IS > field_split_it, SmartPetscObj< AO > ao_map)
Create data structure for handling Schur complement.
Definition: plastic.cpp:1482
virtual MoFEMErrorCode setUp(TS solver)=0
SmartPetscObj< Mat > A
Definition: contact.cpp:983
SmartPetscObj< DM > subDM
field split sub dm
Definition: plastic.cpp:1304
MoFEMErrorCode setUp(KSP solver)
SmartPetscObj< Mat > S
SmartPetscObj< Mat > P
SetUpSchurImpl(MoFEM::Interface &m_field, SmartPetscObj< DM > sub_dm, SmartPetscObj< IS > field_split_is, SmartPetscObj< AO > ao_up)
Definition: plastic.cpp:1273
SmartPetscObj< IS > fieldSplitIS
IS for split Schur block.
Definition: plastic.cpp:1305
virtual ~SetUpSchurImpl()
Definition: plastic.cpp:1284
MoFEMErrorCode postProc()
SmartPetscObj< AO > aoUp
MoFEMErrorCode preProc()
MoFEM::Interface & mField