v0.13.2
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#ifndef EXECUTABLE_DIMENSION
10#define EXECUTABLE_DIMENSION 3
11#endif
12
13#include <MoFEM.hpp>
14#include <MatrixFunction.hpp>
15#include <IntegrationRules.hpp>
16
17using namespace MoFEM;
18
19template <int DIM> struct ElementsAndOps {};
20
21template <> struct ElementsAndOps<2> {
24};
25
26template <> struct ElementsAndOps<3> {
29};
30
31constexpr int SPACE_DIM =
32 EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
33constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
34
35constexpr AssemblyType A = (SCHUR_ASSEMBLE)
36 ? AssemblyType::SCHUR
37 : AssemblyType::PETSC; //< selected assembly type
38constexpr IntegrationType G =
39 IntegrationType::GAUSS; //< selected integration type
40
47
49
50//! [Only used with Hooke equation (linear material model)]
52 GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
55//! [Only used with Hooke equation (linear material model)]
56
57//! [Only used for dynamics]
62//! [Only used for dynamics]
63
64//! [Only used with Hencky/nonlinear material]
66 GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
69//! [Only used with Hencky/nonlinear material]
70
71//! [Essential boundary conditions]
78//! [Essential boundary conditions]
80
84
86using OpForce =
88
93
94PetscBool is_large_strains = PETSC_TRUE;
95PetscBool set_timer = PETSC_FALSE;
96
97double scale = 1.;
98
99double young_modulus = 206913;
100double poisson_ratio = 0.29;
101double rho = 0;
102double sigmaY = 450;
103double H = 129;
104double visH = 0;
105double cn0 = 1;
106double cn1 = 1;
107double zeta = 5e-2;
108double Qinf = 265;
109double b_iso = 16.93;
110
111int order = 2; ///< Order if fixed.
112int geom_order = 2; ///< Order if fixed.
113
114constexpr size_t activ_history_sise = 1;
115
116inline double hardening_exp(double tau) {
117 return std::exp(
118 std::max(static_cast<double>(std::numeric_limits<float>::min_exponent10),
119 -b_iso * tau));
120}
121
122inline double hardening(double tau) {
123 return H * tau + Qinf * (1. - hardening_exp(tau)) + sigmaY;
124}
125
126inline double hardening_dtau(double tau) {
127 auto r = [](auto tau) { return H + Qinf * b_iso * hardening_exp(tau); };
128 constexpr double eps = 1e-12;
129 return std::max(r(tau), eps * r(0));
130}
131
132inline double hardening_dtau2(double tau) {
133 return -(Qinf * (b_iso * b_iso)) * hardening_exp(tau);
134}
135
136#include <HenckyOps.hpp>
137#include <PlasticOps.hpp>
138
139using namespace PlasticOps;
140using namespace HenckyOps;
141struct Example {
142
143 Example(MoFEM::Interface &m_field) : mField(m_field) {}
144
146
147private:
149
155
156 boost::shared_ptr<PlasticOps::CommonData> commonPlasticDataPtr;
157 boost::shared_ptr<HenckyOps::CommonData> commonHenckyDataPtr;
158 boost::shared_ptr<DomainEle> reactionFe;
159
160 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
161 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
162 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
163
164 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
165 boost::shared_ptr<std::vector<unsigned char>> reactionMarker;
166
169 double getScale(const double time) {
170 return scale * MoFEM::TimeScale::getScale(time);
171 };
172 };
173};
174
175//! [Run problem]
180 CHKERR bC();
181 CHKERR OPs();
182 CHKERR tsSolve();
184}
185//! [Run problem]
186
187//! [Set up problem]
191
192 Range domain_ents;
193 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, domain_ents,
194 true);
195 auto get_ents_by_dim = [&](const auto dim) {
196 if (dim == SPACE_DIM) {
197 return domain_ents;
198 } else {
199 Range ents;
200 if (dim == 0)
201 CHKERR mField.get_moab().get_connectivity(domain_ents, ents, true);
202 else
203 CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents, true);
204 return ents;
205 }
206 };
207
208 auto get_base = [&]() {
209 auto domain_ents = get_ents_by_dim(SPACE_DIM);
210 if (domain_ents.empty())
211 CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Empty mesh");
212 const auto type = type_from_handle(domain_ents[0]);
213 switch (type) {
214 case MBQUAD:
216 case MBHEX:
218 case MBTRI:
220 case MBTET:
222 default:
223 CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Element type not handled");
224 }
225 return NOBASE;
226 };
227
228 const auto base = get_base();
229 MOFEM_LOG("WORLD", Sev::inform) << "Base " << ApproximationBaseNames[base];
230
231 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
232 CHKERR simple->addDomainField("EP", L2, base, size_symm);
233 CHKERR simple->addDomainField("TAU", L2, base, 1);
234 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
235
236 CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
237
238 auto ents = get_ents_by_dim(0);
239 ents.merge(get_ents_by_dim(1));
240 // ents.merge(get_ents_by_dim(2));
241 CHKERR simple->setFieldOrder("U", order, &ents);
242 CHKERR simple->setFieldOrder("EP", order - 1);
243 CHKERR simple->setFieldOrder("TAU", order - 2);
244
245 CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
246
247 CHKERR simple->setUp();
248 CHKERR simple->addFieldToEmptyFieldBlocks("U", "TAU");
249
250 auto project_ho_geometry = [&]() {
251 Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
252 return mField.loop_dofs("GEOMETRY", ent_method);
253 };
254 CHKERR project_ho_geometry();
255
257}
258//! [Set up problem]
259
260//! [Create common data]
263
264 auto get_command_line_parameters = [&]() {
266 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-scale", &scale, PETSC_NULL);
267 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho", &rho, PETSC_NULL);
268 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
269 &young_modulus, PETSC_NULL);
270 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
271 &poisson_ratio, PETSC_NULL);
272 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening", &H, PETSC_NULL);
273 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening_viscous", &visH,
274 PETSC_NULL);
275 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-yield_stress", &sigmaY,
276 PETSC_NULL);
277 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn0", &cn0, PETSC_NULL);
278 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn1", &cn1, PETSC_NULL);
279 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-zeta", &zeta, PETSC_NULL);
280 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-Qinf", &Qinf, PETSC_NULL);
281 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-b_iso", &b_iso, PETSC_NULL);
282 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-large_strains",
283 &is_large_strains, PETSC_NULL);
284 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-set_timer", &set_timer,
285 PETSC_NULL);
286
287 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
288 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-geom_order", &geom_order,
289 PETSC_NULL);
290
291 MOFEM_LOG("EXAMPLE", Sev::inform) << "Young modulus " << young_modulus;
292 MOFEM_LOG("EXAMPLE", Sev::inform) << "Poisson ratio " << poisson_ratio;
293 MOFEM_LOG("EXAMPLE", Sev::inform) << "Yield stress " << sigmaY;
294 MOFEM_LOG("EXAMPLE", Sev::inform) << "Hardening " << H;
295 MOFEM_LOG("EXAMPLE", Sev::inform) << "Viscous hardening " << visH;
296 MOFEM_LOG("EXAMPLE", Sev::inform) << "Saturation yield stress " << Qinf;
297 MOFEM_LOG("EXAMPLE", Sev::inform) << "Saturation exponent " << b_iso;
298 MOFEM_LOG("EXAMPLE", Sev::inform) << "cn0 " << cn0;
299 MOFEM_LOG("EXAMPLE", Sev::inform) << "cn1 " << cn1;
300 MOFEM_LOG("EXAMPLE", Sev::inform) << "zeta " << zeta;
301
302 MOFEM_LOG("EXAMPLE", Sev::inform) << "order " << order;
303 MOFEM_LOG("EXAMPLE", Sev::inform) << "geom order " << geom_order;
304
305 PetscBool is_scale = PETSC_TRUE;
306 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_scale", &is_scale,
307 PETSC_NULL);
308 if (is_scale) {
311 rho *= scale;
312 sigmaY *= scale;
313 H *= scale;
314 Qinf *= scale;
315 visH *= scale;
316
317 MOFEM_LOG("EXAMPLE", Sev::inform)
318 << "Scaled Young modulus " << young_modulus;
319 MOFEM_LOG("EXAMPLE", Sev::inform)
320 << "Scaled Poisson ratio " << poisson_ratio;
321 MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Yield stress " << sigmaY;
322 MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Hardening " << H;
323 MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Viscous hardening " << visH;
324 MOFEM_LOG("EXAMPLE", Sev::inform)
325 << "Scaled Saturation yield stress " << Qinf;
326 }
327
329 };
330
331 auto set_matrial_stiffness = [&]() {
338 const double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
339 const double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
340
341 // Plane stress or when 1, plane strain or 3d
342 const double A = (SPACE_DIM == 2)
343 ? 2 * shear_modulus_G /
344 (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
345 : 1;
346
347 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
348 *commonPlasticDataPtr->mDPtr);
349 auto t_D_axiator = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
350 *commonPlasticDataPtr->mDPtr_Axiator);
351 auto t_D_deviator = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
352 *commonPlasticDataPtr->mDPtr_Deviator);
353
354 constexpr double third = boost::math::constants::third<double>();
355 t_D_axiator(i, j, k, l) = A *
356 (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
357 t_kd(i, j) * t_kd(k, l);
358 t_D_deviator(i, j, k, l) =
359 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
360 t_D(i, j, k, l) = t_D_axiator(i, j, k, l) + t_D_deviator(i, j, k, l);
361
363 };
364
365 auto make_d_mat = []() {
366 return boost::make_shared<MatrixDouble>(size_symm * size_symm, 1);
367 };
368
369 commonPlasticDataPtr = boost::make_shared<PlasticOps::CommonData>();
370 commonPlasticDataPtr->mDPtr = make_d_mat();
371 commonPlasticDataPtr->mDPtr_Axiator = make_d_mat();
372 commonPlasticDataPtr->mDPtr_Deviator = make_d_mat();
373
374 commonPlasticDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
375 commonPlasticDataPtr->mStrainPtr = boost::make_shared<MatrixDouble>();
376 commonPlasticDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
377
378 CHKERR get_command_line_parameters();
379 CHKERR set_matrial_stiffness();
380
381 if (is_large_strains) {
382 commonHenckyDataPtr = boost::make_shared<HenckyOps::CommonData>();
383 commonHenckyDataPtr->matGradPtr = commonPlasticDataPtr->mGradPtr;
385 commonHenckyDataPtr->matLogCPlastic =
386 commonPlasticDataPtr->getPlasticStrainPtr();
387 commonPlasticDataPtr->mStrainPtr = commonHenckyDataPtr->getMatLogC();
388 commonPlasticDataPtr->mStressPtr =
389 commonHenckyDataPtr->getMatHenckyStress();
390 }
391
393}
394//! [Create common data]
395
396//! [Boundary condition]
399
401 auto bc_mng = mField.getInterface<BcManager>();
402 auto prb_mng = mField.getInterface<ProblemsManager>();
403
404 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
405 "U", 0, 0);
406 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
407 "U", 1, 1);
408 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
409 "U", 2, 2);
410 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
411 "REMOVE_ALL", "U", 0, 3);
412
413 CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
414 simple->getProblemName(), "U");
415
416 auto &bc_map = bc_mng->getBcMapByBlockName();
417 boundaryMarker = bc_mng->getMergedBlocksMarker(vector<string>{"FIX_"});
418
419 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "REACTION",
420 "U", 0, 3);
421
422 for (auto bc : bc_map)
423 MOFEM_LOG("EXAMPLE", Sev::verbose) << "Marker " << bc.first;
424
425 // OK. We have problem with GMesh, it adding empty characters at the end of
426 // block. So first block is search by regexp. popMarkDOFsOnEntities should
427 // work with regexp.
428 std::string reaction_block_set;
429 for (auto bc : bc_map) {
430 if (bc_mng->checkBlock(bc, "REACTION")) {
431 reaction_block_set = bc.first;
432 break;
433 }
434 }
435
436 if (auto bc = bc_mng->popMarkDOFsOnEntities(reaction_block_set)) {
437 reactionMarker = bc->getBcMarkersPtr();
438
439 // Only take reaction from nodes
440 Range nodes;
441 CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, nodes, true);
442 CHKERR prb_mng->markDofs(simple->getProblemName(), ROW,
443 ProblemsManager::MarkOP::AND, nodes,
445
446 } else {
447 MOFEM_LOG("EXAMPLE", Sev::warning) << "REACTION blockset does not exist";
448 }
449
450 if (!reactionMarker) {
451 MOFEM_LOG("EXAMPLE", Sev::warning) << "REACTION blockset does not exist";
452 }
453
455}
456//! [Boundary condition]
457
458//! [Push operators to pipeline]
461 auto pip = mField.getInterface<PipelineManager>();
463 auto bc_mng = mField.getInterface<BcManager>();
464
465 auto add_domain_base_ops = [&](auto &pipeline) {
467
469 "GEOMETRY");
470
471 pipeline.push_back(new OpCalculateScalarFieldValuesDot(
472 "TAU", commonPlasticDataPtr->getPlasticTauDotPtr()));
473 pipeline.push_back(new OpCalculateScalarFieldValues(
474 "TAU", commonPlasticDataPtr->getPlasticTauPtr()));
476 "EP", commonPlasticDataPtr->getPlasticStrainPtr()));
478 "EP", commonPlasticDataPtr->getPlasticStrainDotPtr()));
480 "U", commonPlasticDataPtr->mGradPtr));
481
483 };
484
485 auto add_domain_stress_ops = [&](auto &pipeline, auto m_D_ptr) {
487
488 if (is_large_strains) {
489
490 if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
491 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
492 "Wrong pointer for grad");
493
494 pipeline.push_back(
496 pipeline.push_back(
498 pipeline.push_back(
500 pipeline.push_back(new OpCalculateHenckyPlasticStress<SPACE_DIM>(
501 "U", commonHenckyDataPtr, m_D_ptr));
502 pipeline.push_back(
504
505 } else {
506 pipeline.push_back(
508 commonPlasticDataPtr->mStrainPtr));
509 pipeline.push_back(
510 new OpPlasticStress("U", commonPlasticDataPtr, m_D_ptr, 1));
511 }
512
513 if (m_D_ptr != commonPlasticDataPtr->mDPtr_Axiator) {
514 pipeline.push_back(
516 pipeline.push_back(new OpCalculatePlasticity(
517 "U", mField, commonPlasticDataPtr, m_D_ptr));
518 }
519
521 };
522
523 auto add_domain_ops_lhs_mechanical = [&](auto &pipeline, auto m_D_ptr) {
525 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
526
527 if (is_large_strains) {
528 pipeline.push_back(
530 pipeline.push_back(
531 new OpKPiola("U", "U", commonHenckyDataPtr->getMatTangent()));
533 "U", "EP", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
534 } else {
535 pipeline.push_back(new OpKCauchy("U", "U", m_D_ptr));
536 pipeline.push_back(new OpCalculatePlasticInternalForceLhs_dEP(
537 "U", "EP", commonPlasticDataPtr, m_D_ptr));
538 }
539
540 pipeline.push_back(new OpUnSetBc("U"));
542 };
543
544 auto add_domain_ops_rhs_mechanical = [&](auto &pipeline) {
546 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
547
549 pipeline, mField, "U", {boost::make_shared<PlasticityTimeScale>()},
550 "BODY_FORCE", Sev::inform);
551
552 // Calculate internal forces
553 if (is_large_strains) {
554 pipeline.push_back(new OpInternalForcePiola(
555 "U", commonHenckyDataPtr->getMatFirstPiolaStress()));
556 } else {
557 pipeline.push_back(
558 new OpInternalForceCauchy("U", commonPlasticDataPtr->mStressPtr));
559 }
560
561 pipeline.push_back(new OpUnSetBc("U"));
563 };
564
565 auto add_domain_ops_lhs_constrain = [&](auto &pipeline, auto m_D_ptr) {
567 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
568
569 if (is_large_strains) {
570 pipeline.push_back(new OpCalculateConstraintsLhs_LogStrain_dU(
571 "TAU", "U", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
572 pipeline.push_back(new OpCalculatePlasticFlowLhs_LogStrain_dU(
573 "EP", "U", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
574 } else {
575 pipeline.push_back(new OpCalculateConstraintsLhs_dU(
576 "TAU", "U", commonPlasticDataPtr, m_D_ptr));
577 pipeline.push_back(new OpCalculatePlasticFlowLhs_dU(
578 "EP", "U", commonPlasticDataPtr, m_D_ptr));
579 }
580
581 pipeline.push_back(new OpCalculatePlasticFlowLhs_dEP(
582 "EP", "EP", commonPlasticDataPtr, m_D_ptr));
583 pipeline.push_back(new OpCalculatePlasticFlowLhs_dTAU(
584 "EP", "TAU", commonPlasticDataPtr, m_D_ptr));
585 pipeline.push_back(new OpCalculateConstraintsLhs_dEP(
586 "TAU", "EP", commonPlasticDataPtr, m_D_ptr));
587 pipeline.push_back(
589
590 pipeline.push_back(new OpUnSetBc("U"));
592 };
593
594 auto add_domain_ops_rhs_constrain = [&](auto &pipeline, auto m_D_ptr) {
596 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
597
598 pipeline.push_back(
600 pipeline.push_back(
602
603 pipeline.push_back(new OpUnSetBc("U"));
605 };
606
607 auto add_boundary_ops_lhs_mechanical = [&](auto &pipeline) {
611 mField, pipeline, simple->getProblemName(), "U");
613 };
614
615 auto add_boundary_ops_rhs_mechanical = [&](auto &pipeline) {
617
619 pipeline, {NOSPACE}, "GEOMETRY");
620
621 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
623 pipeline, mField, "U", {boost::make_shared<PlasticityTimeScale>()},
624 "FORCE", Sev::inform);
625
626 pipeline.push_back(new OpUnSetBc("U"));
627
628 auto u_mat_ptr = boost::make_shared<MatrixDouble>();
629 pipeline.push_back(
630 new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_mat_ptr));
631
634 mField, pipeline, simple->getProblemName(), "U", u_mat_ptr,
635 {boost::make_shared<TimeScale>()}); // note displacements have no
636 // scaling
637
639 };
640
641 CHKERR add_domain_base_ops(pip->getOpDomainLhsPipeline());
642 CHKERR add_domain_stress_ops(pip->getOpDomainLhsPipeline(),
643 commonPlasticDataPtr->mDPtr);
644 CHKERR add_domain_ops_lhs_mechanical(pip->getOpDomainLhsPipeline(),
645 commonPlasticDataPtr->mDPtr);
646 CHKERR add_domain_ops_lhs_constrain(pip->getOpDomainLhsPipeline(),
647 commonPlasticDataPtr->mDPtr);
648 CHKERR
649 add_boundary_ops_lhs_mechanical(pip->getOpBoundaryLhsPipeline());
650
651 CHKERR add_domain_base_ops(pip->getOpDomainRhsPipeline());
652 CHKERR add_domain_stress_ops(pip->getOpDomainRhsPipeline(),
653 commonPlasticDataPtr->mDPtr);
654 CHKERR add_domain_ops_rhs_mechanical(pip->getOpDomainRhsPipeline());
655 CHKERR add_domain_ops_rhs_constrain(pip->getOpDomainRhsPipeline(),
656 commonPlasticDataPtr->mDPtr);
657
658 // Boundary
659 CHKERR
660 add_boundary_ops_rhs_mechanical(pip->getOpBoundaryRhsPipeline());
661
662 auto integration_rule_bc = [](int, int, int ao) { return 2 * ao; };
663
664 auto vol_rule = [](int, int, int ao) { return 2 * ao + geom_order - 1; };
665
666 CHKERR pip->setDomainRhsIntegrationRule(vol_rule);
667 CHKERR pip->setDomainLhsIntegrationRule(vol_rule);
668
669 CHKERR pip->setBoundaryLhsIntegrationRule(integration_rule_bc);
670 CHKERR pip->setBoundaryRhsIntegrationRule(integration_rule_bc);
671
672 auto create_reaction_pipeline = [&](auto &pipeline) {
674
675 if (reactionMarker) {
676
678 "GEOMETRY");
679
680 pipeline.push_back(
682 "U", commonPlasticDataPtr->mGradPtr));
684 "EP", commonPlasticDataPtr->getPlasticStrainPtr()));
685
686 if (is_large_strains) {
687
688 if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
689 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
690 "Wrong pointer for grad");
691
692 pipeline.push_back(
694 pipeline.push_back(
696 pipeline.push_back(
698 pipeline.push_back(new OpCalculateHenckyPlasticStress<SPACE_DIM>(
700 pipeline.push_back(
702
703 } else {
704 pipeline.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
705 "U", commonPlasticDataPtr->mGradPtr,
706 commonPlasticDataPtr->mStrainPtr));
707 pipeline.push_back(new OpPlasticStress(
709 }
710
711 pipeline.push_back(new OpSetBc("U", false, reactionMarker));
712 // Calculate internal force
713 if (is_large_strains) {
714 pipeline.push_back(new OpInternalForcePiola(
715 "U", commonHenckyDataPtr->getMatFirstPiolaStress()));
716 } else {
717 pipeline.push_back(
718 new OpInternalForceCauchy("U", commonPlasticDataPtr->mStressPtr));
719 }
720 pipeline.push_back(new OpUnSetBc("U"));
721 }
722
724 };
725
726 reactionFe = boost::make_shared<DomainEle>(mField);
727 reactionFe->getRuleHook = vol_rule;
728
729 CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
730
732}
733//! [Push operators to pipeline]
734
735//! [Solve]
737 static boost::shared_ptr<SetUpSchur> createSetUpSchur(
738
741
742 );
743 virtual MoFEMErrorCode setUp(KSP solver) = 0;
744
745 virtual MoFEMErrorCode preProc() = 0;
747
748protected:
749 SetUpSchur() = default;
750};
751
754
757 ISManager *is_manager = mField.getInterface<ISManager>();
758
759 auto snes_ctx_ptr = smartGetDMSnesCtx(simple->getDM());
760
761 auto set_section_monitor = [&](auto solver) {
763 SNES snes;
764 CHKERR TSGetSNES(solver, &snes);
765 CHKERR SNESMonitorSet(snes,
766 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
768 (void *)(snes_ctx_ptr.get()), nullptr);
770 };
771
772 auto create_post_process_element = [&]() {
773 auto pp_fe = boost::make_shared<PostProcEle>(mField);
775
777 pp_fe->getOpPtrVector(), {H1}, "GEOMETRY");
778
779 auto x_ptr = boost::make_shared<MatrixDouble>();
780 pp_fe->getOpPtrVector().push_back(
781 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
782 auto u_ptr = boost::make_shared<MatrixDouble>();
783 pp_fe->getOpPtrVector().push_back(
785
786 pp_fe->getOpPtrVector().push_back(
788 "U", commonPlasticDataPtr->mGradPtr));
789 pp_fe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
790 "TAU", commonPlasticDataPtr->getPlasticTauPtr()));
791 pp_fe->getOpPtrVector().push_back(
793 "EP", commonPlasticDataPtr->getPlasticStrainPtr()));
794
795 if (is_large_strains) {
796
797 if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
798 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong pointer for grad");
799
800 pp_fe->getOpPtrVector().push_back(
802 pp_fe->getOpPtrVector().push_back(
804 pp_fe->getOpPtrVector().push_back(
806 pp_fe->getOpPtrVector().push_back(
809 pp_fe->getOpPtrVector().push_back(
811
812 pp_fe->getOpPtrVector().push_back(
813
814 new OpPPMap(
815
816 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
817
818 {},
819
820 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
821
822 {{"GRAD", commonPlasticDataPtr->mGradPtr},
823 {"FIRST_PIOLA", commonHenckyDataPtr->getMatFirstPiolaStress()}},
824
825 {}
826
827 )
828
829 );
830
831 } else {
832 pp_fe->getOpPtrVector().push_back(
834 commonPlasticDataPtr->mStrainPtr));
835 pp_fe->getOpPtrVector().push_back(new OpPlasticStress(
837
838 pp_fe->getOpPtrVector().push_back(
839
840 new OpPPMap(
841
842 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
843
844 {},
845
846 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
847
848 {},
849
850 {{"STRAIN", commonPlasticDataPtr->mStrainPtr},
851 {"STRESS", commonPlasticDataPtr->mStressPtr}}
852
853 )
854
855 );
856 }
857
858 pp_fe->getOpPtrVector().push_back(
859 new OpCalculatePlasticSurface("U", commonPlasticDataPtr));
860
861 pp_fe->getOpPtrVector().push_back(
862
863 new OpPPMap(
864
865 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
866
867 {{"PLASTIC_SURFACE", commonPlasticDataPtr->getPlasticSurfacePtr()},
868 {"PLASTIC_MULTIPLIER", commonPlasticDataPtr->getPlasticTauPtr()}},
869
870 {},
871
872 {},
873
874 {{"PLASTIC_STRAIN", commonPlasticDataPtr->getPlasticStrainPtr()},
875 {"PLASTIC_FLOW", commonPlasticDataPtr->getPlasticFlowPtr()}}
876
877 )
878
879 );
880
881 return pp_fe;
882 };
883
884 auto scatter_create = [&](auto D, auto coeff) {
886 CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
887 ROW, "U", coeff, coeff, is);
888 int loc_size;
889 CHKERR ISGetLocalSize(is, &loc_size);
890 Vec v;
891 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
892 VecScatter scatter;
893 CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
894 return std::make_tuple(SmartPetscObj<Vec>(v),
896 };
897
898 auto set_time_monitor = [&](auto dm, auto solver) {
900 boost::shared_ptr<Monitor> monitor_ptr(
901 new Monitor(dm, create_post_process_element(), reactionFe, uXScatter,
902 uYScatter, uZScatter));
903 boost::shared_ptr<ForcesAndSourcesCore> null;
904 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
905 monitor_ptr, null, null);
907 };
908
909 auto set_fieldsplit_preconditioner = [&](auto solver,
910 boost::shared_ptr<SetUpSchur>
911 &schur_ptr) {
913
914 SNES snes;
915 CHKERR TSGetSNES(solver, &snes);
916 KSP ksp;
917 CHKERR SNESGetKSP(snes, &ksp);
918 PC pc;
919 CHKERR KSPGetPC(ksp, &pc);
920 PetscBool is_pcfs = PETSC_FALSE;
921 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
922
923 // Setup fieldsplit (block) solver - optional: yes/no
924 if (is_pcfs == PETSC_TRUE) {
925
926 auto bc_mng = mField.getInterface<BcManager>();
927 auto name_prb = simple->getProblemName();
928
929 auto create_sub_bc_dm = [&](SmartPetscObj<DM> base_dm,
930 SmartPetscObj<DM> &dm_sub,
931 SmartPetscObj<IS> &is_sub,
932 SmartPetscObj<AO> &ao_sub) {
934
935 dm_sub = createSmartDM(mField.get_comm(), "DMMOFEM");
936 CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "SUB_BC");
937 CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
938 CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
939 CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
940 for (auto f : {"U", "EP", "TAU"}) {
943 }
944 CHKERR DMSetUp(dm_sub);
945
946 CHKERR bc_mng->removeBlockDOFsOnEntities("SUB_BC", "FIX_X", "U", 0, 0);
947 CHKERR bc_mng->removeBlockDOFsOnEntities("SUB_BC", "FIX_Y", "U", 1, 1);
948 CHKERR bc_mng->removeBlockDOFsOnEntities("SUB_BC", "FIX_Z", "U", 2, 2);
949 CHKERR bc_mng->removeBlockDOFsOnEntities("SUB_BC", "FIX_ALL", "U", 0,
950 2);
951
952 auto *prb_ptr = getProblemPtr(dm_sub);
953 if (auto sub_data = prb_ptr->getSubData()) {
954 is_sub = sub_data->getSmartRowIs();
955 ao_sub = sub_data->getSmartRowMap();
956 int is_sub_size;
957 CHKERR ISGetSize(is_sub, &is_sub_size);
958 MOFEM_LOG("EXAMPLE", Sev::inform)
959 << "Field split second block size " << is_sub_size;
960
961 } else {
962 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "No sub data");
963 }
964
966 };
967
968 auto create_sub_u_dm = [&](SmartPetscObj<DM> base_dm,
969 SmartPetscObj<DM> &dm_sub,
970 SmartPetscObj<IS> &is_sub,
971 SmartPetscObj<AO> &ao_sub) {
973 dm_sub = createSmartDM(mField.get_comm(), "DMMOFEM");
974 CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "SUB_U");
975 CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
976 CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
977 CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
978 for (auto f : {"U"}) {
981 }
982 CHKERR DMSetUp(dm_sub);
983
984 auto *prb_ptr = getProblemPtr(dm_sub);
985 if (auto sub_data = prb_ptr->getSubData()) {
986 is_sub = sub_data->getSmartRowIs();
987 ao_sub = sub_data->getSmartRowMap();
988 // ISView(is_sub, PETSC_VIEWER_STDOUT_WORLD);
989 int is_sub_size;
990 CHKERR ISGetSize(is_sub, &is_sub_size);
991 MOFEM_LOG("EXAMPLE", Sev::inform)
992 << "Field split second block size " << is_sub_size;
993
994 } else {
995 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "No sub data");
996 }
997
998 CHKERR DMSetUp(dm_sub);
1000 };
1001
1002 auto create_all_bc_is = [&](SmartPetscObj<IS> &is_all_bc) {
1004 is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_X", "U", 0, 0);
1005 is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_Y", "U", 1, 1, is_all_bc);
1006 is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_Z", "U", 2, 2, is_all_bc);
1007 is_all_bc =
1008 bc_mng->getBlockIS(name_prb, "FIX_ALL", "U", 0, 2, is_all_bc);
1009 int is_all_bc_size;
1010 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
1011 MOFEM_LOG("EXAMPLE", Sev::inform)
1012 << "Field split first block size " << is_all_bc_size;
1014 };
1015
1016 SmartPetscObj<IS> is_all_bc;
1017 SmartPetscObj<DM> dm_bc_sub;
1018 SmartPetscObj<IS> is_bc_sub;
1019 SmartPetscObj<AO> ao_bc_sub;
1020
1021 CHKERR create_all_bc_is(is_all_bc);
1022 CHKERR create_sub_bc_dm(simple->getDM(), dm_bc_sub, is_bc_sub, ao_bc_sub);
1023
1024 SmartPetscObj<IS> is_prb;
1025 CHKERR mField.getInterface<ISManager>()->isCreateProblem(
1026 simple->getProblemName(), ROW, is_prb);
1027
1028 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
1029 is_all_bc); // boundary block
1030 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_bc_sub);
1031
1032 if constexpr (A == AssemblyType::SCHUR) {
1033
1034 SmartPetscObj<IS> is_epp;
1035 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1036 "SUB_BC", ROW, "EP", 0, MAX_DOFS_ON_ENTITY, is_epp);
1037 SmartPetscObj<IS> is_tau;
1038 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1039 "SUB_BC", ROW, "TAU", 0, MAX_DOFS_ON_ENTITY, is_tau);
1040 IS is_union_raw;
1041 CHKERR ISExpand(is_epp, is_tau, &is_union_raw);
1042 SmartPetscObj<IS> is_union(is_union_raw);
1043
1044 SmartPetscObj<DM> dm_u_sub;
1045 SmartPetscObj<IS> is_u_sub;
1046 SmartPetscObj<AO> ao_u_sub;
1047 CHKERR create_sub_u_dm(dm_bc_sub, dm_u_sub, is_u_sub, ao_u_sub);
1048
1049 // Indices has to be map fro very to level, while assembling Schur
1050 // complement.
1051 auto is_up = is_u_sub;
1052 CHKERR AOPetscToApplicationIS(ao_bc_sub, is_up);
1053 auto ao_up = createAOMappingIS(is_u_sub, PETSC_NULL);
1054 schur_ptr =
1055 SetUpSchur::createSetUpSchur(mField, dm_u_sub, is_union, ao_up);
1056 PetscInt n;
1057 KSP *ksps;
1058 CHKERR PCFieldSplitGetSubKSP(pc, &n, &ksps);
1059 CHKERR schur_ptr->setUp(ksps[1]);
1060 }
1061 }
1062
1064 };
1065
1066 auto dm = simple->getDM();
1067 auto D = smartCreateDMVector(dm);
1068 uXScatter = scatter_create(D, 0);
1069 uYScatter = scatter_create(D, 1);
1070 if constexpr (SPACE_DIM == 3)
1071 uZScatter = scatter_create(D, 2);
1072
1073 auto solver = pip->createTSIM();
1074
1075 auto active_pre_lhs = [&]() {
1077 std::fill(commonPlasticDataPtr->activityData.begin(),
1078 commonPlasticDataPtr->activityData.end(), 0);
1080 };
1081
1082 auto active_post_lhs = [&]() {
1084 auto get_iter = [&]() {
1085 SNES snes;
1086 CHK_THROW_MESSAGE(TSGetSNES(solver, &snes), "Can not get SNES");
1087 int iter;
1088 CHK_THROW_MESSAGE(SNESGetIterationNumber(snes, &iter),
1089 "Can not get iter");
1090 return iter;
1091 };
1092
1093 auto iter = get_iter();
1094 if (iter >= 0) {
1095
1096 std::array<int, 5> activity_data;
1097 std::fill(activity_data.begin(), activity_data.end(), 0);
1098 MPI_Allreduce(commonPlasticDataPtr->activityData.data(),
1099 activity_data.data(), activity_data.size(), MPI_INT,
1100 MPI_SUM, mField.get_comm());
1101
1102 int &active_points = activity_data[0];
1103 int &avtive_full_elems = activity_data[1];
1104 int &avtive_elems = activity_data[2];
1105 int &nb_points = activity_data[3];
1106 int &nb_elements = activity_data[4];
1107
1108 if (nb_points) {
1109
1110 double proc_nb_points =
1111 100 * static_cast<double>(active_points) / nb_points;
1112 double proc_nb_active =
1113 100 * static_cast<double>(avtive_elems) / nb_elements;
1114 double proc_nb_full_active = 100;
1115 if (avtive_elems)
1116 proc_nb_full_active =
1117 100 * static_cast<double>(avtive_full_elems) / avtive_elems;
1118
1120 "EXAMPLE", Sev::inform,
1121 "Iter %d nb pts %d nb avtive pts %d (%3.3f\%) nb active elements %d "
1122 "(%3.3f\%) nb full active elems %d (%3.3f\%)",
1123 iter, nb_points, active_points, proc_nb_points, avtive_elems,
1124 proc_nb_active, avtive_full_elems, proc_nb_full_active, iter);
1125 }
1126 }
1127
1129 };
1130
1131 CHKERR TSSetSolution(solver, D);
1132 CHKERR set_section_monitor(solver);
1133 CHKERR set_time_monitor(dm, solver);
1134 CHKERR TSSetSolution(solver, D);
1135 CHKERR TSSetFromOptions(solver);
1136
1137 boost::shared_ptr<SetUpSchur> schur_ptr;
1138 CHKERR set_fieldsplit_preconditioner(solver, schur_ptr);
1139
1140 mField.getInterface<PipelineManager>()->getDomainLhsFE()->preProcessHook =
1141 [&]() {
1143 if (schur_ptr)
1144 CHKERR schur_ptr->preProc();
1145 CHKERR active_pre_lhs();
1147 };
1148 mField.getInterface<PipelineManager>()->getDomainLhsFE()->postProcessHook =
1149 [&]() {
1151 // if (schur_ptr)
1152 // CHKERR schur_ptr->postProc();
1153 CHKERR active_post_lhs();
1155 };
1156
1157 mField.getInterface<PipelineManager>()->getBoundaryLhsFE()->postProcessHook =
1158 [&]() {
1160 if (schur_ptr)
1161 CHKERR schur_ptr->postProc();
1163 };
1164
1165 MOFEM_LOG_CHANNEL("TIMER");
1166 MOFEM_LOG_TAG("TIMER", "timer");
1167 if(set_timer)
1168 BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
1169 MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp";
1170 CHKERR TSSetUp(solver);
1171 MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp <= done";
1172 MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve";
1173 CHKERR TSSolve(solver, NULL);
1174 MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve <= done";
1175
1176 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1177 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1178 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
1179
1181}
1182//! [Solve]
1183
1184static char help[] = "...\n\n";
1185
1186int main(int argc, char *argv[]) {
1187
1188 // Initialisation of MoFEM/PETSc and MOAB data structures
1189 const char param_file[] = "param_file.petsc";
1190 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
1191
1192 // Add logging channel for example
1193 auto core_log = logging::core::get();
1194 core_log->add_sink(
1196 core_log->add_sink(
1198 LogManager::setLog("EXAMPLE");
1199
1200 MOFEM_LOG_TAG("EXAMPLE", "example");
1201
1202 try {
1203
1204 //! [Register MoFEM discrete manager in PETSc]
1205 DMType dm_name = "DMMOFEM";
1206 CHKERR DMRegister_MoFEM(dm_name);
1207 //! [Register MoFEM discrete manager in PETSc
1208
1209 //! [Create MoAB]
1210 moab::Core mb_instance; ///< mesh database
1211 moab::Interface &moab = mb_instance; ///< mesh database interface
1212 //! [Create MoAB]
1213
1214 //! [Create MoFEM]
1215 MoFEM::Core core(moab); ///< finite element database
1216 MoFEM::Interface &m_field = core; ///< finite element database interface
1217 //! [Create MoFEM]
1218
1219 //! [Load mesh]
1220 Simple *simple = m_field.getInterface<Simple>();
1221 CHKERR simple->getOptions();
1222 CHKERR simple->loadFile();
1223 //! [Load mesh]
1224
1225 //! [Example]
1226 Example ex(m_field);
1227 CHKERR ex.runProblem();
1228 //! [Example]
1229 }
1231
1233}
1234
1236
1238 SmartPetscObj<IS> sub_is, SmartPetscObj<AO> sub_ao)
1239 : SetUpSchur(), mField(m_field), subDM(sub_dm), subIS(sub_is),
1240 subAO(sub_ao) {
1241 if (S) {
1244 "Is expected that schur matrix is not allocated. This is "
1245 "possible only is PC is set up twice");
1246 }
1247 }
1248 virtual ~SetUpSchurImpl() { S.reset(); }
1249
1250 MoFEMErrorCode setUp(KSP solver);
1253
1254private:
1256 MoFEMErrorCode setPC(PC pc);
1257
1260
1265};
1266
1269 auto pip = mField.getInterface<PipelineManager>();
1270 PC pc;
1271 CHKERR KSPSetFromOptions(solver);
1272 CHKERR KSPGetPC(solver, &pc);
1273 PetscBool is_pcfs = PETSC_FALSE;
1274 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1275 if (is_pcfs) {
1276 if (S) {
1279 "Is expected that schur matrix is not allocated. This is "
1280 "possible only is PC is set up twice");
1281 }
1284 CHKERR setPC(pc);
1285
1286 } else {
1287 pip->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
1288 pip->getOpBoundaryLhsPipeline().push_back(
1289 new OpSchurAssembleEnd({}, {}, {}, {}, {}));
1290 pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
1291 pip->getOpDomainLhsPipeline().push_back(
1292 new OpSchurAssembleEnd({}, {}, {}, {}, {}));
1293 }
1294
1295 subDM.reset();
1296 subIS.reset();
1297 subAO.reset();
1299}
1300
1303 auto pip = mField.getInterface<PipelineManager>();
1304 // Boundary
1305 pip->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
1306 pip->getOpBoundaryLhsPipeline().push_back(new OpSchurAssembleEnd(
1307
1308 {"EP", "TAU"}, {nullptr, nullptr}, {SmartPetscObj<AO>(), subAO},
1309 {SmartPetscObj<Mat>(), S}, {false, false}
1310
1311 ));
1312 // Domain
1313 pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
1314 pip->getOpDomainLhsPipeline().push_back(new OpSchurAssembleEnd(
1315
1316 {"EP", "TAU"}, {nullptr, nullptr}, {SmartPetscObj<AO>(), subAO},
1317 {SmartPetscObj<Mat>(), S}, {false, false}
1318
1319 ));
1321}
1322
1325 smartPC = SmartPetscObj<PC>(pc, true);
1326 CHKERR PCFieldSplitSetIS(pc, NULL, subIS);
1327 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1329}
1330
1333 if (SetUpSchurImpl::S) {
1334 CHKERR MatZeroEntries(S);
1335 }
1336 MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble Begin";
1338}
1339
1342 MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble End";
1343 if (S) {
1344 CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
1345 CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
1346 }
1348}
1349
1350boost::shared_ptr<SetUpSchur>
1352 SmartPetscObj<DM> sub_dm, SmartPetscObj<IS> is_sub,
1353 SmartPetscObj<AO> ao_sub) {
1354 return boost::shared_ptr<SetUpSchur>(
1355 new SetUpSchurImpl(m_field, sub_dm, is_sub, ao_sub));
1356}
std::string param_file
constexpr double third
#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
static const double eps
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
@ 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
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ NOSPACE
Definition: definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ 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 double shear_modulus_G
Definition: elastic.cpp:60
constexpr double bulk_modulus_K
Definition: elastic.cpp:59
const int dim
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FTensor::Index< 'n', SPACE_DIM > n
constexpr auto t_kd
auto smartCreateDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:983
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:221
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:485
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:444
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:511
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:244
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:277
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:995
boost::ptr_deque< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
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:391
#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.
double D
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'k', 3 > k
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
const FTensor::Tensor2< T, Dim, Dim > Vec
auto f
Definition: HenckyOps.hpp:5
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:1634
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:1044
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:224
auto createSmartDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
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)
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
auto getProblemPtr(DM dm)
get problem pointer from DM
Definition: DMMoFEM.hpp:971
auto smartGetDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:1017
constexpr AssemblyType A
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double young_modulus
Definition: plastic.cpp:99
double Qinf
Definition: plastic.cpp:108
static char help[]
[Solve]
Definition: plastic.cpp:1184
double hardening_dtau2(double tau)
Definition: plastic.cpp:132
constexpr size_t activ_history_sise
Definition: plastic.cpp:114
double rho
Definition: plastic.cpp:101
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:10
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< G >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: plastic.cpp:68
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Only used with Hooke equation (linear material model)]
Definition: plastic.cpp:52
constexpr int SPACE_DIM
Definition: plastic.cpp:31
double visH
Definition: plastic.cpp:104
double poisson_ratio
Definition: plastic.cpp:100
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< G >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:73
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< G >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: plastic.cpp:54
PetscBool set_timer
Definition: plastic.cpp:95
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< G >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: plastic.cpp:75
double scale
Definition: plastic.cpp:97
constexpr auto size_symm
Definition: plastic.cpp:33
double zeta
Definition: plastic.cpp:107
double H
Definition: plastic.cpp:103
double hardening(double tau)
Definition: plastic.cpp:122
constexpr AssemblyType A
Definition: plastic.cpp:35
double cn0
Definition: plastic.cpp:105
int order
Order if fixed.
Definition: plastic.cpp:111
double b_iso
Definition: plastic.cpp:109
PetscBool is_large_strains
Definition: plastic.cpp:94
int geom_order
Order if fixed.
Definition: plastic.cpp:112
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: plastic.cpp:66
double sigmaY
Definition: plastic.cpp:102
double hardening_dtau(double tau)
Definition: plastic.cpp:126
constexpr IntegrationType G
Definition: plastic.cpp:38
double hardening_exp(double tau)
Definition: plastic.cpp:116
double cn1
Definition: plastic.cpp:106
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< G >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
Definition: plastic.cpp:77
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
double getScale(const double time)
Get scaling at a given time.
Definition: plastic.cpp:169
[Example]
Definition: plastic.cpp:141
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
Definition: plastic.cpp:752
FieldApproximationBase base
Definition: plot_base.cpp:68
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: plastic.cpp:161
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:261
Example(MoFEM::Interface &m_field)
Definition: plastic.cpp:143
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:459
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:176
MoFEM::Interface & mField
Definition: plastic.cpp:148
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:162
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: plastic.cpp:164
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:188
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:397
boost::shared_ptr< std::vector< unsigned char > > reactionMarker
Definition: plastic.cpp:165
boost::shared_ptr< DomainEle > reactionFe
Definition: plastic.cpp:158
boost::shared_ptr< PlasticOps::CommonData > commonPlasticDataPtr
Definition: plastic.cpp:156
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:160
boost::shared_ptr< HenckyOps::CommonData > commonHenckyDataPtr
Definition: plastic.cpp:157
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:24
virtual moab::Interface & get_moab()=0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Definition of the displacement bc data structure.
Definition: BCData.hpp:72
Data on single entity (This is passed as argument to DataOperator::doWork)
Essential boundary conditions.
Definition: Essential.hpp:101
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:300
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:346
Assembly methods.
Definition: Natural.hpp:57
Get value at integration points for scalar field.
Calculate symmetric tensor field rates ant integratio pts.
Calculate symmetric tensor field values at integration pts.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Enforce essential constrains on lhs.
Definition: Essential.hpp:71
Enforce essential constrains on rhs.
Definition: Essential.hpp:55
Post post-proc data at points from hash maps.
Scale base functions by inverses of measure of element.
Clear Schur complement internal data.
Definition: Schur.hpp:22
Assemble Schur complement.
Definition: Schur.hpp:77
Set indices on entities on finite element.
PipelineManager interface.
MoFEM::VolumeElementForcesAndSourcesCore VolEle
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]
Definition: plastic.cpp:736
SetUpSchur()=default
virtual MoFEMErrorCode setUp(KSP solver)=0
virtual MoFEMErrorCode postProc()=0
virtual MoFEMErrorCode preProc()=0
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, SmartPetscObj< DM >, SmartPetscObj< IS >, SmartPetscObj< AO >)
Definition: plastic.cpp:1351
SmartPetscObj< PC > smartPC
Definition: plastic.cpp:1259
SmartPetscObj< DM > subDM
Definition: plastic.cpp:1262
MoFEMErrorCode setUp(KSP solver)
Definition: plastic.cpp:1267
SmartPetscObj< Mat > S
Definition: plastic.cpp:1258
SetUpSchurImpl(MoFEM::Interface &m_field, SmartPetscObj< DM > sub_dm, SmartPetscObj< IS > sub_is, SmartPetscObj< AO > sub_ao)
Definition: plastic.cpp:1237
SmartPetscObj< AO > subAO
Definition: plastic.cpp:1264
virtual ~SetUpSchurImpl()
Definition: plastic.cpp:1248
MoFEMErrorCode postProc()
Definition: plastic.cpp:1340
MoFEMErrorCode setPC(PC pc)
Definition: plastic.cpp:1323
SmartPetscObj< IS > subIS
Definition: plastic.cpp:1263
MoFEMErrorCode setOperator()
Definition: plastic.cpp:1301
MoFEMErrorCode preProc()
Definition: plastic.cpp:1331
MoFEM::Interface & mField
Definition: plastic.cpp:1261