v0.15.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// #undef 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> {
27 using DomainEle = PipelineManager::FaceEle;
28 using BoundaryEle = PipelineManager::EdgeEle;
29 using SideEle = PipelineManager::ElementsAndOpsByDim<2>::FaceSideEle;
30 static constexpr FieldSpace CONTACT_SPACE = HCURL;
31};
32
33template <> struct ElementsAndOps<3> {
34 using DomainEle = PipelineManager::VolEle;
35 using BoundaryEle = PipelineManager::FaceEle;
36 using SideEle = PipelineManager::ElementsAndOpsByDim<3>::FaceSideEle;
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::BLOCK_SCHUR
46 : AssemblyType::PETSC; //< selected assembly type
47constexpr IntegrationType IT =
48 IntegrationType::GAUSS; //< selected integration type
49
50constexpr FieldSpace ElementsAndOps<2>::CONTACT_SPACE;
51constexpr FieldSpace ElementsAndOps<3>::CONTACT_SPACE;
52constexpr FieldSpace CONTACT_SPACE = ElementsAndOps<SPACE_DIM>::CONTACT_SPACE;
53
54using EntData = EntitiesFieldData::EntData;
55using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle;
56using DomainEleOp = DomainEle::UserDataOperator;
57using BoundaryEle = ElementsAndOps<SPACE_DIM>::BoundaryEle;
58using BoundaryEleOp = BoundaryEle::UserDataOperator;
59using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;
60using SkinPostProcEle = PostProcBrokenMeshInMoab<BoundaryEle>;
61using SideEle = ElementsAndOps<SPACE_DIM>::SideEle;
62using SetPtsData = FieldEvaluatorInterface::SetPtsData;
63
64inline double iso_hardening_exp(double tau, double b_iso) {
65 return std::exp(
66 std::max(static_cast<double>(std::numeric_limits<float>::min_exponent10),
67 -b_iso * tau));
68}
69
70/**
71 * Isotropic hardening
72 */
73inline double iso_hardening(double tau, double H, double Qinf, double b_iso,
74 double sigmaY) {
75 return H * tau + Qinf * (1. - iso_hardening_exp(tau, b_iso)) + sigmaY;
76}
77
78inline double iso_hardening_dtau(double tau, double H, double Qinf,
79 double b_iso) {
80 auto r = [&](auto tau) {
81 return H + Qinf * b_iso * iso_hardening_exp(tau, b_iso);
82 };
83 constexpr double eps = 1e-12;
84 return std::max(r(tau), eps * r(0));
85}
86
87/**
88 * Kinematic hardening
89 */
90template <typename T, int DIM>
91inline auto
92kinematic_hardening(FTensor::Tensor2_symmetric<T, DIM> &t_plastic_strain,
93 double C1_k) {
94 FTensor::Index<'i', DIM> i;
95 FTensor::Index<'j', DIM> j;
96 FTensor::Tensor2_symmetric<double, DIM> t_alpha;
97 if (C1_k < std::numeric_limits<double>::epsilon()) {
98 t_alpha(i, j) = 0;
99 return t_alpha;
100 }
101 t_alpha(i, j) = C1_k * t_plastic_strain(i, j);
102 return t_alpha;
103}
104
105template <int DIM>
106inline auto kinematic_hardening_dplastic_strain(double C1_k) {
107 FTensor::Index<'i', DIM> i;
108 FTensor::Index<'j', DIM> j;
109 FTensor::Index<'k', DIM> k;
110 FTensor::Index<'l', DIM> l;
111 FTensor::Ddg<double, DIM, DIM> t_diff;
112 constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
113 t_diff(i, j, k, l) = C1_k * (t_kd(i, k) ^ t_kd(j, l)) / 4.;
114 return t_diff;
115}
116
117PetscBool is_large_strains = PETSC_TRUE; ///< Large strains
118PetscBool set_timer = PETSC_FALSE; ///< Set timer
119PetscBool do_eval_field = PETSC_FALSE; ///< Evaluate field
120
121int atom_test = 0; ///< Atom test
122
123double scale = 1.;
124
125double young_modulus = 206913; ///< Young modulus
126double poisson_ratio = 0.29; ///< Poisson ratio
127double sigmaY = 450; ///< Yield stress
128double H = 129; ///< Hardening
129double visH = 0; ///< Viscous hardening
130double zeta = 5e-2; ///< Viscous hardening
131double Qinf = 265; ///< Saturation yield stress
132double b_iso = 16.93; ///< Saturation exponent
133double C1_k = 0; ///< Kinematic hardening
134
135double cn0 = 1;
136double cn1 = 1;
137
138int order = 2; ///< Order displacement
139int tau_order = order - 2; ///< Order of tau files
140int ep_order = order - 1; ///< Order of ep files
141int geom_order = 2; ///< Order if fixed.
142
143PetscBool is_quasi_static = PETSC_TRUE;
144double rho = 0.0;
145double alpha_damping = 0;
146
147#include <HenckyOps.hpp>
148#include <PlasticOps.hpp>
149#include <PlasticNaturalBCs.hpp>
150
151#ifdef ADD_CONTACT
152 #ifdef ENABLE_PYTHON_BINDING
153 #include <boost/python.hpp>
154 #include <boost/python/def.hpp>
155 #include <boost/python/numpy.hpp>
156namespace bp = boost::python;
157namespace np = boost::python::numpy;
158 #endif
159
160namespace ContactOps {
161
162double cn_contact = 0.1;
163
164}; // namespace ContactOps
165
166 #include <ContactOps.hpp>
167#endif // ADD_CONTACT
168
169using DomainRhsBCs = NaturalBC<DomainEleOp>::Assembly<AT>::LinearForm<IT>;
170using OpDomainRhsBCs =
171 DomainRhsBCs::OpFlux<PlasticOps::DomainBCs, 1, SPACE_DIM>;
172using BoundaryRhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::LinearForm<IT>;
173using OpBoundaryRhsBCs =
174 BoundaryRhsBCs::OpFlux<PlasticOps::BoundaryBCs, 1, SPACE_DIM>;
175using BoundaryLhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::BiLinearForm<IT>;
176using OpBoundaryLhsBCs =
177 BoundaryLhsBCs::OpFlux<PlasticOps::BoundaryBCs, 1, SPACE_DIM>;
178
179using namespace PlasticOps;
180using namespace HenckyOps;
181
182namespace PlasticOps {
183
184template <int FE_DIM, int PROBLEM_DIM, int SPACE_DIM> struct AddHOOps;
185
186template <> struct AddHOOps<2, 3, 3> {
187 AddHOOps() = delete;
188 static MoFEMErrorCode
189 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
190 std::vector<FieldSpace> space, std::string geom_field_name);
191};
192
193template <> struct AddHOOps<1, 2, 2> {
194 AddHOOps() = delete;
195 static MoFEMErrorCode
196 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
197 std::vector<FieldSpace> space, std::string geom_field_name);
198};
199
200template <> struct AddHOOps<3, 3, 3> {
201 AddHOOps() = delete;
202 static MoFEMErrorCode
203 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
204 std::vector<FieldSpace> space, std::string geom_field_name);
205};
206
207template <> struct AddHOOps<2, 2, 2> {
208 AddHOOps() = delete;
209 static MoFEMErrorCode
210 add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
211 std::vector<FieldSpace> space, std::string geom_field_name);
212};
213
214} // namespace PlasticOps
215
216struct Example {
217
218 Example(MoFEM::Interface &m_field) : mField(m_field) {}
219
220 MoFEMErrorCode runProblem();
221
222 enum { VOL, COUNT };
223 static inline std::array<double, 2> meshVolumeAndCount = {0, 0};
224
225private:
226 MoFEM::Interface &mField;
227
228 MoFEMErrorCode setupProblem();
229 MoFEMErrorCode createCommonData();
230 MoFEMErrorCode bC();
231 MoFEMErrorCode OPs();
232 MoFEMErrorCode tsSolve();
233 MoFEMErrorCode testOperators();
234
235 boost::shared_ptr<DomainEle> reactionFe;
236
237 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
238 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
239 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
240
241 struct ScaledTimeScale : public MoFEM::TimeScale {
242 using MoFEM::TimeScale::TimeScale;
243 double getScale(const double time) {
244 return scale * MoFEM::TimeScale::getScale(time);
245 };
246 };
247
248#ifdef ADD_CONTACT
249 #ifdef ENABLE_PYTHON_BINDING
250 boost::shared_ptr<ContactOps::SDFPython> sdfPythonPtr;
251 #endif
252#endif // ADD_CONTACT
253};
254
255//! [Run problem]
256MoFEMErrorCode Example::runProblem() {
257 MoFEMFunctionBegin;
258 CHKERR createCommonData();
259 CHKERR setupProblem();
260 CHKERR bC();
261 CHKERR OPs();
262 PetscBool test_ops = PETSC_FALSE;
263 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test_operators", &test_ops,
264 PETSC_NULLPTR);
265 if (test_ops == PETSC_FALSE) {
266 CHKERR tsSolve();
267 } else {
268 CHKERR testOperators();
269 }
270 MoFEMFunctionReturn(0);
271}
272//! [Run problem]
273
274//! [Set up problem]
275MoFEMErrorCode Example::setupProblem() {
276 MoFEMFunctionBegin;
277 Simple *simple = mField.getInterface<Simple>();
278
279 Range domain_ents;
280 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, domain_ents,
281 true);
282 auto get_ents_by_dim = [&](const auto dim) {
283 if (dim == SPACE_DIM) {
284 return domain_ents;
285 } else {
286 Range ents;
287 if (dim == 0)
288 CHKERR mField.get_moab().get_connectivity(domain_ents, ents, true);
289 else
290 CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents, true);
291 return ents;
292 }
293 };
294
295 auto get_base = [&]() {
296 auto domain_ents = get_ents_by_dim(SPACE_DIM);
297 if (domain_ents.empty())
298 CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Empty mesh");
299 const auto type = type_from_handle(domain_ents[0]);
300 switch (type) {
301 case MBQUAD:
302 return DEMKOWICZ_JACOBI_BASE;
303 case MBHEX:
304 return DEMKOWICZ_JACOBI_BASE;
305 case MBTRI:
306 return AINSWORTH_LEGENDRE_BASE;
307 case MBTET:
308 return AINSWORTH_LEGENDRE_BASE;
309 default:
310 CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Element type not handled");
311 }
312 return NOBASE;
313 };
314
315 const auto base = get_base();
316 MOFEM_LOG("PLASTICITY", Sev::inform)
317 << "Base " << ApproximationBaseNames[base];
318
319 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
320 CHKERR simple->addDomainField("EP", L2, base, size_symm);
321 CHKERR simple->addDomainField("TAU", L2, base, 1);
322 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
323
324 CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
325
326 PetscBool order_edge = PETSC_FALSE;
327 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-order_edge", &order_edge,
328 PETSC_NULLPTR);
329 PetscBool order_face = PETSC_FALSE;
330 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-order_face", &order_face,
331 PETSC_NULLPTR);
332 PetscBool order_volume = PETSC_FALSE;
333 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-order_volume", &order_volume,
334 PETSC_NULLPTR);
335
336 if (order_edge || order_face || order_volume) {
337
338 MOFEM_LOG("PLASTICITY", Sev::inform) << "Order edge " << order_edge
339 ? "true"
340 : "false";
341 MOFEM_LOG("PLASTICITY", Sev::inform) << "Order face " << order_face
342 ? "true"
343 : "false";
344 MOFEM_LOG("PLASTICITY", Sev::inform) << "Order volume " << order_volume
345 ? "true"
346 : "false";
347
348 auto ents = get_ents_by_dim(0);
349 if (order_edge)
350 ents.merge(get_ents_by_dim(1));
351 if (order_face)
352 ents.merge(get_ents_by_dim(2));
353 if (order_volume)
354 ents.merge(get_ents_by_dim(3));
355 CHKERR simple->setFieldOrder("U", order, &ents);
356 } else {
357 CHKERR simple->setFieldOrder("U", order);
358 }
359 CHKERR simple->setFieldOrder("EP", ep_order);
360 CHKERR simple->setFieldOrder("TAU", tau_order);
361
362 CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
363
364#ifdef ADD_CONTACT
365 CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
366 SPACE_DIM);
367 CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
368 SPACE_DIM);
369
370 auto get_skin = [&]() {
371 Range body_ents;
372 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
373 Skinner skin(&mField.get_moab());
374 Range skin_ents;
375 CHKERR skin.find_skin(0, body_ents, false, skin_ents);
376 return skin_ents;
377 };
378
379 auto filter_blocks = [&](auto skin) {
380 bool is_contact_block = true;
381 Range contact_range;
382 for (auto m :
383 mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
384
385 (boost::format("%s(.*)") % "CONTACT").str()
386
387 ))
388
389 ) {
390 is_contact_block =
391 true; ///< blocs interation is collective, so that is set irrespective
392 ///< if there are entities in given rank or not in the block
393 MOFEM_LOG("CONTACT", Sev::inform)
394 << "Find contact block set: " << m->getName();
395 auto meshset = m->getMeshset();
396 Range contact_meshset_range;
397 CHKERR mField.get_moab().get_entities_by_dimension(
398 meshset, SPACE_DIM - 1, contact_meshset_range, true);
399
400 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
401 contact_meshset_range);
402 contact_range.merge(contact_meshset_range);
403 }
404 if (is_contact_block) {
405 MOFEM_LOG("SYNC", Sev::inform)
406 << "Nb entities in contact surface: " << contact_range.size();
407 MOFEM_LOG_SYNCHRONISE(mField.get_comm());
408 skin = intersect(skin, contact_range);
409 }
410 return skin;
411 };
412
413 auto filter_true_skin = [&](auto skin) {
414 Range boundary_ents;
415 ParallelComm *pcomm =
416 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
417 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
418 PSTATUS_NOT, -1, &boundary_ents);
419 return boundary_ents;
420 };
421
422 auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
423 CHKERR simple->setFieldOrder("SIGMA", 0);
424 CHKERR simple->setFieldOrder("SIGMA", order - 1, &boundary_ents);
425#endif
426
427 CHKERR simple->setUp();
428 CHKERR simple->addFieldToEmptyFieldBlocks("U", "TAU");
429
430 auto project_ho_geometry = [&]() {
431 Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
432 return mField.loop_dofs("GEOMETRY", ent_method);
433 };
434 PetscBool project_geometry = PETSC_TRUE;
435 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-project_geometry",
436 &project_geometry, PETSC_NULLPTR);
437 if (project_geometry) {
438 CHKERR project_ho_geometry();
439 }
440
441 auto get_volume = [&]() {
442 using VolOp = DomainEle::UserDataOperator;
443 auto *op_ptr = new VolOp(NOSPACE, VolOp::OPSPACE);
444 std::array<double, 2> volume_and_count;
445 op_ptr->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
446 EntityType type,
447 EntitiesFieldData::EntData &data) {
448 MoFEMFunctionBegin;
449 auto op_ptr = static_cast<VolOp *>(base_op_ptr);
450 volume_and_count[VOL] += op_ptr->getMeasure();
451 volume_and_count[COUNT] += 1;
452 // in necessary at integration over Gauss points.
453 MoFEMFunctionReturn(0);
454 };
455 volume_and_count = {0, 0};
456 auto fe = boost::make_shared<DomainEle>(mField);
457 fe->getOpPtrVector().push_back(op_ptr);
458
459 auto dm = simple->getDM();
460 CHK_THROW_MESSAGE(
461 DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), fe),
462 "cac volume");
463 std::array<double, 2> tot_volume_and_count;
464 MPI_Allreduce(volume_and_count.data(), tot_volume_and_count.data(),
465 volume_and_count.size(), MPI_DOUBLE, MPI_SUM,
466 mField.get_comm());
467 return tot_volume_and_count;
468 };
469
470 meshVolumeAndCount = get_volume();
471 MOFEM_LOG("PLASTICITY", Sev::inform)
472 << "Mesh volume " << meshVolumeAndCount[VOL] << " nb. of elements "
473 << meshVolumeAndCount[COUNT];
474
475 MoFEMFunctionReturn(0);
476}
477//! [Set up problem]
478
479//! [Create common data]
480MoFEMErrorCode Example::createCommonData() {
481 MoFEMFunctionBegin;
482
483 auto get_command_line_parameters = [&]() {
484 MoFEMFunctionBegin;
485
486 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-scale", &scale, PETSC_NULLPTR);
487 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-young_modulus",
488 &young_modulus, PETSC_NULLPTR);
489 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-poisson_ratio",
490 &poisson_ratio, PETSC_NULLPTR);
491 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-hardening", &H, PETSC_NULLPTR);
492 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-hardening_viscous", &visH,
493 PETSC_NULLPTR);
494 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-yield_stress", &sigmaY,
495 PETSC_NULLPTR);
496 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-cn0", &cn0, PETSC_NULLPTR);
497 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-cn1", &cn1, PETSC_NULLPTR);
498 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-zeta", &zeta, PETSC_NULLPTR);
499 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-Qinf", &Qinf, PETSC_NULLPTR);
500 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-b_iso", &b_iso, PETSC_NULLPTR);
501 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-C1_k", &C1_k, PETSC_NULLPTR);
502 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-large_strains",
503 &is_large_strains, PETSC_NULLPTR);
504 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-set_timer", &set_timer,
505 PETSC_NULLPTR);
506 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-atom_test", &atom_test,
507 PETSC_NULLPTR);
508
509 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
510 PetscBool tau_order_is_set; ///< true if tau order is set
511 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-tau_order", &tau_order,
512 &tau_order_is_set);
513 PetscBool ep_order_is_set; ///< true if tau order is set
514 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-ep_order", &ep_order,
515 &ep_order_is_set);
516 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-geom_order", &geom_order,
517 PETSC_NULLPTR);
518
519 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-rho", &rho, PETSC_NULLPTR);
520 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-alpha_damping",
521 &alpha_damping, PETSC_NULLPTR);
522
523 MOFEM_LOG("PLASTICITY", Sev::inform) << "Young modulus " << young_modulus;
524 MOFEM_LOG("PLASTICITY", Sev::inform) << "Poisson ratio " << poisson_ratio;
525 MOFEM_LOG("PLASTICITY", Sev::inform) << "Yield stress " << sigmaY;
526 MOFEM_LOG("PLASTICITY", Sev::inform) << "Hardening " << H;
527 MOFEM_LOG("PLASTICITY", Sev::inform) << "Viscous hardening " << visH;
528 MOFEM_LOG("PLASTICITY", Sev::inform) << "Saturation yield stress " << Qinf;
529 MOFEM_LOG("PLASTICITY", Sev::inform) << "Saturation exponent " << b_iso;
530 MOFEM_LOG("PLASTICITY", Sev::inform) << "Kinematic hardening " << C1_k;
531 MOFEM_LOG("PLASTICITY", Sev::inform) << "cn0 " << cn0;
532 MOFEM_LOG("PLASTICITY", Sev::inform) << "cn1 " << cn1;
533 MOFEM_LOG("PLASTICITY", Sev::inform) << "zeta " << zeta;
534
535 if (tau_order_is_set == PETSC_FALSE)
536 tau_order = order - 2;
537 if (ep_order_is_set == PETSC_FALSE)
538 ep_order = order - 1;
539
540 MOFEM_LOG("PLASTICITY", Sev::inform) << "Approximation order " << order;
541 MOFEM_LOG("PLASTICITY", Sev::inform)
542 << "Ep approximation order " << ep_order;
543 MOFEM_LOG("PLASTICITY", Sev::inform)
544 << "Tau approximation order " << tau_order;
545 MOFEM_LOG("PLASTICITY", Sev::inform)
546 << "Geometry approximation order " << geom_order;
547
548 MOFEM_LOG("PLASTICITY", Sev::inform) << "Density " << rho;
549 MOFEM_LOG("PLASTICITY", Sev::inform) << "alpha_damping " << alpha_damping;
550
551 PetscBool is_scale = PETSC_TRUE;
552 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-is_scale", &is_scale,
553 PETSC_NULLPTR);
554 if (is_scale) {
555 scale /= young_modulus;
556 }
557
558 MOFEM_LOG("PLASTICITY", Sev::inform) << "Scale " << scale;
559
560#ifdef ADD_CONTACT
561 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-cn_contact",
562 &ContactOps::cn_contact, PETSC_NULLPTR);
563 MOFEM_LOG("CONTACT", Sev::inform)
564 << "cn_contact " << ContactOps::cn_contact;
565#endif // ADD_CONTACT
566
567 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-quasi_static",
568 &is_quasi_static, PETSC_NULLPTR);
569 MOFEM_LOG("PLASTICITY", Sev::inform)
570 << "Is quasi static: " << (is_quasi_static ? "true" : "false");
571
572 MoFEMFunctionReturn(0);
573 };
574
575 CHKERR get_command_line_parameters();
576
577#ifdef ADD_CONTACT
578 #ifdef ENABLE_PYTHON_BINDING
579 auto file_exists = [](std::string myfile) {
580 std::ifstream file(myfile.c_str());
581 if (file) {
582 return true;
583 }
584 return false;
585 };
586 char sdf_file_name[255] = "sdf.py";
587 CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, "-sdf_file",
588 sdf_file_name, 255, PETSC_NULLPTR);
589
590 if (file_exists(sdf_file_name)) {
591 MOFEM_LOG("CONTACT", Sev::inform) << sdf_file_name << " file found";
592 sdfPythonPtr = boost::make_shared<ContactOps::SDFPython>();
593 CHKERR sdfPythonPtr->sdfInit(sdf_file_name);
594 ContactOps::sdfPythonWeakPtr = sdfPythonPtr;
595 } else {
596 MOFEM_LOG("CONTACT", Sev::warning) << sdf_file_name << " file NOT found";
597 }
598 #endif
599#endif // ADD_CONTACT
600
601 MoFEMFunctionReturn(0);
602}
603//! [Create common data]
604
605//! [Boundary condition]
606MoFEMErrorCode Example::bC() {
607 MoFEMFunctionBegin;
608
609 auto simple = mField.getInterface<Simple>();
610 auto bc_mng = mField.getInterface<BcManager>();
611
612 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
613 "U", 0, 0);
614 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
615 "U", 1, 1);
616 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
617 "U", 2, 2);
618 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
619 "REMOVE_ALL", "U", 0, 3);
620
621#ifdef ADD_CONTACT
622 for (auto b : {"FIX_X", "REMOVE_X"})
623 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
624 "SIGMA", 0, 0, false, true);
625 for (auto b : {"FIX_Y", "REMOVE_Y"})
626 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
627 "SIGMA", 1, 1, false, true);
628 for (auto b : {"FIX_Z", "REMOVE_Z"})
629 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
630 "SIGMA", 2, 2, false, true);
631 for (auto b : {"FIX_ALL", "REMOVE_ALL"})
632 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
633 "SIGMA", 0, 3, false, true);
634 CHKERR bc_mng->removeBlockDOFsOnEntities(
635 simple->getProblemName(), "NO_CONTACT", "SIGMA", 0, 3, false, true);
636#endif
637
638 CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
639 simple->getProblemName(), "U");
640
641 auto &bc_map = bc_mng->getBcMapByBlockName();
642 for (auto bc : bc_map)
643 MOFEM_LOG("PLASTICITY", Sev::verbose) << "Marker " << bc.first;
644
645 MoFEMFunctionReturn(0);
646}
647//! [Boundary condition]
648
649//! [Push operators to pipeline]
650MoFEMErrorCode Example::OPs() {
651 MoFEMFunctionBegin;
652 auto pip_mng = mField.getInterface<PipelineManager>();
653
654 auto integration_rule_bc = [](int, int, int ao) { return 2 * ao; };
655
656 auto vol_rule = [](int, int, int ao) { return 2 * ao + geom_order - 1; };
657
658 auto add_boundary_ops_lhs_mechanical = [&](auto &pip) {
659 MoFEMFunctionBegin;
660
661 CHKERR PlasticOps::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
662 pip, {HDIV}, "GEOMETRY");
663 pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
664
665 // Add Natural BCs to LHS
666 CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
667 pip, mField, "U", Sev::inform);
668
669#ifdef ADD_CONTACT
670 auto simple = mField.getInterface<Simple>();
671 CHKERR
672 ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
673 pip, "SIGMA", "U");
674 CHKERR
675 ContactOps::opFactoryBoundaryToDomainLhs<SPACE_DIM, AT, IT, DomainEle>(
676 mField, pip, simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",
677 vol_rule);
678#endif // ADD_CONTACT
679
680 MoFEMFunctionReturn(0);
681 };
682
683 auto add_boundary_ops_rhs_mechanical = [&](auto &pip) {
684 MoFEMFunctionBegin;
685
686 CHKERR PlasticOps::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
687 pip, {HDIV}, "GEOMETRY");
688 pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
689
690 // Add Natural BCs to RHS
691 CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
692 pip, mField, "U", {boost::make_shared<ScaledTimeScale>()}, Sev::inform);
693
694#ifdef ADD_CONTACT
695 CHKERR ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, IT, BoundaryEleOp>(
696 pip, "SIGMA", "U");
697#endif // ADD_CONTACT
698
699 MoFEMFunctionReturn(0);
700 };
701
702 auto add_domain_ops_lhs = [this](auto &pip) {
703 MoFEMFunctionBegin;
704 CHKERR PlasticOps::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
705 pip, {H1, HDIV}, "GEOMETRY");
706
707 if (is_quasi_static == PETSC_FALSE) {
708
709 //! [Only used for dynamics]
710 using OpMass = FormsIntegrators<DomainEleOp>::Assembly<AT>::BiLinearForm<
711 GAUSS>::OpMass<1, SPACE_DIM>;
712 //! [Only used for dynamics]
713
714 auto get_inertia_and_mass_damping = [this](const double, const double,
715 const double) {
716 auto *pip = mField.getInterface<PipelineManager>();
717 auto &fe_domain_lhs = pip->getDomainLhsFE();
718 return (rho / scale) * fe_domain_lhs->ts_aa +
719 (alpha_damping / scale) * fe_domain_lhs->ts_a;
720 };
721 pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
722 }
723
724 CHKERR PlasticOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
725 mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
726
727 MoFEMFunctionReturn(0);
728 };
729
730 auto add_domain_ops_rhs = [this](auto &pip) {
731 MoFEMFunctionBegin;
732
733 CHKERR PlasticOps::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
734 pip, {H1, HDIV}, "GEOMETRY");
735
736 CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
737 pip, mField, "U",
738 {boost::make_shared<ScaledTimeScale>("body_force_hist.txt")},
739 Sev::inform);
740
741 // only in case of dynamics
742 if (is_quasi_static == PETSC_FALSE) {
743
744 //! [Only used for dynamics]
745 using OpInertiaForce = FormsIntegrators<DomainEleOp>::Assembly<
746 AT>::LinearForm<IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
747 //! [Only used for dynamics]
748
749 auto mat_acceleration = boost::make_shared<MatrixDouble>();
750 pip.push_back(new OpCalculateVectorFieldValuesDotDot<SPACE_DIM>(
751 "U", mat_acceleration));
752 pip.push_back(
753 new OpInertiaForce("U", mat_acceleration, [](double, double, double) {
754 return rho / scale;
755 }));
756 if (alpha_damping > 0) {
757 auto mat_velocity = boost::make_shared<MatrixDouble>();
758 pip.push_back(
759 new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", mat_velocity));
760 pip.push_back(
761 new OpInertiaForce("U", mat_velocity, [](double, double, double) {
762 return alpha_damping / scale;
763 }));
764 }
765 }
766
767 CHKERR PlasticOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
768 mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
769
770#ifdef ADD_CONTACT
771 CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
772 pip, "SIGMA", "U");
773#endif // ADD_CONTACT
774
775 MoFEMFunctionReturn(0);
776 };
777
778 CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
779 CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
780
781 // Boundary
782 CHKERR add_boundary_ops_lhs_mechanical(pip_mng->getOpBoundaryLhsPipeline());
783 CHKERR add_boundary_ops_rhs_mechanical(pip_mng->getOpBoundaryRhsPipeline());
784
785 CHKERR pip_mng->setDomainRhsIntegrationRule(vol_rule);
786 CHKERR pip_mng->setDomainLhsIntegrationRule(vol_rule);
787
788 CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_bc);
789 CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_bc);
790
791 auto create_reaction_pipeline = [&](auto &pip) {
792 MoFEMFunctionBegin;
793 CHKERR PlasticOps::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
794 pip, {H1}, "GEOMETRY");
795 CHKERR PlasticOps::opFactoryDomainReactions<SPACE_DIM, AT, IT, DomainEleOp>(
796 mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
797 MoFEMFunctionReturn(0);
798 };
799
800 reactionFe = boost::make_shared<DomainEle>(mField);
801 reactionFe->getRuleHook = vol_rule;
802 CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
803 reactionFe->postProcessHook =
804 EssentialPreProcReaction<DisplacementCubitBcData>(mField, reactionFe);
805
806 MoFEMFunctionReturn(0);
807}
808//! [Push operators to pipeline]
809
810//! [Solve]
811struct SetUpSchur {
812
813 /**
814 * @brief Create data structure for handling Schur complement
815 *
816 * @param m_field
817 * @param sub_dm Schur complement sub dm
818 * @param field_split_it IS of Schur block
819 * @param ao_map AO map from sub dm to main problem
820 * @return boost::shared_ptr<SetUpSchur>
821 */
822 static boost::shared_ptr<SetUpSchur> createSetUpSchur(
823
824 MoFEM::Interface &m_field, SmartPetscObj<DM> sub_dm,
825 SmartPetscObj<IS> field_split_it, SmartPetscObj<AO> ao_map
826
827 );
828 virtual MoFEMErrorCode setUp(TS solver) = 0;
829
830protected:
831 SetUpSchur() = default;
832};
833
834MoFEMErrorCode Example::tsSolve() {
835 MoFEMFunctionBegin;
836
837 Simple *simple = mField.getInterface<Simple>();
838 PipelineManager *pip_mng = mField.getInterface<PipelineManager>();
839 ISManager *is_manager = mField.getInterface<ISManager>();
840
841 auto snes_ctx_ptr = getDMSnesCtx(simple->getDM());
842
843 auto set_section_monitor = [&](auto solver) {
844 MoFEMFunctionBegin;
845 SNES snes;
846 CHKERR TSGetSNES(solver, &snes);
847 CHKERR SNESMonitorSet(snes,
848 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
849 void *))MoFEMSNESMonitorFields,
850 (void *)(snes_ctx_ptr.get()), nullptr);
851 MoFEMFunctionReturn(0);
852 };
853
854 auto create_post_process_elements = [&]() {
855 auto push_vol_ops = [this](auto &pip) {
856 CHKERR PlasticOps::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
857 pip, {H1, HDIV}, "GEOMETRY");
858
859 auto [common_plastic_ptr, common_hencky_ptr] =
860 PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
861 mField, "MAT_PLASTIC", pip, "U", "EP", "TAU", 1., Sev::inform);
862
863 if (common_hencky_ptr) {
864 if (common_plastic_ptr->mGradPtr != common_hencky_ptr->matGradPtr)
865 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong pointer for grad");
866 }
867
868 return std::make_pair(common_plastic_ptr, common_hencky_ptr);
869 };
870
871 auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
872 MoFEMFunctionBegin;
873
874 auto &pip = pp_fe->getOpPtrVector();
875
876 auto [common_plastic_ptr, common_hencky_ptr] = p;
877
878 using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
879
880 auto x_ptr = boost::make_shared<MatrixDouble>();
881 pip.push_back(
882 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
883 auto u_ptr = boost::make_shared<MatrixDouble>();
884 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
885
886 if (is_large_strains) {
887
888 pip.push_back(
889
890 new OpPPMap(
891
892 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
893
894 {{"PLASTIC_SURFACE",
895 common_plastic_ptr->getPlasticSurfacePtr()},
896 {"PLASTIC_MULTIPLIER",
897 common_plastic_ptr->getPlasticTauPtr()}},
898
899 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
900
901 {{"GRAD", common_hencky_ptr->matGradPtr},
902 {"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()}},
903
904 {{"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()},
905 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
906 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
907
908 )
909
910 );
911
912 } else {
913
914 pip.push_back(
915
916 new OpPPMap(
917
918 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
919
920 {{"PLASTIC_SURFACE",
921 common_plastic_ptr->getPlasticSurfacePtr()},
922 {"PLASTIC_MULTIPLIER",
923 common_plastic_ptr->getPlasticTauPtr()}},
924
925 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
926
927 {},
928
929 {{"STRAIN", common_plastic_ptr->mStrainPtr},
930 {"STRESS", common_plastic_ptr->mStressPtr},
931 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
932 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
933
934 )
935
936 );
937 }
938
939 MoFEMFunctionReturn(0);
940 };
941
942 PetscBool post_proc_vol;
943 PetscBool post_proc_skin;
944
945 if constexpr (SPACE_DIM == 2) {
946 post_proc_vol = PETSC_TRUE;
947 post_proc_skin = PETSC_FALSE;
948 } else {
949 post_proc_vol = PETSC_FALSE;
950 post_proc_skin = PETSC_TRUE;
951 }
952 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_vol", &post_proc_vol,
953 PETSC_NULLPTR);
954 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_skin",
955 &post_proc_skin, PETSC_NULLPTR);
956
957 auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops,
958 post_proc_vol]() {
959 if (post_proc_vol == PETSC_FALSE)
960 return boost::shared_ptr<PostProcEle>();
961 auto pp_fe = boost::make_shared<PostProcEle>(mField);
962 CHK_MOAB_THROW(
963 push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
964 "push_vol_post_proc_ops");
965 return pp_fe;
966 };
967
968 auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops,
969 post_proc_skin]() {
970 if (post_proc_skin == PETSC_FALSE)
971 return boost::shared_ptr<SkinPostProcEle>();
972
973 auto simple = mField.getInterface<Simple>();
974 auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
975 auto op_side = new OpLoopSide<SideEle>(mField, simple->getDomainFEName(),
976 SPACE_DIM, Sev::verbose);
977 pp_fe->getOpPtrVector().push_back(op_side);
978 CHK_MOAB_THROW(push_vol_post_proc_ops(
979 pp_fe, push_vol_ops(op_side->getOpPtrVector())),
980 "push_vol_post_proc_ops");
981 return pp_fe;
982 };
983
984 return std::make_pair(vol_post_proc(), skin_post_proc());
985 };
986
987 auto scatter_create = [&](auto D, auto coeff) {
988 SmartPetscObj<IS> is;
989 CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
990 ROW, "U", coeff, coeff, is);
991 int loc_size;
992 CHKERR ISGetLocalSize(is, &loc_size);
993 Vec v;
994 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
995 VecScatter scatter;
996 CHKERR VecScatterCreate(D, is, v, PETSC_NULLPTR, &scatter);
997 return std::make_tuple(SmartPetscObj<Vec>(v),
998 SmartPetscObj<VecScatter>(scatter));
999 };
1000
1001 boost::shared_ptr<SetPtsData> field_eval_data;
1002 boost::shared_ptr<MatrixDouble> u_field_ptr;
1003
1004 std::array<double, 3> field_eval_coords{0.0, 0.0, 0.0};
1005 int coords_dim = 3;
1006 CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
1007 field_eval_coords.data(), &coords_dim,
1008 &do_eval_field);
1009
1010 boost::shared_ptr<std::map<std::string, boost::shared_ptr<VectorDouble>>>
1011 scalar_field_ptrs = boost::make_shared<
1012 std::map<std::string, boost::shared_ptr<VectorDouble>>>();
1013 boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
1014 vector_field_ptrs = boost::make_shared<
1015 std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
1016 boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
1017 sym_tensor_field_ptrs = boost::make_shared<
1018 std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
1019 boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
1020 tensor_field_ptrs = boost::make_shared<
1021 std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
1022
1023 if (do_eval_field) {
1024 auto u_field_ptr = boost::make_shared<MatrixDouble>();
1025 field_eval_data =
1026 mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
1027
1028 CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree<SPACE_DIM>(
1029 field_eval_data, simple->getDomainFEName());
1030
1031 field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
1032 auto no_rule = [](int, int, int) { return -1; };
1033 auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
1034 field_eval_fe_ptr->getRuleHook = no_rule;
1035
1036 CHKERR PlasticOps::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
1037 field_eval_fe_ptr->getOpPtrVector(), {H1, HDIV}, "GEOMETRY");
1038
1039 auto [common_plastic_ptr, common_hencky_ptr] =
1040 PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
1041 mField, "MAT_PLASTIC", field_eval_fe_ptr->getOpPtrVector(), "U",
1042 "EP", "TAU", 1., Sev::inform);
1043
1044 field_eval_fe_ptr->getOpPtrVector().push_back(
1045 new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_field_ptr));
1046
1047 if ((common_plastic_ptr) && (common_hencky_ptr) && (scalar_field_ptrs)) {
1048 if (is_large_strains) {
1049 scalar_field_ptrs->insert(
1050 {"PLASTIC_SURFACE", common_plastic_ptr->getPlasticSurfacePtr()});
1051 scalar_field_ptrs->insert(
1052 {"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()});
1053 vector_field_ptrs->insert({"U", u_field_ptr});
1054 sym_tensor_field_ptrs->insert(
1055 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()});
1056 sym_tensor_field_ptrs->insert(
1057 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()});
1058 sym_tensor_field_ptrs->insert(
1059 {"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()});
1060 tensor_field_ptrs->insert({"GRAD", common_hencky_ptr->matGradPtr});
1061 tensor_field_ptrs->insert(
1062 {"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()});
1063 } else {
1064 scalar_field_ptrs->insert(
1065 {"PLASTIC_SURFACE", common_plastic_ptr->getPlasticSurfacePtr()});
1066 scalar_field_ptrs->insert(
1067 {"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()});
1068 vector_field_ptrs->insert({"U", u_field_ptr});
1069 sym_tensor_field_ptrs->insert(
1070 {"STRAIN", common_plastic_ptr->mStrainPtr});
1071 sym_tensor_field_ptrs->insert(
1072 {"STRESS", common_plastic_ptr->mStressPtr});
1073 sym_tensor_field_ptrs->insert(
1074 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()});
1075 sym_tensor_field_ptrs->insert(
1076 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()});
1077 }
1078 }
1079 }
1080
1081 auto test_monitor_ptr = boost::make_shared<FEMethod>();
1082
1083 auto set_time_monitor = [&](auto dm, auto solver) {
1084 MoFEMFunctionBegin;
1085 boost::shared_ptr<Monitor<SPACE_DIM>> monitor_ptr(new Monitor<SPACE_DIM>(
1086 dm, create_post_process_elements(), reactionFe, uXScatter, uYScatter,
1087 uZScatter, field_eval_coords, field_eval_data, scalar_field_ptrs,
1088 vector_field_ptrs, sym_tensor_field_ptrs, tensor_field_ptrs));
1089 boost::shared_ptr<ForcesAndSourcesCore> null;
1090
1091 test_monitor_ptr->postProcessHook = [&]() {
1092 MoFEMFunctionBegin;
1093
1094 if (atom_test && fabs(test_monitor_ptr->ts_t - 0.5) < 1e-12 &&
1095 test_monitor_ptr->ts_step == 25) {
1096
1097 if (scalar_field_ptrs->at("PLASTIC_MULTIPLIER")->size()) {
1098 auto t_tau =
1099 getFTensor0FromVec(*scalar_field_ptrs->at("PLASTIC_MULTIPLIER"));
1100 MOFEM_LOG("PlasticSync", Sev::inform) << "Eval point tau: " << t_tau;
1101
1102 if (atom_test == 1 && fabs(t_tau - 0.688861) > 1e-5) {
1103 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1104 "atom test %d failed: wrong plastic multiplier value",
1105 atom_test);
1106 }
1107 }
1108
1109 if (vector_field_ptrs->at("U")->size1()) {
1110 FTensor::Index<'i', SPACE_DIM> i;
1111 auto t_disp =
1112 getFTensor1FromMat<SPACE_DIM>(*vector_field_ptrs->at("U"));
1113 MOFEM_LOG("PlasticSync", Sev::inform) << "Eval point U: " << t_disp;
1114
1115 if (atom_test == 1 && fabs(t_disp(0) - 0.25 / 2.) > 1e-5 ||
1116 fabs(t_disp(1) + 0.0526736) > 1e-5) {
1117 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1118 "atom test %d failed: wrong displacement value",
1119 atom_test);
1120 }
1121 }
1122
1123 if (sym_tensor_field_ptrs->at("PLASTIC_STRAIN")->size1()) {
1124 auto t_plastic_strain = getFTensor2SymmetricFromMat<SPACE_DIM>(
1125 *sym_tensor_field_ptrs->at("PLASTIC_STRAIN"));
1126 MOFEM_LOG("PlasticSync", Sev::inform)
1127 << "Eval point EP: " << t_plastic_strain;
1128
1129 if (atom_test == 1 &&
1130 fabs(t_plastic_strain(0, 0) - 0.221943) > 1e-5 ||
1131 fabs(t_plastic_strain(0, 1)) > 1e-5 ||
1132 fabs(t_plastic_strain(1, 1) + 0.110971) > 1e-5) {
1133 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1134 "atom test %d failed: wrong plastic strain value",
1135 atom_test);
1136 }
1137 }
1138
1139 if (tensor_field_ptrs->at("FIRST_PIOLA")->size1()) {
1140 auto t_piola_stress = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
1141 *tensor_field_ptrs->at("FIRST_PIOLA"));
1142 MOFEM_LOG("PlasticSync", Sev::inform)
1143 << "Eval point Piola stress: " << t_piola_stress;
1144
1145 if (atom_test == 1 && fabs((t_piola_stress(0, 0) - 198.775) /
1146 t_piola_stress(0, 0)) > 1e-5 ||
1147 fabs(t_piola_stress(0, 1)) + fabs(t_piola_stress(1, 0)) +
1148 fabs(t_piola_stress(1, 1)) >
1149 1e-5) {
1150 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1151 "atom test %d failed: wrong Piola stress value",
1152 atom_test);
1153 }
1154 }
1155 }
1156
1157 MOFEM_LOG_SYNCHRONISE(mField.get_comm());
1158 MoFEMFunctionReturn(0);
1159 };
1160
1161 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
1162 monitor_ptr, null, test_monitor_ptr);
1163
1164 MoFEMFunctionReturn(0);
1165 };
1166
1167 auto set_schur_pc = [&](auto solver,
1168 boost::shared_ptr<SetUpSchur> &schur_ptr) {
1169 MoFEMFunctionBeginHot;
1170
1171 auto name_prb = simple->getProblemName();
1172
1173 // create sub dm for Schur complement
1174 auto create_schur_dm = [&](SmartPetscObj<DM> base_dm,
1175 SmartPetscObj<DM> &dm_sub) {
1176 MoFEMFunctionBegin;
1177 dm_sub = createDM(mField.get_comm(), "DMMOFEM");
1178 CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "SCHUR");
1179 CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
1180 CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
1181 CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
1182 for (auto f : {"U"}) {
1183 CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
1184 CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
1185 }
1186 CHKERR DMSetUp(dm_sub);
1187
1188 MoFEMFunctionReturn(0);
1189 };
1190
1191 auto create_block_dm = [&](SmartPetscObj<DM> base_dm,
1192 SmartPetscObj<DM> &dm_sub) {
1193 MoFEMFunctionBegin;
1194 dm_sub = createDM(mField.get_comm(), "DMMOFEM");
1195 CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "BLOCK");
1196 CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
1197 CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
1198 CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
1199#ifdef ADD_CONTACT
1200 for (auto f : {"SIGMA", "EP", "TAU"}) {
1201 CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
1202 CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
1203 }
1204#else
1205 for (auto f : {"EP", "TAU"}) {
1206 CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
1207 CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
1208 }
1209#endif
1210 CHKERR DMSetUp(dm_sub);
1211 MoFEMFunctionReturn(0);
1212 };
1213
1214 // Create nested (sub BC) Schur DM
1215 if constexpr (AT == AssemblyType::BLOCK_SCHUR) {
1216
1217 SmartPetscObj<DM> dm_schur;
1218 CHKERR create_schur_dm(simple->getDM(), dm_schur);
1219 SmartPetscObj<DM> dm_block;
1220 CHKERR create_block_dm(simple->getDM(), dm_block);
1221
1222#ifdef ADD_CONTACT
1223
1224 auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
1225 auto block_mat_data = createBlockMatStructure(
1226 simple->getDM(),
1227
1228 {
1229
1230 {simple->getDomainFEName(),
1231
1232 {{"U", "U"},
1233 {"SIGMA", "SIGMA"},
1234 {"U", "SIGMA"},
1235 {"SIGMA", "U"},
1236 {"EP", "EP"},
1237 {"TAU", "TAU"},
1238 {"U", "EP"},
1239 {"EP", "U"},
1240 {"EP", "TAU"},
1241 {"TAU", "EP"},
1242 {"TAU", "U"}
1243
1244 }},
1245
1246 {simple->getBoundaryFEName(),
1247
1248 {{"SIGMA", "SIGMA"}, {"U", "SIGMA"}, {"SIGMA", "U"}
1249
1250 }}
1251
1252 }
1253
1254 );
1255
1256 return createSchurNestedMatrixStruture(
1257
1258 {dm_schur, dm_block}, block_mat_data,
1259
1260 {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, true
1261
1262 );
1263 };
1264
1265#else
1266
1267 auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
1268 auto block_mat_data =
1269 createBlockMatStructure(simple->getDM(),
1270
1271 {{simple->getDomainFEName(),
1272
1273 {{"U", "U"},
1274 {"EP", "EP"},
1275 {"TAU", "TAU"},
1276 {"U", "EP"},
1277 {"EP", "U"},
1278 {"EP", "TAU"},
1279 {"TAU", "U"},
1280 {"TAU", "EP"}
1281
1282 }}}
1283
1284 );
1285
1286 return createSchurNestedMatrixStruture(
1287
1288 {dm_schur, dm_block}, block_mat_data,
1289
1290 {"EP", "TAU"}, {nullptr, nullptr}, false
1291
1292 );
1293 };
1294
1295#endif
1296
1297 auto nested_mat_data = get_nested_mat_data(dm_schur, dm_block);
1298 CHKERR DMMoFEMSetNestSchurData(simple->getDM(), nested_mat_data);
1299
1300 auto block_is = getDMSubData(dm_block)->getSmartRowIs();
1301 auto ao_schur = getDMSubData(dm_schur)->getSmartRowMap();
1302
1303 // Indices has to be map fro very to level, while assembling Schur
1304 // complement.
1305 schur_ptr =
1306 SetUpSchur::createSetUpSchur(mField, dm_schur, block_is, ao_schur);
1307 CHKERR schur_ptr->setUp(solver);
1308 }
1309
1310 MoFEMFunctionReturnHot(0);
1311 };
1312
1313 auto dm = simple->getDM();
1314 auto D = createDMVector(dm);
1315 auto DD = vectorDuplicate(D);
1316 uXScatter = scatter_create(D, 0);
1317 uYScatter = scatter_create(D, 1);
1318 if constexpr (SPACE_DIM == 3)
1319 uZScatter = scatter_create(D, 2);
1320
1321 auto create_solver = [pip_mng]() {
1322 if (is_quasi_static == PETSC_TRUE)
1323 return pip_mng->createTSIM();
1324 else
1325 return pip_mng->createTSIM2();
1326 };
1327
1328 auto solver = create_solver();
1329
1330 auto active_pre_lhs = []() {
1331 MoFEMFunctionBegin;
1332 std::fill(PlasticOps::CommonData::activityData.begin(),
1333 PlasticOps::CommonData::activityData.end(), 0);
1334 MoFEMFunctionReturn(0);
1335 };
1336
1337 auto active_post_lhs = [&]() {
1338 MoFEMFunctionBegin;
1339 auto get_iter = [&]() {
1340 SNES snes;
1341 CHK_THROW_MESSAGE(TSGetSNES(solver, &snes), "Can not get SNES");
1342 int iter;
1343 CHK_THROW_MESSAGE(SNESGetIterationNumber(snes, &iter),
1344 "Can not get iter");
1345 return iter;
1346 };
1347
1348 auto iter = get_iter();
1349 if (iter >= 0) {
1350
1351 std::array<int, 5> activity_data;
1352 std::fill(activity_data.begin(), activity_data.end(), 0);
1353 MPI_Allreduce(PlasticOps::CommonData::activityData.data(),
1354 activity_data.data(), activity_data.size(), MPI_INT,
1355 MPI_SUM, mField.get_comm());
1356
1357 int &active_points = activity_data[0];
1358 int &avtive_full_elems = activity_data[1];
1359 int &avtive_elems = activity_data[2];
1360 int &nb_points = activity_data[3];
1361 int &nb_elements = activity_data[4];
1362
1363 if (nb_points) {
1364
1365 double proc_nb_points =
1366 100 * static_cast<double>(active_points) / nb_points;
1367 double proc_nb_active =
1368 100 * static_cast<double>(avtive_elems) / nb_elements;
1369 double proc_nb_full_active = 100;
1370 if (avtive_elems)
1371 proc_nb_full_active =
1372 100 * static_cast<double>(avtive_full_elems) / avtive_elems;
1373
1374 MOFEM_LOG_C("PLASTICITY", Sev::inform,
1375 "Iter %d nb pts %d nb active pts %d (%3.3f\%) nb active "
1376 "elements %d "
1377 "(%3.3f\%) nb full active elems %d (%3.3f\%)",
1378 iter, nb_points, active_points, proc_nb_points,
1379 avtive_elems, proc_nb_active, avtive_full_elems,
1380 proc_nb_full_active, iter);
1381 }
1382 }
1383
1384 MoFEMFunctionReturn(0);
1385 };
1386
1387 auto add_active_dofs_elem = [&](auto dm) {
1388 MoFEMFunctionBegin;
1389 auto fe_pre_proc = boost::make_shared<FEMethod>();
1390 fe_pre_proc->preProcessHook = active_pre_lhs;
1391 auto fe_post_proc = boost::make_shared<FEMethod>();
1392 fe_post_proc->postProcessHook = active_post_lhs;
1393 auto ts_ctx_ptr = getDMTsCtx(dm);
1394 ts_ctx_ptr->getPreProcessIJacobian().push_front(fe_pre_proc);
1395 ts_ctx_ptr->getPostProcessIJacobian().push_back(fe_post_proc);
1396 MoFEMFunctionReturn(0);
1397 };
1398
1399 auto set_essential_bc = [&](auto dm, auto solver) {
1400 MoFEMFunctionBegin;
1401 // This is low level pushing finite elements (pipelines) to solver
1402
1403 auto pre_proc_ptr = boost::make_shared<FEMethod>();
1404 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
1405 auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
1406
1407 // Add boundary condition scaling
1408 auto disp_time_scale = boost::make_shared<TimeScale>();
1409
1410 auto get_bc_hook_rhs = [this, pre_proc_ptr, disp_time_scale]() {
1411 EssentialPreProc<DisplacementCubitBcData> hook(mField, pre_proc_ptr,
1412 {disp_time_scale}, false);
1413 return hook;
1414 };
1415 pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
1416
1417 auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
1418 MoFEMFunctionBegin;
1419 CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
1420 mField, post_proc_rhs_ptr, nullptr, Sev::verbose)();
1421 CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(
1422 mField, post_proc_rhs_ptr, 1.)();
1423 MoFEMFunctionReturn(0);
1424 };
1425 auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
1426 return EssentialPostProcLhs<DisplacementCubitBcData>(
1427 mField, post_proc_lhs_ptr, 1.);
1428 };
1429 post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
1430
1431 auto ts_ctx_ptr = getDMTsCtx(dm);
1432 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
1433 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
1434 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
1435 post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
1436 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
1437
1438 MoFEMFunctionReturn(0);
1439 };
1440
1441 auto B = createDMMatrix(dm);
1442 if (is_quasi_static == PETSC_FALSE) {
1443 CHKERR TSSetIJacobian(solver, B, B, PETSC_NULLPTR, PETSC_NULLPTR);
1444 } else {
1445 CHKERR TSSetI2Jacobian(solver, B, B, PETSC_NULLPTR, PETSC_NULLPTR);
1446 }
1447 if (is_quasi_static == PETSC_TRUE) {
1448 CHKERR TSSetSolution(solver, D);
1449 } else {
1450 CHKERR TS2SetSolution(solver, D, DD);
1451 }
1452 CHKERR set_section_monitor(solver);
1453 CHKERR set_time_monitor(dm, solver);
1454 CHKERR TSSetFromOptions(solver);
1455
1456 CHKERR add_active_dofs_elem(dm);
1457 boost::shared_ptr<SetUpSchur> schur_ptr;
1458 CHKERR set_schur_pc(solver, schur_ptr);
1459 CHKERR set_essential_bc(dm, solver);
1460
1461 MOFEM_LOG_CHANNEL("TIMER");
1462 MOFEM_LOG_TAG("TIMER", "timer");
1463 if (set_timer)
1464 BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
1465 MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp";
1466 CHKERR TSSetUp(solver);
1467 MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp <= done";
1468 MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve";
1469 CHKERR TSSolve(solver, NULL);
1470 MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve <= done";
1471
1472 MoFEMFunctionReturn(0);
1473}
1474//! [Solve]
1475
1476//! [TestOperators]
1477MoFEMErrorCode Example::testOperators() {
1478 MoFEMFunctionBegin;
1479
1480 // get operators tester
1481 auto simple = mField.getInterface<Simple>();
1482 auto opt = mField.getInterface<OperatorsTester>(); // get interface to
1483 // OperatorsTester
1484 auto pip = mField.getInterface<PipelineManager>(); // get interface to
1485 // pipeline manager
1486
1487 constexpr double eps = 1e-9;
1488
1489 auto x = opt->setRandomFields(simple->getDM(), {
1490
1491 {"U", {-1e-6, 1e-6}},
1492
1493 {"EP", {-1e-6, 1e-6}},
1494
1495 {"TAU", {0, 1e-4}}
1496
1497 });
1498
1499 auto dot_x_plastic_active =
1500 opt->setRandomFields(simple->getDM(), {
1501
1502 {"U", {-1, 1}},
1503
1504 {"EP", {-1, 1}},
1505
1506 {"TAU", {0.1, 1}}
1507
1508 });
1509 auto diff_x_plastic_active =
1510 opt->setRandomFields(simple->getDM(), {
1511
1512 {"U", {-1, 1}},
1513
1514 {"EP", {-1, 1}},
1515
1516 {"TAU", {-1, 1}}
1517
1518 });
1519
1520 auto dot_x_elastic =
1521 opt->setRandomFields(simple->getDM(), {
1522
1523 {"U", {-1, 1}},
1524
1525 {"EP", {-1, 1}},
1526
1527 {"TAU", {-1, -0.1}}
1528
1529 });
1530 auto diff_x_elastic =
1531 opt->setRandomFields(simple->getDM(), {
1532
1533 {"U", {-1, 1}},
1534
1535 {"EP", {-1, 1}},
1536
1537 {"TAU", {-1, 1}}
1538
1539 });
1540
1541 auto test_domain_ops = [&](auto fe_name, auto lhs_pipeline, auto rhs_pipeline,
1542 auto dot_x, auto diff_x) {
1543 MoFEMFunctionBegin;
1544
1545 auto diff_res = opt->checkCentralFiniteDifference(
1546 simple->getDM(), fe_name, rhs_pipeline, lhs_pipeline, x, dot_x,
1547 SmartPetscObj<Vec>(), diff_x, 0, 0.5, eps);
1548
1549 // Calculate norm of difference between directional derivative calculated
1550 // from finite difference, and tangent matrix.
1551 double fnorm;
1552 CHKERR VecNorm(diff_res, NORM_2, &fnorm);
1553 MOFEM_LOG_C("PLASTICITY", Sev::inform,
1554 "Test consistency of tangent matrix %3.4e", fnorm);
1555
1556 constexpr double err = 1e-5;
1557 if (fnorm > err)
1558 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1559 "Norm of directional derivative too large err = %3.4e", fnorm);
1560
1561 MoFEMFunctionReturn(0);
1562 };
1563
1564 MOFEM_LOG("PLASTICITY", Sev::inform) << "Elastic active";
1565 CHKERR test_domain_ops(simple->getDomainFEName(), pip->getDomainLhsFE(),
1566 pip->getDomainRhsFE(), dot_x_elastic, diff_x_elastic);
1567
1568 MOFEM_LOG("PLASTICITY", Sev::inform) << "Plastic active";
1569 CHKERR test_domain_ops(simple->getDomainFEName(), pip->getDomainLhsFE(),
1570 pip->getDomainRhsFE(), dot_x_plastic_active,
1571 diff_x_plastic_active);
1572
1573 MoFEMFunctionReturn(0);
1574};
1575
1576//! [TestOperators]
1577
1578static char help[] = "...\n\n";
1579
1580int main(int argc, char *argv[]) {
1581
1582#ifdef ADD_CONTACT
1583 #ifdef ENABLE_PYTHON_BINDING
1584 Py_Initialize();
1585 np::initialize();
1586 #endif
1587#endif // ADD_CONTACT
1588
1589 // Initialisation of MoFEM/PETSc and MOAB data structures
1590 const char param_file[] = "param_file.petsc";
1591 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
1592
1593 // Add logging channel for example
1594 auto core_log = logging::core::get();
1595 core_log->add_sink(
1596 LogManager::createSink(LogManager::getStrmWorld(), "PLASTICITY"));
1597 core_log->add_sink(
1598 LogManager::createSink(LogManager::getStrmWorld(), "TIMER"));
1599 LogManager::setLog("PLASTICITY");
1600 MOFEM_LOG_TAG("PLASTICITY", "Plasticity");
1601
1602#ifdef ADD_CONTACT
1603 core_log->add_sink(
1604 LogManager::createSink(LogManager::getStrmWorld(), "CONTACT"));
1605 LogManager::setLog("CONTACT");
1606 MOFEM_LOG_TAG("CONTACT", "Contact");
1607#endif // ADD_CONTACT
1608
1609 core_log->add_sink(
1610 LogManager::createSink(LogManager::getStrmSync(), "PlasticSync"));
1611 LogManager::setLog("PlasticSync");
1612 MOFEM_LOG_TAG("PlasticSync", "PlasticSync");
1613
1614 try {
1615
1616 //! [Register MoFEM discrete manager in PETSc]
1617 DMType dm_name = "DMMOFEM";
1618 CHKERR DMRegister_MoFEM(dm_name);
1619 //! [Register MoFEM discrete manager in PETSc
1620
1621 //! [Create MoAB]
1622 moab::Core mb_instance; ///< mesh database
1623 moab::Interface &moab = mb_instance; ///< mesh database interface
1624 //! [Create MoAB]
1625
1626 //! [Create MoFEM]
1627 MoFEM::Core core(moab); ///< finite element database
1628 MoFEM::Interface &m_field = core; ///< finite element database interface
1629 //! [Create MoFEM]
1630
1631 //! [Load mesh]
1632 Simple *simple = m_field.getInterface<Simple>();
1633 CHKERR simple->getOptions();
1634 CHKERR simple->loadFile();
1635 //! [Load mesh]
1636
1637 //! [Example]
1638 Example ex(m_field);
1639 CHKERR ex.runProblem();
1640 //! [Example]
1641 }
1642 CATCH_ERRORS;
1643
1644 CHKERR MoFEM::Core::Finalize();
1645
1646#ifdef ADD_CONTACT
1647 #ifdef ENABLE_PYTHON_BINDING
1648 if (Py_FinalizeEx() < 0) {
1649 exit(120);
1650 }
1651 #endif
1652#endif // ADD_CONTACT
1653
1654 return 0;
1655}
1656
1657struct SetUpSchurImpl : public SetUpSchur {
1658
1659 SetUpSchurImpl(MoFEM::Interface &m_field, SmartPetscObj<DM> sub_dm,
1660 SmartPetscObj<IS> field_split_is, SmartPetscObj<AO> ao_up)
1661 : SetUpSchur(), mField(m_field), subDM(sub_dm),
1662 fieldSplitIS(field_split_is), aoSchur(ao_up) {
1663 if (S) {
1664 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
1665 "Is expected that schur matrix is not "
1666 "allocated. This is "
1667 "possible only is PC is set up twice");
1668 }
1669 }
1670 virtual ~SetUpSchurImpl() { S.reset(); }
1671
1672 MoFEMErrorCode setUp(TS solver);
1673 MoFEMErrorCode preProc();
1674 MoFEMErrorCode postProc();
1675
1676private:
1677 SmartPetscObj<Mat> S;
1678
1679 MoFEM::Interface &mField;
1680 SmartPetscObj<DM> subDM; ///< field split sub dm
1681 SmartPetscObj<IS> fieldSplitIS; ///< IS for split Schur block
1682 SmartPetscObj<AO> aoSchur; ///> main DM to subDM
1683};
1684
1685MoFEMErrorCode SetUpSchurImpl::setUp(TS solver) {
1686 MoFEMFunctionBegin;
1687 auto simple = mField.getInterface<Simple>();
1688 auto pip_mng = mField.getInterface<PipelineManager>();
1689
1690 SNES snes;
1691 CHKERR TSGetSNES(solver, &snes);
1692 KSP ksp;
1693 CHKERR SNESGetKSP(snes, &ksp);
1694 CHKERR KSPSetFromOptions(ksp);
1695
1696 PC pc;
1697 CHKERR KSPGetPC(ksp, &pc);
1698 PetscBool is_pcfs = PETSC_FALSE;
1699 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1700 if (is_pcfs) {
1701 if (S) {
1702 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
1703 "Is expected that schur matrix is not "
1704 "allocated. This is "
1705 "possible only is PC is set up twice");
1706 }
1707
1708 S = createDMMatrix(subDM);
1709 CHKERR MatSetBlockSize(S, SPACE_DIM);
1710
1711 // Set DM to use shell block matrix
1712 DM solver_dm;
1713 CHKERR TSGetDM(solver, &solver_dm);
1714 CHKERR DMSetMatType(solver_dm, MATSHELL);
1715
1716 auto ts_ctx_ptr = getDMTsCtx(solver_dm);
1717 auto A = createDMBlockMat(simple->getDM());
1718 auto P = createDMNestSchurMat(simple->getDM());
1719
1720 if (is_quasi_static == PETSC_TRUE) {
1721 auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a,
1722 Mat A, Mat B, void *ctx) {
1723 return TsSetIJacobian(ts, t, u, u_t, a, B, A, ctx);
1724 };
1725 CHKERR TSSetIJacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
1726 } else {
1727 auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, Vec utt,
1728 PetscReal a, PetscReal aa, Mat A, Mat B,
1729 void *ctx) {
1730 return TsSetI2Jacobian(ts, t, u, u_t, utt, a, aa, B, A, ctx);
1731 };
1732 CHKERR TSSetI2Jacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
1733 }
1734 CHKERR KSPSetOperators(ksp, A, P);
1735
1736 auto set_ops = [&]() {
1737 MoFEMFunctionBegin;
1738 auto pip_mng = mField.getInterface<PipelineManager>();
1739
1740#ifndef ADD_CONTACT
1741 // Boundary
1742 pip_mng->getOpBoundaryLhsPipeline().push_front(
1743 createOpSchurAssembleBegin());
1744 pip_mng->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd(
1745
1746 {"EP", "TAU"}, {nullptr, nullptr}, aoSchur, S, false, false
1747
1748 ));
1749 // Domain
1750 pip_mng->getOpDomainLhsPipeline().push_front(
1751 createOpSchurAssembleBegin());
1752 pip_mng->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd(
1753
1754 {"EP", "TAU"}, {nullptr, nullptr}, aoSchur, S, false, false
1755
1756 ));
1757#else
1758
1759 double eps_stab = 1e-4;
1760 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-eps_stab", &eps_stab,
1761 PETSC_NULLPTR);
1762
1763 using B = FormsIntegrators<BoundaryEleOp>::Assembly<
1764 BLOCK_PRECONDITIONER_SCHUR>::BiLinearForm<IT>;
1765 using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
1766
1767 // Boundary
1768 pip_mng->getOpBoundaryLhsPipeline().push_front(
1769 createOpSchurAssembleBegin());
1770 pip_mng->getOpBoundaryLhsPipeline().push_back(
1771 new OpMassStab("SIGMA", "SIGMA", [eps_stab](double, double, double) {
1772 return eps_stab;
1773 }));
1774 pip_mng->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd(
1775
1776 {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, aoSchur, S,
1777 false, false
1778
1779 ));
1780 // Domain
1781 pip_mng->getOpDomainLhsPipeline().push_front(
1782 createOpSchurAssembleBegin());
1783 pip_mng->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd(
1784
1785 {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, aoSchur, S,
1786 false, false
1787
1788 ));
1789#endif // ADD_CONTACT
1790 MoFEMFunctionReturn(0);
1791 };
1792
1793 auto set_assemble_elems = [&]() {
1794 MoFEMFunctionBegin;
1795 auto schur_asmb_pre_proc = boost::make_shared<FEMethod>();
1796 schur_asmb_pre_proc->preProcessHook = [this]() {
1797 MoFEMFunctionBegin;
1798 CHKERR MatZeroEntries(S);
1799 MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble Begin";
1800 MoFEMFunctionReturn(0);
1801 };
1802 auto schur_asmb_post_proc = boost::make_shared<FEMethod>();
1803
1804 schur_asmb_post_proc->postProcessHook = [this, schur_asmb_post_proc]() {
1805 MoFEMFunctionBegin;
1806 MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble End";
1807
1808 // Apply essential constrains to Schur complement
1809 CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
1810 CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
1811 CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
1812 mField, schur_asmb_post_proc, 1, S, aoSchur)();
1813
1814 MoFEMFunctionReturn(0);
1815 };
1816 auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
1817 ts_ctx_ptr->getPreProcessIJacobian().push_front(schur_asmb_pre_proc);
1818 ts_ctx_ptr->getPostProcessIJacobian().push_front(schur_asmb_post_proc);
1819 MoFEMFunctionReturn(0);
1820 };
1821
1822 auto set_pc = [&]() {
1823 MoFEMFunctionBegin;
1824 CHKERR PCFieldSplitSetIS(pc, NULL, fieldSplitIS);
1825 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1826 MoFEMFunctionReturn(0);
1827 };
1828
1829 auto set_diagonal_pc = [&]() {
1830 MoFEMFunctionBegin;
1831 KSP *subksp;
1832 CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULLPTR, &subksp);
1833 auto get_pc = [](auto ksp) {
1834 PC pc_raw;
1835 CHKERR KSPGetPC(ksp, &pc_raw);
1836 return SmartPetscObj<PC>(pc_raw,
1837 true); // bump reference
1838 };
1839 CHKERR setSchurA00MatSolvePC(get_pc(subksp[0]));
1840 CHKERR PetscFree(subksp);
1841 MoFEMFunctionReturn(0);
1842 };
1843
1844 CHKERR set_ops();
1845 CHKERR set_pc();
1846 CHKERR set_assemble_elems();
1847
1848 CHKERR TSSetUp(solver);
1849 CHKERR KSPSetUp(ksp);
1850 CHKERR set_diagonal_pc();
1851
1852 } else {
1853 pip_mng->getOpBoundaryLhsPipeline().push_front(
1854 createOpSchurAssembleBegin());
1855 pip_mng->getOpBoundaryLhsPipeline().push_back(
1856 createOpSchurAssembleEnd({}, {}));
1857 pip_mng->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
1858 pip_mng->getOpDomainLhsPipeline().push_back(
1859 createOpSchurAssembleEnd({}, {}));
1860 }
1861
1862 // fieldSplitIS.reset();
1863 // aoSchur.reset();
1864 MoFEMFunctionReturn(0);
1865}
1866
1867boost::shared_ptr<SetUpSchur>
1868SetUpSchur::createSetUpSchur(MoFEM::Interface &m_field,
1869 SmartPetscObj<DM> sub_dm, SmartPetscObj<IS> is_sub,
1870 SmartPetscObj<AO> ao_up) {
1871 return boost::shared_ptr<SetUpSchur>(
1872 new SetUpSchurImpl(m_field, sub_dm, is_sub, ao_up));
1873}
1874
1875namespace PlasticOps {
1876
1877MoFEMErrorCode AddHOOps<2, 3, 3>::add(
1878 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
1879 std::vector<FieldSpace> spaces, std::string geom_field_name) {
1880 MoFEMFunctionBegin;
1881 CHKERR MoFEM::AddHOOps<2, 3, 3>::add(pipeline, spaces, geom_field_name);
1882 MoFEMFunctionReturn(0);
1883}
1884
1885MoFEMErrorCode AddHOOps<1, 2, 2>::add(
1886 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
1887 std::vector<FieldSpace> spaces, std::string geom_field_name) {
1888 MoFEMFunctionBegin;
1889 CHKERR MoFEM::AddHOOps<1, 2, 2>::add(pipeline, spaces, geom_field_name);
1890 MoFEMFunctionReturn(0);
1891}
1892
1893template <int FE_DIM, int PROBLEM_DIM, int SPACE_DIM>
1894MoFEMErrorCode
1895scaleL2(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
1896 std::string geom_field_name) {
1897 MoFEMFunctionBegin;
1898
1899 auto jac_ptr = boost::make_shared<MatrixDouble>();
1900 auto det_ptr = boost::make_shared<VectorDouble>();
1901 pipeline.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
1902 geom_field_name, jac_ptr));
1903 pipeline.push_back(new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, nullptr));
1904
1905 auto scale_ptr = boost::make_shared<double>(1.);
1906 auto scale = Example::meshVolumeAndCount[0] /
1907 Example::meshVolumeAndCount[1]; // average volume of elements
1908 using OP = ForcesAndSourcesCore::UserDataOperator;
1909 auto op_scale = new OP(NOSPACE, OP::OPSPACE);
1910 op_scale->doWorkRhsHook = [scale_ptr, det_ptr,
1911 scale](DataOperator *base_op_ptr, int, EntityType,
1912 EntitiesFieldData::EntData &) {
1913 *scale_ptr = scale / det_ptr->size(); // distribute average element size
1914 // over integration points
1915 return 0;
1916 };
1917 pipeline.push_back(op_scale);
1918
1919 for (auto base : {AINSWORTH_LEGENDRE_BASE, AINSWORTH_LOBATTO_BASE,
1920 DEMKOWICZ_JACOBI_BASE}) {
1921 pipeline.push_back(
1922 new OpScaleBaseBySpaceInverseOfMeasure(L2, base, det_ptr, scale_ptr));
1923 }
1924
1925 MoFEMFunctionReturn(0);
1926}
1927
1928MoFEMErrorCode AddHOOps<3, 3, 3>::add(
1929 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
1930 std::vector<FieldSpace> spaces, std::string geom_field_name) {
1931 MoFEMFunctionBegin;
1932 constexpr bool scale_l2 = false;
1933 if (scale_l2) {
1934 CHKERR scaleL2<3, 3, 3>(pipeline, geom_field_name);
1935 }
1936 CHKERR MoFEM::AddHOOps<3, 3, 3>::add(pipeline, spaces, geom_field_name,
1937 nullptr, nullptr, nullptr);
1938 MoFEMFunctionReturn(0);
1939}
1940
1941MoFEMErrorCode AddHOOps<2, 2, 2>::add(
1942 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
1943 std::vector<FieldSpace> spaces, std::string geom_field_name) {
1944 MoFEMFunctionBegin;
1945 constexpr bool scale_l2 = false;
1946 if (scale_l2) {
1947 CHKERR scaleL2<2, 2, 2>(pipeline, geom_field_name);
1948 }
1949 CHKERR MoFEM::AddHOOps<2, 2, 2>::add(pipeline, spaces, geom_field_name,
1950 nullptr, nullptr, nullptr);
1951 MoFEMFunctionReturn(0);
1952}
1953
1954} // namespace PlasticOps
filter_true_skin
static auto filter_true_skin(MoFEM::Interface &m_field, Range &&skin)
Definition EshelbianPlasticity.cpp:166
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition LogManager.hpp:353
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition LogManager.hpp:319
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
SPACE_DIM
constexpr int SPACE_DIM
[Define dimension]
Definition adjoint.cpp:22
poisson_ratio
constexpr double poisson_ratio
Definition adjoint.cpp:74
young_modulus
constexpr double young_modulus
Definition adjoint.cpp:73
main
int main()
Definition adol-c_atom.cpp:46
a
constexpr double a
Definition approx_sphere.cpp:30
eps
static const double eps
Definition check_base_functions_derivatives_on_tet.cpp:11
SPACE_DIM
constexpr int SPACE_DIM
Definition child_and_parent.cpp:16
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition child_and_parent.cpp:39
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
[Define dimension]
Definition child_and_parent.cpp:34
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition Kronecker_Delta.hpp:49
FTensor::Tensor2_symmetric
Definition Tensor2_symmetric_value.hpp:14
SCHUR_ASSEMBLE
#define SCHUR_ASSEMBLE
Definition contact.cpp:18
ROW
@ ROW
Definition definitions.h:136
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition definitions.h:385
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition definitions.h:62
NOBASE
@ NOBASE
Definition definitions.h:59
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition definitions.h:612
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
FieldSpace
FieldSpace
approximation spaces
Definition definitions.h:82
L2
@ L2
field with C-1 continuity
Definition definitions.h:88
H1
@ H1
continuous field
Definition definitions.h:85
NOSPACE
@ NOSPACE
Definition definitions.h:83
HCURL
@ HCURL
field with continuous tangents
Definition definitions.h:86
HDIV
@ HDIV
field with continuous normal traction
Definition definitions.h:87
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition definitions.h:228
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition definitions.h:592
MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition definitions.h:33
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
ApproximationBaseNames
static const char *const ApproximationBaseNames[]
Definition definitions.h:72
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
order
constexpr int order
Definition dg_projection.cpp:18
OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
Definition dynamic_first_order_con_law.cpp:61
t_kd
constexpr auto t_kd
Definition free_surface.cpp:139
MoFEM::DMMoFEMCreateSubDM
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:215
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:488
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition DMMoFEM.cpp:450
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition DMMoFEM.cpp:238
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
MoFEM::DMoFEMLoopFiniteElements
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:576
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1102
MoFEM::DMMoFEMAddSubFieldCol
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition DMMoFEM.cpp:280
MoFEM::createDMMatrix
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:1059
MoFEM::IntegrationType
IntegrationType
Form integrator integration types.
Definition FormsIntegrators.hpp:136
MoFEM::AssemblyType
AssemblyType
[Storage and set boundary conditions]
Definition FormsIntegrators.hpp:104
MoFEM::GAUSS
@ GAUSS
Definition FormsIntegrators.hpp:136
MoFEM::PETSC
@ PETSC
Definition FormsIntegrators.hpp:105
MoFEM::BLOCK_PRECONDITIONER_SCHUR
@ BLOCK_PRECONDITIONER_SCHUR
Definition FormsIntegrators.hpp:109
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition LogManager.cpp:392
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition LogManager.hpp:347
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition LogManager.hpp:292
MoFEM::CoreInterface::loop_dofs
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.
MoFEM::MeshsetsManager::getCubitMeshsetPtr
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
Definition MeshsetsManager.cpp:589
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
D
double D
Definition initial_diffusion.cpp:44
v
const double v
phase velocity of light in medium (cm/ns)
Definition initial_diffusion.cpp:40
l
FTensor::Index< 'l', 3 > l
Definition matrix_function.cpp:21
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:20
ContactOps::opFactoryDomainRhs
MoFEMErrorCode opFactoryDomainRhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
Definition ContactOps.hpp:1184
ContactOps::cn_contact
double cn_contact
Definition contact.cpp:99
ContactOps::opFactoryBoundaryToDomainLhs
MoFEMErrorCode opFactoryBoundaryToDomainLhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string fe_domain_name, std::string sigma, std::string u, std::string geom, ForcesAndSourcesCore::RuleHookFun rule, bool is_axisymmetric=false)
Definition ContactOps.hpp:1263
ContactOps::opFactoryBoundaryRhs
MoFEMErrorCode opFactoryBoundaryRhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
Definition ContactOps.hpp:1341
ContactOps::opFactoryBoundaryLhs
MoFEMErrorCode opFactoryBoundaryLhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
Definition ContactOps.hpp:1317
HenckyOps::eps
constexpr double eps
Definition HenckyOps.hpp:13
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEM::type_from_handle
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition Templates.hpp:1932
MoFEM::DMMoFEMTSSetMonitor
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:1046
MoFEM::TsSetIJacobian
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:169
MoFEM::getDMTsCtx
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition DMMoFEM.hpp:1144
MoFEM::createOpSchurAssembleEnd
OpSchurAssembleBase * createOpSchurAssembleEnd(std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range > > field_ents, SmartPetscObj< AO > ao, SmartPetscObj< Mat > schur, bool sym_schur, bool symm_op)
Construct a new Op Schur Assemble End object.
Definition Schur.cpp:2585
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition DeprecatedPetsc.hpp:142
MoFEM::setSchurA00MatSolvePC
MoFEMErrorCode setSchurA00MatSolvePC(SmartPetscObj< PC > pc)
Set PC for A00 block.
Definition Schur.cpp:2627
MoFEM::MoFEMSNESMonitorFields
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *ctx)
Sens monitor printing residual field by field.
Definition SnesCtx.cpp:592
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition DeprecatedPetsc.hpp:182
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition DeprecatedPetsc.hpp:162
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition PetscSmartObj.hpp:223
MoFEM::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition DeprecatedPetsc.hpp:192
MoFEM::getDMSubData
auto getDMSubData(DM dm)
Get sub problem data structure.
Definition DMMoFEM.hpp:1160
MoFEM::getFTensor1FromMat
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
Definition Templates.hpp:238
MoFEM::getFTensor2FromMat
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
Definition Templates.hpp:278
MoFEM::TsSetI2Jacobian
PetscErrorCode TsSetI2Jacobian(TS ts, PetscReal t, Vec u, Vec u_t, Vec u_tt, PetscReal a, PetscReal aa, Mat A, Mat B, void *ctx)
Calculation Jacobian for second order PDE in time.
Definition TsCtx.cpp:519
MoFEM::createBlockMatStructure
boost::shared_ptr< BlockStructure > createBlockMatStructure(DM dm, SchurFEOpsFEandFields schur_fe_op_vec)
Create a Mat Diag Blocks object.
Definition Schur.cpp:1082
MoFEM::createSchurNestedMatrixStruture
boost::shared_ptr< NestSchurData > createSchurNestedMatrixStruture(std::pair< SmartPetscObj< DM >, SmartPetscObj< DM > > dms, boost::shared_ptr< BlockStructure > block_mat_data_ptr, std::vector< std::string > fields_names, std::vector< boost::shared_ptr< Range > > field_ents, bool add_preconditioner_block)
Get the Schur Nest Mat Array object.
Definition Schur.cpp:2343
MoFEM::getFTensor2SymmetricFromMat
static auto getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
Definition Templates.hpp:334
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition Templates.hpp:135
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition DeprecatedPetsc.hpp:172
MoFEM::DMMoFEMSetNestSchurData
MoFEMErrorCode DMMoFEMSetNestSchurData(DM dm, boost::shared_ptr< NestSchurData >)
Definition DMMoFEM.cpp:1554
MoFEM::getDMSnesCtx
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition DMMoFEM.hpp:1130
MoFEM::createDMNestSchurMat
auto createDMNestSchurMat(DM dm)
Definition DMMoFEM.hpp:1086
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition PetscSmartObj.hpp:143
MoFEM::createOpSchurAssembleBegin
OpSchurAssembleBase * createOpSchurAssembleBegin()
Definition Schur.cpp:2580
MoFEM::createDMBlockMat
auto createDMBlockMat(DM dm)
Definition DMMoFEM.hpp:1079
PlasticOps::createCommonPlasticOps
auto createCommonPlasticOps(MoFEM::Interface &m_field, std::string block_name, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string u, std::string ep, std::string tau, double scale, Sev sev)
Definition PlasticOps.hpp:426
PlasticOps::opFactoryDomainReactions
MoFEMErrorCode opFactoryDomainReactions(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
Definition PlasticOps.hpp:630
PlasticOps::opFactoryDomainRhs
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
Definition PlasticOps.hpp:500
PlasticOps::scaleL2
MoFEMErrorCode scaleL2(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string geom_field_name)
Definition plastic.cpp:1895
PlasticOps::opFactoryDomainLhs
MoFEMErrorCode opFactoryDomainLhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
Definition PlasticOps.hpp:547
A
constexpr AssemblyType A
Definition operators_tests.cpp:30
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition photon_diffusion.cpp:29
young_modulus
double young_modulus
Young modulus.
Definition plastic.cpp:125
AT
constexpr AssemblyType AT
Definition plastic.cpp:44
C1_k
double C1_k
Kinematic hardening.
Definition plastic.cpp:133
Qinf
double Qinf
Saturation yield stress.
Definition plastic.cpp:131
IT
constexpr IntegrationType IT
Definition plastic.cpp:47
help
static char help[]
[TestOperators]
Definition plastic.cpp:1578
rho
double rho
Definition plastic.cpp:144
atom_test
int atom_test
Atom test.
Definition plastic.cpp:121
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition plastic.cpp:13
do_eval_field
PetscBool do_eval_field
Evaluate field.
Definition plastic.cpp:119
is_quasi_static
PetscBool is_quasi_static
Definition plastic.cpp:143
alpha_damping
double alpha_damping
Definition plastic.cpp:145
SPACE_DIM
constexpr int SPACE_DIM
Definition plastic.cpp:40
visH
double visH
Viscous hardening.
Definition plastic.cpp:129
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition plastic.cpp:126
kinematic_hardening
auto kinematic_hardening(FTensor::Tensor2_symmetric< T, DIM > &t_plastic_strain, double C1_k)
Definition plastic.cpp:92
set_timer
PetscBool set_timer
Set timer.
Definition plastic.cpp:118
iso_hardening_dtau
double iso_hardening_dtau(double tau, double H, double Qinf, double b_iso)
Definition plastic.cpp:78
scale
double scale
Definition plastic.cpp:123
size_symm
constexpr auto size_symm
Definition plastic.cpp:42
zeta
double zeta
Viscous hardening.
Definition plastic.cpp:130
H
double H
Hardening.
Definition plastic.cpp:128
tau_order
int tau_order
Order of tau files.
Definition plastic.cpp:139
iso_hardening_exp
double iso_hardening_exp(double tau, double b_iso)
Definition plastic.cpp:64
cn0
double cn0
Definition plastic.cpp:135
order
int order
Order displacement.
Definition plastic.cpp:138
b_iso
double b_iso
Saturation exponent.
Definition plastic.cpp:132
is_large_strains
PetscBool is_large_strains
Large strains.
Definition plastic.cpp:117
geom_order
int geom_order
Order if fixed.
Definition plastic.cpp:141
sigmaY
double sigmaY
Yield stress.
Definition plastic.cpp:127
iso_hardening
double iso_hardening(double tau, double H, double Qinf, double b_iso, double sigmaY)
Definition plastic.cpp:73
kinematic_hardening_dplastic_strain
auto kinematic_hardening_dplastic_strain(double C1_k)
Definition plastic.cpp:106
SideEle
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition plastic.cpp:61
ep_order
int ep_order
Order of ep files.
Definition plastic.cpp:140
cn1
double cn1
Definition plastic.cpp:136
CONTACT_SPACE
constexpr FieldSpace CONTACT_SPACE
Definition plastic.cpp:52
t
constexpr double t
plate stiffness
Definition plate.cpp:58
OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition seepage.cpp:55
m
FTensor::Index< 'm', 3 > m
Definition shallow_wave.cpp:80
ElementsAndOps< 2 >::SideEle
PipelineManager::ElementsAndOpsByDim< 2 >::FaceSideEle SideEle
Definition plastic.cpp:29
ElementsAndOps< 3 >::SideEle
PipelineManager::ElementsAndOpsByDim< 3 >::FaceSideEle SideEle
Definition plastic.cpp:36
Example::ScaledTimeScale::getScale
double getScale(const double time)
Get scaling at given time.
Definition plastic.cpp:243
Example
[Example]
Definition plastic.cpp:216
Example::uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition plastic.cpp:239
Example::uXScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition plastic.cpp:237
Example::testOperators
MoFEMErrorCode testOperators()
[Solve]
Definition plastic.cpp:1477
Example::tsSolve
MoFEMErrorCode tsSolve()
Definition plastic.cpp:834
Example::base
FieldApproximationBase base
Definition plot_base.cpp:68
Example::createCommonData
MoFEMErrorCode createCommonData()
[Set up problem]
Definition plastic.cpp:480
Example::meshVolumeAndCount
static std::array< double, 2 > meshVolumeAndCount
Definition plastic.cpp:223
Example::uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition plastic.cpp:238
Example::Example
Example(MoFEM::Interface &m_field)
Definition plastic.cpp:218
Example::OPs
MoFEMErrorCode OPs()
[Boundary condition]
Definition plastic.cpp:650
Example::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition plastic.cpp:256
Example::mField
MoFEM::Interface & mField
Definition plastic.cpp:226
Example::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition plastic.cpp:275
Example::bC
MoFEMErrorCode bC()
[Create common data]
Definition plastic.cpp:606
Example::reactionFe
boost::shared_ptr< DomainEle > reactionFe
Definition plastic.cpp:235
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition HODataOperators.hpp:417
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition BcManager.hpp:29
MoFEM::CommInterface
Managing BitRefLevels.
Definition CommInterface.hpp:21
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition Core.hpp:82
MoFEM::CoreTmp< 0 >::Initialize
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
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:118
MoFEM::DataOperator
base operator to do operations at Gauss Pt. level
Definition DataOperators.hpp:24
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition BCData.hpp:72
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition EntitiesFieldData.hpp:128
MoFEM::EssentialPostProcLhs
Class (Function) to enforce essential constrains on the left hand side diagonal.
Definition Essential.hpp:33
MoFEM::EssentialPostProcRhs
Class (Function) to enforce essential constrains on the right hand side diagonal.
Definition Essential.hpp:41
MoFEM::EssentialPreProcReaction
Class (Function) to calculate residual side diagonal.
Definition Essential.hpp:49
MoFEM::EssentialPreProc
Class (Function) to enforce essential constrains.
Definition Essential.hpp:25
MoFEM::FieldEvaluatorInterface
Field evaluator interface.
Definition FieldEvaluator.hpp:21
MoFEM::ForcesAndSourcesCore::UserDataOperator::getMeasure
double getMeasure() const
get measure of element
Definition ForcesAndSourcesCore.hpp:1289
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPSPACE
@ OPSPACE
operator do Work is execute on space data
Definition ForcesAndSourcesCore.hpp:567
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition ISManager.hpp:23
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition LogManager.cpp:293
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition LogManager.cpp:339
MoFEM::LogManager::getStrmSync
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
Definition LogManager.cpp:343
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition MeshsetsManager.hpp:104
MoFEM::NaturalBC::Assembly
Assembly methods.
Definition Natural.hpp:65
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Definition UserDataOperators.hpp:1576
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition UserDataOperators.hpp:598
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor field rank 1, i.e. vector field.
Definition UserDataOperators.hpp:470
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition ForcesAndSourcesCore.hpp:1304
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition PostProcBrokenMeshInMoabBase.hpp:709
MoFEM::OpScaleBaseBySpaceInverseOfMeasure
Scale base functions by inverses of measure of element.
Definition HODataOperators.hpp:391
MoFEM::OperatorsTester
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
Definition OperatorsTester.hpp:21
MoFEM::PipelineManager
PipelineManager interface.
Definition PipelineManager.hpp:24
MoFEM::PipelineManager::VolEle
MoFEM::VolumeElementForcesAndSourcesCore VolEle
Definition PipelineManager.hpp:34
MoFEM::PipelineManager::getDomainLhsFE
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Definition PipelineManager.hpp:451
MoFEM::PipelineManager::FaceEle
MoFEM::FaceElementForcesAndSourcesCore FaceEle
Definition PipelineManager.hpp:35
MoFEM::PipelineManager::EdgeEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Definition PipelineManager.hpp:36
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition Projection10NodeCoordsOnField.hpp:24
MoFEM::Simple
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
MoFEM::SmartPetscObj
intrusive_ptr for managing petsc objects
Definition PetscSmartObj.hpp:85
MoFEM::TimeScale
Force scale operator for reading two columns.
Definition ScalingMethod.hpp:32
MoFEM::TimeScale::getScale
double getScale(const double time)
Get scaling at a given time.
Definition ScalingMethod.cpp:137
MoFEM::TimeScale::TimeScale
TimeScale(std::string file_name="", bool error_if_file_not_given=false, ScalingFun def_scaling_fun=[](double time) { return time;})
TimeScale constructor.
Definition ScalingMethod.cpp:22
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93
PlasticOps::CommonData::activityData
static std::array< int, 5 > activityData
Definition PlasticOps.hpp:107
SetUpSchurImpl::subDM
SmartPetscObj< DM > subDM
field split sub dm
Definition plastic.cpp:1680
SetUpSchurImpl::S
SmartPetscObj< Mat > S
Definition test_broken_space.cpp:529
SetUpSchurImpl::aoSchur
SmartPetscObj< AO > aoSchur
Definition plastic.cpp:1682
SetUpSchurImpl::SetUpSchurImpl
SetUpSchurImpl(MoFEM::Interface &m_field, SmartPetscObj< DM > sub_dm, SmartPetscObj< IS > field_split_is, SmartPetscObj< AO > ao_up)
Definition plastic.cpp:1659
SetUpSchurImpl::setUp
MoFEMErrorCode setUp(SmartPetscObj< KSP >)
Definition test_broken_space.cpp:532
SetUpSchurImpl::fieldSplitIS
SmartPetscObj< IS > fieldSplitIS
IS for split Schur block.
Definition plastic.cpp:1681
SetUpSchurImpl::~SetUpSchurImpl
virtual ~SetUpSchurImpl()
Definition plastic.cpp:1670
SetUpSchurImpl::postProc
MoFEMErrorCode postProc()
SetUpSchurImpl::preProc
MoFEMErrorCode preProc()
SetUpSchurImpl::mField
MoFEM::Interface & mField
Definition test_broken_space.cpp:528
SetUpSchur
[Push operators to pipeline]
Definition test_broken_space.cpp:44
SetUpSchur::SetUpSchur
SetUpSchur()=default
SetUpSchur::createSetUpSchur
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
Definition test_broken_space.cpp:828
SetUpSchur::setUp
virtual MoFEMErrorCode setUp(TS solver)=0
AT
constexpr AssemblyType AT
Definition test_broken_space.cpp:24
VolOp
VolEle::UserDataOperator VolOp
Definition test_cache_on_entities.cpp:23
Generated by Doxygen 1.12.0 and hosted at University of Glasgow