v0.13.1
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/* This file is part of MoFEM.
10 * MoFEM is free software: you can redistribute it and/or modify it under
11 * the terms of the GNU Lesser General Public License as published by the
12 * Free Software Foundation, either version 3 of the License, or (at your
13 * option) any later version.
14 *
15 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18 * License for more details.
19 *
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
22
23#ifndef EXECUTABLE_DIMENSION
24#define EXECUTABLE_DIMENSION 3
25#endif
26
27#include <MoFEM.hpp>
28#include <MatrixFunction.hpp>
29#include <IntegrationRules.hpp>
30
31using namespace MoFEM;
32
33template <int DIM> struct ElementsAndOps {};
34
35template <> struct ElementsAndOps<2> {
41};
42
43template <> struct ElementsAndOps<3> {
49};
50
51constexpr int SPACE_DIM =
52 EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
53
54constexpr EntityType boundary_ent = SPACE_DIM == 3 ? MBTRI : MBEDGE;
61
64
65//! [Body force]
67 GAUSS>::OpSource<1, SPACE_DIM>;
68//! [Body force]
69
70//! [Only used with Hooke equation (linear material model)]
72 GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
75//! [Only used with Hooke equation (linear material model)]
76
77//! [Only used for dynamics]
82//! [Only used for dynamics]
83
84//! [Only used with Hencky/nonlinear material]
86 GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
89//! [Only used with Hencky/nonlinear material]
90
91//! [Essential boundary conditions]
98//! [Essential boundary conditions]
100
101PetscBool is_large_strains = PETSC_TRUE;
102
103double scale = 1.;
104
105double young_modulus = 206913;
106double poisson_ratio = 0.29;
107double rho = 0;
108double sigmaY = 450;
109double H = 129;
110double visH = 0;
111double cn = 1;
112double Qinf = 265;
113double b_iso = 16.93;
114int order = 2;
115
116inline long double hardening(long double tau, double temp) {
117 return H * tau + Qinf * (1. - std::exp(-b_iso * tau)) + sigmaY;
118}
119
120inline long double hardening_dtau(long double tau, double temp) {
121 return H + Qinf * b_iso * std::exp(-b_iso * tau);
122}
123
124#include <HenckyOps.hpp>
125#include <PlasticOps.hpp>
126#include <OpPostProcElastic.hpp>
127
128using namespace PlasticOps;
129using namespace HenckyOps;
130struct Example {
131
132 Example(MoFEM::Interface &m_field) : mField(m_field) {}
133
135
136private:
138
144
145 boost::shared_ptr<PlasticOps::CommonData> commonPlasticDataPtr;
146 boost::shared_ptr<HenckyOps::CommonData> commonHenckyDataPtr;
147 boost::shared_ptr<PostProcEle> postProcFe;
148 boost::shared_ptr<DomainEle> reactionFe;
149 boost::shared_ptr<DomainEle> feAxiatorRhs;
150 boost::shared_ptr<DomainEle> feAxiatorLhs;
151
152 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
153 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
154 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
155
156 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
157 boost::shared_ptr<std::vector<unsigned char>> reactionMarker;
158
159 std::vector<FTensor::Tensor1<double, 3>> bodyForces;
160};
161
162//! [Run problem]
167 CHKERR bC();
168 CHKERR OPs();
169 CHKERR tsSolve();
171}
172//! [Run problem]
173
174//! [Set up problem]
178 // Add field
180 constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
181 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
182 CHKERR simple->addDomainField("TAU", L2, base, 1);
183 CHKERR simple->addDomainField("EP", L2, base, size_symm);
184 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
185 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
186 CHKERR simple->setFieldOrder("U", order);
187 CHKERR simple->setFieldOrder("TAU", order - 1);
188 CHKERR simple->setFieldOrder("EP", order - 1);
189 CHKERR simple->setUp();
191}
192//! [Set up problem]
193
194//! [Create common data]
197
198 auto get_command_line_parameters = [&]() {
200 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-scale", &scale, PETSC_NULL);
201 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho", &rho, PETSC_NULL);
202 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
203 &young_modulus, PETSC_NULL);
204 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
205 &poisson_ratio, PETSC_NULL);
206 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening", &H, PETSC_NULL);
207 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening_viscous", &visH,
208 PETSC_NULL);
209 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-yield_stress", &sigmaY,
210 PETSC_NULL);
211 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn", &cn, PETSC_NULL);
212 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-Qinf", &Qinf, PETSC_NULL);
213 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-b_iso", &b_iso, PETSC_NULL);
214
215 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-large_strains",
216 &is_large_strains, PETSC_NULL);
217
218 MOFEM_LOG("EXAMPLE", Sev::inform) << "Young modulus " << young_modulus;
219 MOFEM_LOG("EXAMPLE", Sev::inform) << "Poisson ratio " << poisson_ratio;
220 MOFEM_LOG("EXAMPLE", Sev::inform) << "Yield stress " << sigmaY;
221 MOFEM_LOG("EXAMPLE", Sev::inform) << "Hardening " << H;
222 MOFEM_LOG("EXAMPLE", Sev::inform) << "Viscous hardening " << visH;
223 MOFEM_LOG("EXAMPLE", Sev::inform) << "Saturation yield stress " << Qinf;
224 MOFEM_LOG("EXAMPLE", Sev::inform) << "Saturation exponent " << b_iso;
225 MOFEM_LOG("EXAMPLE", Sev::inform) << "cn " << cn;
226
227 PetscBool is_scale = PETSC_TRUE;
228 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_scale", &is_scale,
229 PETSC_NULL);
230 if (is_scale) {
233 rho *= scale;
234 sigmaY *= scale;
235 H *= scale;
236 Qinf *= scale;
237 visH *= scale;
238
239 MOFEM_LOG("EXAMPLE", Sev::inform)
240 << "Scaled Young modulus " << young_modulus;
241 MOFEM_LOG("EXAMPLE", Sev::inform)
242 << "Scaled Poisson ratio " << poisson_ratio;
243 MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Yield stress " << sigmaY;
244 MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Hardening " << H;
245 MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Viscous hardening " << visH;
246 MOFEM_LOG("EXAMPLE", Sev::inform)
247 << "Scaled Saturation yield stress " << Qinf;
248 }
249
251 };
252
253 auto set_matrial_stiffness = [&]() {
260 const double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
261 const double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
262
263 // Plane stress or when 1, plane strain or 3d
264 const double A = (SPACE_DIM == 2)
265 ? 2 * shear_modulus_G /
266 (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
267 : 1;
268
269 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
270 *commonPlasticDataPtr->mDPtr);
271 auto t_D_axiator = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
272 *commonPlasticDataPtr->mDPtr_Axiator);
273 auto t_D_deviator = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
274 *commonPlasticDataPtr->mDPtr_Deviator);
275
276 constexpr double third = boost::math::constants::third<double>();
277 t_D_axiator(i, j, k, l) = A *
278 (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
279 t_kd(i, j) * t_kd(k, l);
280 t_D_deviator(i, j, k, l) =
281 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
282 t_D(i, j, k, l) = t_D_axiator(i, j, k, l) + t_D_deviator(i, j, k, l);
283
285 };
286
287 commonPlasticDataPtr = boost::make_shared<PlasticOps::CommonData>();
288 constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
289 commonPlasticDataPtr->mDPtr = boost::make_shared<MatrixDouble>();
290 commonPlasticDataPtr->mDPtr->resize(size_symm * size_symm, 1);
291 commonPlasticDataPtr->mDPtr_Axiator = boost::make_shared<MatrixDouble>();
292 commonPlasticDataPtr->mDPtr_Axiator->resize(size_symm * size_symm, 1);
293 commonPlasticDataPtr->mDPtr_Deviator = boost::make_shared<MatrixDouble>();
294 commonPlasticDataPtr->mDPtr_Deviator->resize(size_symm * size_symm, 1);
295
296 commonPlasticDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
297 commonPlasticDataPtr->mStrainPtr = boost::make_shared<MatrixDouble>();
298 commonPlasticDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
299
300 CHKERR get_command_line_parameters();
301 CHKERR set_matrial_stiffness();
302
303 if (is_large_strains) {
304 commonHenckyDataPtr = boost::make_shared<HenckyOps::CommonData>();
305 commonHenckyDataPtr->matGradPtr = commonPlasticDataPtr->mGradPtr;
307 commonHenckyDataPtr->matLogCPlastic =
308 commonPlasticDataPtr->getPlasticStrainPtr();
309 commonPlasticDataPtr->mStrainPtr = commonHenckyDataPtr->getMatLogC();
310 commonPlasticDataPtr->mStressPtr =
311 commonHenckyDataPtr->getMatHenckyStress();
312 }
313
315}
316//! [Create common data]
317
318//! [Boundary condition]
321
323 auto bc_mng = mField.getInterface<BcManager>();
324 auto prb_mng = mField.getInterface<ProblemsManager>();
325
326 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
327 "U", 0, 0);
328 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
329 "U", 1, 1);
330 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
331 "U", 2, 2);
332 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
333 "REMOVE_ALL", "U", 0, 3);
334
335 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_X", "U",
336 0, 0);
337 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_Y", "U",
338 1, 1);
339 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_Z", "U",
340 2, 2);
341 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
342 "U", 0, 3);
343
344 auto &bc_map = bc_mng->getBcMapByBlockName();
345 boundaryMarker = bc_mng->getMergedBlocksMarker(vector<string>{"FIX_"});
346
347 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "REACTION",
348 "U", 0, 3);
349
350 for (auto bc : bc_map)
351 MOFEM_LOG("EXAMPLE", Sev::verbose) << "Marker " << bc.first;
352
353 // OK. We have problem with GMesh, it adding empty characters at the end of
354 // block. So first block is search by regexp. popMarkDOFsOnEntities should
355 // work with regexp.
356 std::string reaction_block_set;
357 for (auto bc : bc_map) {
358 if (bc_mng->checkBlock(bc, "REACTION")) {
359 reaction_block_set = bc.first;
360 break;
361 }
362 }
363
364 if (auto bc = bc_mng->popMarkDOFsOnEntities(reaction_block_set)) {
365 reactionMarker = bc->getBcMarkersPtr();
366
367 // Only take reaction from nodes
368 Range nodes;
369 CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, nodes, true);
370 CHKERR prb_mng->markDofs(simple->getProblemName(), ROW,
371 ProblemsManager::MarkOP::AND, nodes,
373
374 } else {
375 MOFEM_LOG("EXAMPLE", Sev::warning) << "REACTION blockset does not exist";
376 }
377
378 if (!reactionMarker) {
379 MOFEM_LOG("EXAMPLE", Sev::warning) << "REACTION blockset does not exist";
380 }
381
383}
384//! [Boundary condition]
385
386//! [Push operators to pipeline]
389 auto pipeline_mng = mField.getInterface<PipelineManager>();
391 auto bc_mng = mField.getInterface<BcManager>();
392
393 feAxiatorLhs = boost::make_shared<DomainEle>(mField);
394 feAxiatorRhs = boost::make_shared<DomainEle>(mField);
395 auto integration_rule_axiator = [](int, int, int approx_order) {
396 return 2 * (approx_order - 1);
397 };
398 feAxiatorLhs->getRuleHook = integration_rule_axiator;
399 feAxiatorRhs->getRuleHook = integration_rule_axiator;
400
401 auto integration_rule_deviator = [](int o_row, int o_col, int approx_order) {
402 return 2 * (approx_order - 1);
403 };
404 auto integration_rule_bc = [](int, int, int approx_order) {
405 return 2 * approx_order;
406 };
407
408 auto add_domain_base_ops = [&](auto &pipeline) {
410
411 auto det_ptr = boost::make_shared<VectorDouble>();
412 auto jac_ptr = boost::make_shared<MatrixDouble>();
413 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
414 pipeline.push_back(new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
415 pipeline.push_back(
416 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
417 pipeline.push_back(
418 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
419 pipeline.push_back(new OpSetHOWeights(det_ptr));
420
421 pipeline.push_back(new OpCalculateScalarFieldValuesDot(
422 "TAU", commonPlasticDataPtr->getPlasticTauDotPtr()));
424 "EP", commonPlasticDataPtr->getPlasticStrainPtr()));
426 "EP", commonPlasticDataPtr->getPlasticStrainDotPtr()));
428 "U", commonPlasticDataPtr->mGradPtr));
429 pipeline.push_back(new OpCalculateScalarFieldValues(
430 "TAU", commonPlasticDataPtr->getPlasticTauPtr()));
431
433 };
434
435 auto add_domain_stress_ops = [&](auto &pipeline, auto m_D_ptr) {
437
438 if (is_large_strains) {
439
440 if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
441 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
442 "Wrong pointer for grad");
443
444 pipeline.push_back(
446 pipeline.push_back(
448 pipeline.push_back(
450 pipeline.push_back(new OpCalculateHenckyPlasticStress<SPACE_DIM>(
451 "U", commonHenckyDataPtr, m_D_ptr));
452 pipeline.push_back(
454
455 } else {
456 pipeline.push_back(
458 commonPlasticDataPtr->mStrainPtr));
459 pipeline.push_back(
460 new OpPlasticStress("U", commonPlasticDataPtr, m_D_ptr, 1));
461 }
462
463 if (m_D_ptr != commonPlasticDataPtr->mDPtr_Axiator)
464 pipeline.push_back(
466
468 };
469
470 auto add_domain_ops_lhs_mechanical = [&](auto &pipeline, auto m_D_ptr) {
472 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
473
474 if (is_large_strains) {
475 pipeline.push_back(
477 pipeline.push_back(
478 new OpKPiola("U", "U", commonHenckyDataPtr->getMatTangent()));
480 "U", "EP", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
481 } else {
482 pipeline.push_back(new OpKCauchy("U", "U", m_D_ptr));
483 pipeline.push_back(new OpCalculatePlasticInternalForceLhs_dEP(
484 "U", "EP", commonPlasticDataPtr, m_D_ptr));
485 }
486
487 pipeline.push_back(new OpUnSetBc("U"));
489 };
490
491 auto add_domain_ops_rhs_mechanical = [&](auto &pipeline) {
493 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
494
496 const std::string block_name = "BODY_FORCE";
497 if (it->getName().compare(0, block_name.size(), block_name) == 0) {
498 std::vector<double> attr;
499 CHKERR it->getAttributes(attr);
500 if (attr.size() == 3) {
501 bodyForces.push_back(
502 FTensor::Tensor1<double, 3>{attr[0], attr[1], attr[2]});
503 } else {
504 SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
505 "Should be three atributes in BODYFORCE blockset, but is %d",
506 attr.size());
507 }
508 }
509 }
510
511 auto get_body_force = [this](const double, const double, const double) {
512 auto *pipeline_mng = mField.getInterface<PipelineManager>();
515 t_source(i) = 0;
516 auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
517 const auto time = fe_domain_rhs->ts_t;
518 // hardcoded gravity load
519 for (auto &t_b : bodyForces)
520 t_source(i) += (scale * t_b(i)) * time;
521 return t_source;
522 };
523
524 pipeline.push_back(new OpBodyForce("U", get_body_force));
525
526 // Calculate internal forece
527 if (is_large_strains) {
528 pipeline.push_back(new OpInternalForcePiola(
529 "U", commonHenckyDataPtr->getMatFirstPiolaStress()));
530 } else {
531 pipeline.push_back(
532 new OpInternalForceCauchy("U", commonPlasticDataPtr->mStressPtr));
533 }
534
535 pipeline.push_back(new OpUnSetBc("U"));
537 };
538
539 auto add_domain_ops_lhs_constrain = [&](auto &pipeline, auto m_D_ptr) {
541 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
542
543 if (is_large_strains) {
544 pipeline.push_back(
546 pipeline.push_back(new OpCalculateContrainsLhs_LogStrain_dU(
547 "TAU", "U", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
548 pipeline.push_back(new OpCalculatePlasticFlowLhs_LogStrain_dU(
549 "EP", "U", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
550 } else {
551 pipeline.push_back(new OpCalculatePlasticFlowLhs_dU(
552 "EP", "U", commonPlasticDataPtr, m_D_ptr));
553 pipeline.push_back(new OpCalculateContrainsLhs_dU(
554 "TAU", "U", commonPlasticDataPtr, m_D_ptr));
555 }
556
557 pipeline.push_back(new OpCalculatePlasticFlowLhs_dEP(
558 "EP", "EP", commonPlasticDataPtr, m_D_ptr));
559 pipeline.push_back(
561 pipeline.push_back(new OpCalculateContrainsLhs_dEP(
562 "TAU", "EP", commonPlasticDataPtr, m_D_ptr));
563 pipeline.push_back(
565
566 pipeline.push_back(new OpUnSetBc("U"));
568 };
569
570 auto add_domain_ops_rhs_constrain = [&](auto &pipeline) {
572 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
573
574 pipeline.push_back(
576 pipeline.push_back(
578
580 };
581
582 auto add_boundary_ops_lhs_mechanical = [&](auto &pipeline) {
584 auto &bc_map = mField.getInterface<BcManager>()->getBcMapByBlockName();
585 for (auto bc : bc_map) {
586 if (bc_mng->checkBlock(bc, "FIX_")){
587 pipeline.push_back(
588 new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
589 pipeline.push_back(new OpBoundaryMass(
590 "U", "U", [](double, double, double) { return 1.; },
591 bc.second->getBcEntsPtr()));
592 pipeline.push_back(new OpUnSetBc("U"));
593 }
594 }
596 };
597
598 auto add_boundary_ops_rhs_mechanical = [&](auto &pipeline) {
600
601 auto get_time = [&](double, double, double) {
602 auto *pipeline_mng = mField.getInterface<PipelineManager>();
603 auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
604 return fe_domain_rhs->ts_t;
605 };
606
607 auto get_time_scaled = [&](double, double, double) {
608 auto *pipeline_mng = mField.getInterface<PipelineManager>();
609 auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
610 return fe_domain_rhs->ts_t * scale;
611 };
612
613 auto get_minus_time = [&](double, double, double) {
614 auto *pipeline_mng = mField.getInterface<PipelineManager>();
615 auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
616 return -fe_domain_rhs->ts_t;
617 };
618
619 auto time_scaled = [&](double, double, double) {
620 auto *pipeline_mng = mField.getInterface<PipelineManager>();
621 auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
622 return -fe_domain_rhs->ts_t;
623 };
624
625 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
626
628 if (it->getName().compare(0, 5, "FORCE") == 0) {
629 Range force_edges;
630 std::vector<double> attr_vec;
631 CHKERR it->getMeshsetIdEntitiesByDimension(
632 mField.get_moab(), SPACE_DIM - 1, force_edges, true);
633 it->getAttributes(attr_vec);
634 if (attr_vec.size() < SPACE_DIM)
635 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
636 "Wrong size of boundary attributes vector. Set right block "
637 "size attributes.");
638 auto force_vec_ptr = boost::make_shared<MatrixDouble>(SPACE_DIM, 1);
639 std::copy(&attr_vec[0], &attr_vec[SPACE_DIM],
640 force_vec_ptr->data().begin());
641 pipeline.push_back(
642 new OpBoundaryVec("U", force_vec_ptr, get_time_scaled,
643 boost::make_shared<Range>(force_edges)));
644 }
645 }
646
647 pipeline.push_back(new OpUnSetBc("U"));
648
649 auto u_mat_ptr = boost::make_shared<MatrixDouble>();
650 pipeline.push_back(
651 new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_mat_ptr));
652
653 for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName()) {
654 if (bc_mng->checkBlock(bc, "FIX_")) {
655 pipeline.push_back(
656 new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
657 auto attr_vec = boost::make_shared<MatrixDouble>(SPACE_DIM, 1);
658 attr_vec->clear();
659 if (bc.second->bcAttributes.size() < SPACE_DIM)
660 SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
661 "Wrong size of boundary attributes vector. Set right block "
662 "size attributes. Size of attributes %d",
663 bc.second->bcAttributes.size());
664 std::copy(&bc.second->bcAttributes[0],
665 &bc.second->bcAttributes[SPACE_DIM],
666 attr_vec->data().begin());
667
668 pipeline.push_back(new OpBoundaryVec("U", attr_vec, time_scaled,
669 bc.second->getBcEntsPtr()));
670 pipeline.push_back(new OpBoundaryInternal(
671 "U", u_mat_ptr, [](double, double, double) { return 1.; },
672 bc.second->getBcEntsPtr()));
673
674 pipeline.push_back(new OpUnSetBc("U"));
675 }
676 }
677
679 };
680
681 // Axiator
682 CHKERR add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
683 CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainLhsPipeline(),
684 commonPlasticDataPtr->mDPtr_Deviator);
685 CHKERR add_domain_ops_lhs_mechanical(pipeline_mng->getOpDomainLhsPipeline(),
686 commonPlasticDataPtr->mDPtr_Deviator);
687 CHKERR add_domain_ops_lhs_constrain(pipeline_mng->getOpDomainLhsPipeline(),
688 commonPlasticDataPtr->mDPtr_Deviator);
689 CHKERR add_boundary_ops_lhs_mechanical(
690 pipeline_mng->getOpBoundaryLhsPipeline());
691
692 CHKERR add_domain_base_ops(feAxiatorLhs->getOpPtrVector());
693 CHKERR add_domain_stress_ops(feAxiatorLhs->getOpPtrVector(),
694 commonPlasticDataPtr->mDPtr_Axiator);
695 CHKERR add_domain_ops_lhs_mechanical(feAxiatorLhs->getOpPtrVector(),
696 commonPlasticDataPtr->mDPtr_Axiator);
697
698 CHKERR add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
699 CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainRhsPipeline(),
700 commonPlasticDataPtr->mDPtr_Deviator);
701 CHKERR add_domain_ops_rhs_mechanical(pipeline_mng->getOpDomainRhsPipeline());
702 CHKERR add_domain_ops_rhs_constrain(pipeline_mng->getOpDomainRhsPipeline());
703 CHKERR add_boundary_ops_rhs_mechanical(
704 pipeline_mng->getOpBoundaryRhsPipeline());
705
706 CHKERR add_domain_base_ops(feAxiatorRhs->getOpPtrVector());
707 CHKERR add_domain_stress_ops(feAxiatorRhs->getOpPtrVector(),
708 commonPlasticDataPtr->mDPtr_Axiator);
709 CHKERR add_domain_ops_rhs_mechanical(feAxiatorRhs->getOpPtrVector());
710
711 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule_deviator);
712 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule_deviator);
713
714 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule_bc);
715 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule_bc);
716
717 auto create_reaction_pipeline = [&](auto &pipeline) {
719
720 if (reactionMarker) {
721
722 auto det_ptr = boost::make_shared<VectorDouble>();
723 auto jac_ptr = boost::make_shared<MatrixDouble>();
724 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
725 pipeline.push_back(new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
726 pipeline.push_back(
727 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
728 pipeline.push_back(
729 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
730 pipeline.push_back(new OpSetHOWeights(det_ptr));
731
732 pipeline.push_back(
734 "U", commonPlasticDataPtr->mGradPtr));
736 "EP", commonPlasticDataPtr->getPlasticStrainPtr()));
737
738 if (is_large_strains) {
739
740 if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
741 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
742 "Wrong pointer for grad");
743
744 pipeline.push_back(
746 pipeline.push_back(
748 pipeline.push_back(
750 pipeline.push_back(new OpCalculateHenckyPlasticStress<SPACE_DIM>(
752 pipeline.push_back(
754
755 } else {
756 pipeline.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
757 "U", commonPlasticDataPtr->mGradPtr,
758 commonPlasticDataPtr->mStrainPtr));
759 pipeline.push_back(new OpPlasticStress(
761 }
762
763 pipeline.push_back(new OpSetBc("U", false, reactionMarker));
764 // Calculate internal forece
765 if (is_large_strains) {
766 pipeline.push_back(new OpInternalForcePiola(
767 "U", commonHenckyDataPtr->getMatFirstPiolaStress()));
768 } else {
769 pipeline.push_back(
770 new OpInternalForceCauchy("U", commonPlasticDataPtr->mStressPtr));
771 }
772 pipeline.push_back(new OpUnSetBc("U"));
773 }
774
776 };
777
778 reactionFe = boost::make_shared<DomainEle>(mField);
779 reactionFe->getRuleHook = integration_rule_deviator;
780
781 CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
782
784}
785//! [Push operators to pipeline]
786
787//! [Solve]
790
793 ISManager *is_manager = mField.getInterface<ISManager>();
794
795 auto snes_ctx_ptr = smartGetDMSnesCtx(simple->getDM());
796
797 auto set_section_monitor = [&](auto solver) {
799 SNES snes;
800 CHKERR TSGetSNES(solver, &snes);
801 CHKERR SNESMonitorSet(snes,
802 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
804 (void *)(snes_ctx_ptr.get()), nullptr);
806 };
807
808 auto create_post_process_element = [&]() {
810 postProcFe = boost::make_shared<PostProcEle>(mField);
811 postProcFe->generateReferenceElementMesh();
812
813 auto det_ptr = boost::make_shared<VectorDouble>();
814 auto jac_ptr = boost::make_shared<MatrixDouble>();
815 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
816 postProcFe->getOpPtrVector().push_back(
817 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
818 postProcFe->getOpPtrVector().push_back(
819 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
820 postProcFe->getOpPtrVector().push_back(
821 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
822
823 postProcFe->getOpPtrVector().push_back(
825 "U", commonPlasticDataPtr->mGradPtr));
826 postProcFe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
827 "TAU", commonPlasticDataPtr->getPlasticTauPtr()));
828 postProcFe->getOpPtrVector().push_back(
830 "EP", commonPlasticDataPtr->getPlasticStrainPtr()));
831
832 if (is_large_strains) {
833
834 if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
835 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
836 "Wrong pointer for grad");
837
838 postProcFe->getOpPtrVector().push_back(
840 postProcFe->getOpPtrVector().push_back(
842 postProcFe->getOpPtrVector().push_back(
844 postProcFe->getOpPtrVector().push_back(
847 postProcFe->getOpPtrVector().push_back(
849 postProcFe->getOpPtrVector().push_back(new OpPostProcHencky<SPACE_DIM>(
850 "U", postProcFe->postProcMesh, postProcFe->mapGaussPts,
852
853 } else {
854 postProcFe->getOpPtrVector().push_back(
856 commonPlasticDataPtr->mStrainPtr));
857 postProcFe->getOpPtrVector().push_back(new OpPlasticStress(
859 postProcFe->getOpPtrVector().push_back(
861 "U", postProcFe->postProcMesh, postProcFe->mapGaussPts,
862 commonPlasticDataPtr->mStrainPtr,
863 commonPlasticDataPtr->mStressPtr));
864 }
865
866 postProcFe->getOpPtrVector().push_back(
868 postProcFe->getOpPtrVector().push_back(
869 new OpPostProcPlastic("U", postProcFe->postProcMesh,
870 postProcFe->mapGaussPts, commonPlasticDataPtr));
871 postProcFe->addFieldValuesPostProc("U");
873 };
874
875 auto scatter_create = [&](auto D, auto coeff) {
877 CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
878 ROW, "U", coeff, coeff, is);
879 int loc_size;
880 CHKERR ISGetLocalSize(is, &loc_size);
881 Vec v;
882 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
883 VecScatter scatter;
884 CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
885 return std::make_tuple(SmartPetscObj<Vec>(v),
887 };
888
889 auto set_time_monitor = [&](auto dm, auto solver) {
891 boost::shared_ptr<Monitor> monitor_ptr(new Monitor(
893 boost::shared_ptr<ForcesAndSourcesCore> null;
894 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
895 monitor_ptr, null, null);
897 };
898
899 auto set_fieldsplit_preconditioner = [&](auto solver) {
901
902 SNES snes;
903 CHKERR TSGetSNES(solver, &snes);
904 KSP ksp;
905 CHKERR SNESGetKSP(snes, &ksp);
906 PC pc;
907 CHKERR KSPGetPC(ksp, &pc);
908 PetscBool is_pcfs = PETSC_FALSE;
909 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
910
911 // Setup fieldsplit (block) solver - optional: yes/no
912 if (is_pcfs == PETSC_TRUE) {
913
914 auto bc_mng = mField.getInterface<BcManager>();
915 auto name_prb = simple->getProblemName();
916 auto is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_X", "U", 0, 0);
917 is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_Y", "U", 1, 1, is_all_bc);
918 is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_Z", "U", 2, 2, is_all_bc);
919 is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_ALL", "U", 0, 2, is_all_bc);
920
921 int is_all_bc_size;
922 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
923 MOFEM_LOG("EXAMPLE", Sev::inform)
924 << "Field split block size " << is_all_bc_size;
925
926 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
927 is_all_bc); // boundary block
928 }
929
931 };
932
933 auto dm = simple->getDM();
934 auto D = smartCreateDMVector(dm);
935
936 boost::shared_ptr<FEMethod> null;
937 CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), feAxiatorLhs,
938 null, null);
939 CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(), feAxiatorRhs,
940 null, null);
941
942 CHKERR create_post_process_element();
943 uXScatter = scatter_create(D, 0);
944 uYScatter = scatter_create(D, 1);
945 if (SPACE_DIM == 3)
946 uZScatter = scatter_create(D, 2);
947
948 auto solver = pipeline_mng->createTSIM();
949
950 CHKERR TSSetSolution(solver, D);
951 CHKERR set_section_monitor(solver);
952 CHKERR set_time_monitor(dm, solver);
953 CHKERR TSSetSolution(solver, D);
954 CHKERR TSSetFromOptions(solver);
955 CHKERR set_fieldsplit_preconditioner(solver);
956 CHKERR TSSetUp(solver);
957 CHKERR TSSolve(solver, NULL);
958
959 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
960 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
961 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
962
964}
965//! [Solve]
966
967static char help[] = "...\n\n";
968
969int main(int argc, char *argv[]) {
970
971 // Initialisation of MoFEM/PETSc and MOAB data structures
972 const char param_file[] = "param_file.petsc";
974
975 // Add logging channel for example
976 auto core_log = logging::core::get();
977 core_log->add_sink(
978 LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
979 LogManager::setLog("EXAMPLE");
980 MOFEM_LOG_TAG("EXAMPLE", "example");
981
982 try {
983
984 //! [Register MoFEM discrete manager in PETSc]
985 DMType dm_name = "DMMOFEM";
986 CHKERR DMRegister_MoFEM(dm_name);
987 //! [Register MoFEM discrete manager in PETSc
988
989 //! [Create MoAB]
990 moab::Core mb_instance; ///< mesh database
991 moab::Interface &moab = mb_instance; ///< mesh database interface
992 //! [Create MoAB]
993
994 //! [Create MoFEM]
995 MoFEM::Core core(moab); ///< finite element database
996 MoFEM::Interface &m_field = core; ///< finite element database insterface
997 //! [Create MoFEM]
998
999 //! [Load mesh]
1000 Simple *simple = m_field.getInterface<Simple>();
1001 CHKERR simple->getOptions();
1002 CHKERR simple->loadFile();
1003 //! [Load mesh]
1004
1005 //! [Example]
1006 Example ex(m_field);
1007 CHKERR ex.runProblem();
1008 //! [Example]
1009 }
1011
1013}
std::string param_file
constexpr double third
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEM::FaceElementForcesAndSourcesCore FaceEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
@ ROW
Definition: definitions.h:136
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
FieldApproximationBase
approximation base
Definition: definitions.h:71
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ BLOCKSET
Definition: definitions.h:161
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
@ MOFEM_INVALID_DATA
Definition: definitions.h:49
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
constexpr double shear_modulus_G
Definition: elastic.cpp:66
constexpr double bulk_modulus_K
Definition: elastic.cpp:65
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
constexpr auto t_kd
void temp(int x, int y=10)
Definition: simple.cpp:4
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
Definition: DMMMoFEM.cpp:759
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:482
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
Definition: DMMMoFEM.cpp:812
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:977
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:311
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:342
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
FTensor::Index< 'i', SPACE_DIM > i
double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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: DMMMoFEM.cpp:1015
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:236
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)
CoreTmp< 0 > Core
Definition: Core.hpp:1096
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
auto smartGetDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:999
const double D
diffusivity
double A
double young_modulus
Definition: plastic.cpp:105
int main(int argc, char *argv[])
Definition: plastic.cpp:969
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: plastic.cpp:88
EntitiesFieldData::EntData EntData
Definition: plastic.cpp:55
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: plastic.cpp:86
double Qinf
Definition: plastic.cpp:112
static char help[]
[Solve]
Definition: plastic.cpp:967
double rho
Definition: plastic.cpp:107
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: plastic.cpp:74
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:24
constexpr int SPACE_DIM
Definition: plastic.cpp:51
double visH
Definition: plastic.cpp:110
double poisson_ratio
Definition: plastic.cpp:106
double scale
Definition: plastic.cpp:103
long double hardening(long double tau, double temp)
Definition: plastic.cpp:116
double H
Definition: plastic.cpp:109
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: plastic.cpp:95
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Body force]
Definition: plastic.cpp:72
long double hardening_dtau(long double tau, double temp)
Definition: plastic.cpp:120
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
Definition: plastic.cpp:97
int order
Definition: plastic.cpp:114
double b_iso
Definition: plastic.cpp:113
PetscBool is_large_strains
Definition: plastic.cpp:101
double sigmaY
Definition: plastic.cpp:108
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:93
double cn
Definition: plastic.cpp:111
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
[Body force]
Definition: plastic.cpp:67
constexpr EntityType boundary_ent
Definition: plastic.cpp:54
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:41
[Example]
Definition: plastic.cpp:130
boost::shared_ptr< DomainEle > feAxiatorRhs
Definition: plastic.cpp:149
boost::shared_ptr< DomainEle > reactionFe
Definition: plastic.cpp:148
boost::shared_ptr< std::vector< unsigned char > > reactionMarker
Definition: plastic.cpp:157
boost::shared_ptr< PostProcEle > postProcFe
Definition: plastic.cpp:147
boost::shared_ptr< DomainEle > feAxiatorLhs
Definition: plastic.cpp:150
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
Definition: plastic.cpp:788
FieldApproximationBase base
Definition: plot_base.cpp:83
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: plastic.cpp:153
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:195
std::vector< FTensor::Tensor1< double, 3 > > bodyForces
Definition: plastic.cpp:159
Example(MoFEM::Interface &m_field)
Definition: plastic.cpp:132
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:387
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:163
MoFEM::Interface & mField
Definition: plastic.cpp:137
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:154
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: plastic.cpp:156
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:175
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:319
boost::shared_ptr< PlasticOps::CommonData > commonPlasticDataPtr
Definition: plastic.cpp:145
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:152
boost::shared_ptr< HenckyOps::CommonData > commonHenckyDataPtr
Definition: plastic.cpp:146
Definition: HenckyOps.hpp:512
Simple interface for fast problem set-up.
Definition: BcManager.hpp:27
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Definition: BcManager.hpp:129
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:33
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.
Scale base functions by inverses of measure of element.
Set indices on entities on finite element.
Set inverse jacobian to base functions.
Set inverse jacobian to base functions.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:35
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: PlasticOps.hpp:279
Postprocess on face.
Post processing.
[Class definition]