v0.13.1
thermo_plastic.cpp
Go to the documentation of this file.
1/**
2 * \file thermo_plastic.cpp
3 * \example thermo_plastic.cpp
4 *
5 * Thermo plasticity
6 *
7 */
8
9
10
11#ifndef EXECUTABLE_DIMENSION
12#define EXECUTABLE_DIMENSION 3
13#endif
14
15#include <MoFEM.hpp>
16#include <MatrixFunction.hpp>
17#include <IntegrationRules.hpp>
18
19using namespace MoFEM;
20
21template <int DIM> struct ElementsAndOps {};
22
23template <> struct ElementsAndOps<2> {
26};
27
28template <> struct ElementsAndOps<3> {
31};
32
33constexpr int SPACE_DIM =
34 EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
35
36constexpr EntityType boundary_ent = SPACE_DIM == 3 ? MBTRI : MBEDGE;
43
46
47//! [Body force]
49 GAUSS>::OpSource<1, SPACE_DIM>;
50//! [Body force]
51
52//! [Only used with Hooke equation (linear material model)]
54 GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
57//! [Only used with Hooke equation (linear material model)]
58
59//! [Only used for dynamics]
64//! [Only used for dynamics]
65
66//! [Only used with Hencky/nonlinear material]
68 GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
71//! [Only used with Hencky/nonlinear material]
72
73//! [Essential boundary conditions]
80//! [Essential boundary conditions]
82
83// Thermal operators
84/**
85 * @brief Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)
86 *
87 */
90
91/**
92 * @brief Integrate Lhs div of base of flux time base of temperature (FLUX x T)
93 * and transpose of it, i.e. (T x FLAX)
94 *
95 */
97 GAUSS>::OpMixDivTimesScalar<SPACE_DIM>;
98
99/**
100 * @brief Integrate Lhs base of temerature times (heat capacity) times base of
101 * temperature (T x T)
102 *
103 */
106
107/**
108 * @brief Integrating Rhs flux base (1/k) flux (FLUX)
109 */
111 GAUSS>::OpBaseTimesVector<3, 3, 1>;
112
113/**
114 * @brief Integrate Rhs div flux base times temperature (T)
115 *
116 */
118 GAUSS>::OpMixDivTimesU<3, 1, 2>;
119
120/**
121 * @brief Integrate Rhs base of temerature time heat capacity times heat rate
122 * (T)
123 *
124 */
127
128/**
129 * @brief Integrate Rhs base of temerature times divergenc of flux (T)
130 *
131 */
133
135 GAUSS>::OpSource<1, 1>;
138
139double scale = 1.;
140
141double young_modulus = 206913;
142double poisson_ratio = 0.29;
143double coeff_expansion = 10e-6;
144double ref_temp = 0.0;
145double rho = 0;
146double sigmaY = 450;
147double H = 129;
148double visH = 0;
149double cn = 1;
150double Qinf = 0; // 265;
151double b_iso = 16.93;
153 16.2; // Force / (time temerature ) or Power /
154 // (length temperature) // Time unit is hour. force unit kN
156 5961.6; // length^2/(time^2 temerature) // length is millimeter time is hour
157double omega_0 = 2e-3;
158double omega_h = 2e-3;
159double omega_inf = 0;
161int order = 2;
162
164double amplitude_cycle = 0.5;
165double amplitude_shift = 0.5;
166double phase_shift = 0.8;
167
168inline long double hardening(long double tau, double temp) {
169 return H * tau * (1 - omega_h * temp) +
170 Qinf * (1. - std::exp(-b_iso * tau)) * (1 - omega_inf * temp) +
171 sigmaY * (1 - omega_0 * temp);
172}
173
174inline long double hardening_dtau(long double tau, double temp) {
175 return H * (1 - omega_h * temp) +
176 Qinf * b_iso * std::exp(-b_iso * tau) * (1 - omega_inf * temp);
177}
178
179inline long double hardening_dtemp(long double tau, double temp) {
180 return -H * tau * omega_h - Qinf * (1. - std::exp(-b_iso * tau)) * omega_inf -
181 sigmaY * omega_0;
182}
183
184#include <HenckyOps.hpp>
185#include <PlasticOps.hpp>
186using namespace PlasticOps;
187#include <PlasticThermalOps.hpp>
188
189using namespace HenckyOps;
190struct Example {
191
192 Example(MoFEM::Interface &m_field) : mField(m_field) {}
193
195
196private:
198
204
205 boost::shared_ptr<PlasticThermalOps::CommonData> commonPlasticDataPtr;
206 boost::shared_ptr<HenckyOps::CommonData> commonHenckyDataPtr;
207 boost::shared_ptr<PostProcEle> postProcFe;
208 boost::shared_ptr<DomainEle> reactionFe;
209 boost::shared_ptr<DomainEle> feAxiatorRhs;
210 boost::shared_ptr<DomainEle> feAxiatorLhs;
211 boost::shared_ptr<DomainEle> feThermalRhs;
212 boost::shared_ptr<DomainEle> feThermalLhs;
213
214 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
215 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
216 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
217
218 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
219 boost::shared_ptr<std::vector<unsigned char>> reactionMarker;
220
221 std::vector<FTensor::Tensor1<double, 3>> bodyForces;
222
223 struct BcTempFun {
224 BcTempFun(double v, FEMethod &fe) : valTemp(v), fE(fe) {}
225 double operator()(const double, const double, const double) {
226 return -valTemp; // * fE.ts_t;
227 }
228
229 private:
230 double valTemp;
232 };
233 std::vector<BcTempFun> bcTemperatureFunctions;
234};
235
236//! [Run problem]
241 CHKERR bC();
242 CHKERR OPs();
243 CHKERR tsSolve();
245}
246//! [Run problem]
247
248//! [Set up problem]
252 // Add field
254 constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
255 // Mechanical fields
256 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
257 CHKERR simple->addDomainField("TAU", L2, base, 1);
258 CHKERR simple->addDomainField("EP", L2, base, size_symm);
259 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
260 // Temperature
261 const auto flux_space = (SPACE_DIM == 2) ? HCURL : HDIV;
262 CHKERR simple->addDomainField("T", L2, AINSWORTH_LEGENDRE_BASE, 1);
263 CHKERR simple->addDomainField("FLUX", flux_space, DEMKOWICZ_JACOBI_BASE, 1);
264 CHKERR simple->addBoundaryField("FLUX", flux_space, DEMKOWICZ_JACOBI_BASE, 1);
265
266 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
267 CHKERR simple->setFieldOrder("U", order);
268 CHKERR simple->setFieldOrder("TAU", order - 1);
269 CHKERR simple->setFieldOrder("EP", order - 1);
270 CHKERR simple->setFieldOrder("FLUX", order);
271 CHKERR simple->setFieldOrder("T", order - 1);
272 CHKERR simple->setUp();
274}
275//! [Set up problem]
276
277//! [Create common data]
280
281 auto get_command_line_parameters = [&]() {
283 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-scale", &scale, PETSC_NULL);
284 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho", &rho, PETSC_NULL);
285 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
286 &young_modulus, PETSC_NULL);
287 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
288 &poisson_ratio, PETSC_NULL);
289 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-coeff_expansion",
290 &coeff_expansion, PETSC_NULL);
291 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-ref_temp", &ref_temp,
292 PETSC_NULL);
293 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening", &H, PETSC_NULL);
294 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening_viscous", &visH,
295 PETSC_NULL);
296 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-yield_stress", &sigmaY,
297 PETSC_NULL);
298 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn", &cn, PETSC_NULL);
299 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-Qinf", &Qinf, PETSC_NULL);
300 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-b_iso", &b_iso, PETSC_NULL);
301
302 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-capacity", &heat_capacity,
303 PETSC_NULL);
304 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-conductivity",
305 &heat_conductivity, PETSC_NULL);
306 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-omega_0", &omega_0,
307 PETSC_NULL);
308 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-omega_h", &omega_h,
309 PETSC_NULL);
310 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-omega_inf", &omega_inf,
311 PETSC_NULL);
312 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-fraction_of_dissipation",
313 &fraction_of_dissipation, PETSC_NULL);
314 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-number_of_cycles_in_total_time",
315 &number_of_cycles_in_total_time, PETSC_NULL);
316 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-amplitude_cycle",
317 &amplitude_cycle, PETSC_NULL);
318 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-phase_shift", &phase_shift,
319 PETSC_NULL);
320 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-amplitude_shift",
321 &amplitude_shift, PETSC_NULL);
322 MOFEM_LOG("EXAMPLE", Sev::inform) << "Young modulus " << young_modulus;
323 MOFEM_LOG("EXAMPLE", Sev::inform) << "Poisson ratio " << poisson_ratio;
324 MOFEM_LOG("EXAMPLE", Sev::inform) << "Coeff_expansion " << coeff_expansion;
325 MOFEM_LOG("EXAMPLE", Sev::inform) << "Reference_temperature " << ref_temp;
326 MOFEM_LOG("EXAMPLE", Sev::inform) << "Yield stress " << sigmaY;
327 MOFEM_LOG("EXAMPLE", Sev::inform) << "Hardening " << H;
328 MOFEM_LOG("EXAMPLE", Sev::inform) << "Viscous hardening " << visH;
329 MOFEM_LOG("EXAMPLE", Sev::inform) << "Saturation yield stress " << Qinf;
330 MOFEM_LOG("EXAMPLE", Sev::inform) << "Saturation exponent " << b_iso;
331 MOFEM_LOG("EXAMPLE", Sev::inform) << "cn " << cn;
332 MOFEM_LOG("EXAMPLE", Sev::inform)
333 << "fraction_of_dissipation " << fraction_of_dissipation;
334
335 PetscBool is_scale = PETSC_TRUE;
336 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_scale", &is_scale,
337 PETSC_NULL);
338 if (is_scale) {
341 rho *= scale;
342 sigmaY *= scale;
343 H *= scale;
344 Qinf *= scale;
345 visH *= scale;
347
348 MOFEM_LOG("EXAMPLE", Sev::inform)
349 << "Scaled Young modulus " << young_modulus;
350 MOFEM_LOG("EXAMPLE", Sev::inform)
351 << "Scaled Poisson ratio " << poisson_ratio;
352 MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Yield stress " << sigmaY;
353 MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Hardening " << H;
354 MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Viscous hardening " << visH;
355 MOFEM_LOG("EXAMPLE", Sev::inform)
356 << "Scaled Saturation yield stress " << Qinf;
357 MOFEM_LOG("EXAMPLE", Sev::inform)
358 << "Scaled heat capacity " << heat_capacity;
359 }
360
362 };
363
364 auto set_matrial_stiffness = [&]() {
371 const double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
372 const double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
373
374 // Plane stress or when 1, plane strain or 3d
375 const double A = (SPACE_DIM == 2)
376 ? 2 * shear_modulus_G /
377 (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
378 : 1;
379
380 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
381 *commonPlasticDataPtr->mDPtr);
382 auto t_D_axiator = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
383 *commonPlasticDataPtr->mDPtr_Axiator);
384 auto t_D_deviator = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
385 *commonPlasticDataPtr->mDPtr_Deviator);
386
387 constexpr double third = boost::math::constants::third<double>();
388 t_D_axiator(i, j, k, l) = A *
389 (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
390 t_kd(i, j) * t_kd(k, l);
391 t_D_deviator(i, j, k, l) =
392 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
393 t_D(i, j, k, l) = t_D_axiator(i, j, k, l) + t_D_deviator(i, j, k, l);
394
396 };
397
398 commonPlasticDataPtr = boost::make_shared<PlasticThermalOps::CommonData>();
399 constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
400 commonPlasticDataPtr->mDPtr = boost::make_shared<MatrixDouble>();
401 commonPlasticDataPtr->mDPtr->resize(size_symm * size_symm, 1);
402 commonPlasticDataPtr->mDPtr_Axiator = boost::make_shared<MatrixDouble>();
403 commonPlasticDataPtr->mDPtr_Axiator->resize(size_symm * size_symm, 1);
404 commonPlasticDataPtr->mDPtr_Deviator = boost::make_shared<MatrixDouble>();
405 commonPlasticDataPtr->mDPtr_Deviator->resize(size_symm * size_symm, 1);
406
407 commonPlasticDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
408 commonPlasticDataPtr->mStrainPtr = boost::make_shared<MatrixDouble>();
409 commonPlasticDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
410
411 CHKERR get_command_line_parameters();
412 CHKERR set_matrial_stiffness();
413
415}
416//! [Create common data]
417
418//! [Boundary condition]
421
423 auto bc_mng = mField.getInterface<BcManager>();
424 auto prb_mng = mField.getInterface<ProblemsManager>();
425
426 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
427 "U", 0, 0);
428 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
429 "U", 1, 1);
430 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
431 "U", 2, 2);
432 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
433 "REMOVE_ALL", "U", 0, 3);
434 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
435 "ZERO_FLUX", "FLUX", 0, 1);
436
437 PetscBool zero_fix_skin_flux = PETSC_TRUE;
438 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-fix_skin_flux",
439 &zero_fix_skin_flux, PETSC_NULL);
440 if (zero_fix_skin_flux) {
441 Range faces;
442 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, faces);
443 Skinner skin(&mField.get_moab());
444 Range skin_ents;
445 CHKERR skin.find_skin(0, faces, false, skin_ents);
446 ParallelComm *pcomm =
447 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
448 if (pcomm == NULL)
449 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
450 "Communicator not created");
451 CHKERR pcomm->filter_pstatus(skin_ents,
452 PSTATUS_SHARED | PSTATUS_MULTISHARED,
453 PSTATUS_NOT, -1, &skin_ents);
454 CHKERR prb_mng->removeDofsOnEntities(simple->getProblemName(), "FLUX",
455 skin_ents, 0, 1);
456 }
457
458 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_X", "U",
459 0, 0);
460 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_Y", "U",
461 1, 1);
462 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_Z", "U",
463 2, 2);
464 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
465 "U", 0, 3);
466
467 auto &bc_map = bc_mng->getBcMapByBlockName();
468 boundaryMarker = bc_mng->getMergedBlocksMarker(vector<string>{"FIX_"});
469
470 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "REACTION",
471 "U", 0, 3);
472 for (auto b : bc_map) {
473 MOFEM_LOG("EXAMPLE", Sev::verbose) << "Marker " << b.first;
474 }
475
476 // OK. We have problem with GMesh, it adding empty characters at the end of
477 // block. So first block is search by regexp. popMarkDOFsOnEntities should
478 // work with regexp.
479 std::string reaction_block_set;
480 for (auto bc : bc_map) {
481 if (bc_mng->checkBlock(bc, "REACTION")) {
482 reaction_block_set = bc.first;
483 break;
484 }
485 }
486
487 if (auto bc = bc_mng->popMarkDOFsOnEntities(reaction_block_set)) {
488 reactionMarker = bc->getBcMarkersPtr();
489
490 // Only take reaction from nodes
491 Range nodes;
492 CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, nodes, true);
493 CHKERR prb_mng->markDofs(simple->getProblemName(), ROW,
494 ProblemsManager::MarkOP::AND, nodes,
496
497 } else {
498 MOFEM_LOG("EXAMPLE", Sev::warning) << "REACTION blockset does not exist";
499 }
500
501 if (!reactionMarker) {
502 MOFEM_LOG("EXAMPLE", Sev::warning) << "REACTION blockset does not exist";
503 }
504
506}
507//! [Boundary condition]
508
509//! [Push operators to pipeline]
512 auto pipeline_mng = mField.getInterface<PipelineManager>();
514 auto bc_mng = mField.getInterface<BcManager>();
515
516 auto integration_rule_deviator = [](int o_row, int o_col, int approx_order) {
517 return 2 * (approx_order - 1);
518 };
519 auto integration_rule_bc = [](int, int, int approx_order) {
520 return 2 * approx_order;
521 };
522
523 feAxiatorRhs = boost::make_shared<DomainEle>(mField);
524 feAxiatorLhs = boost::make_shared<DomainEle>(mField);
525 auto integration_rule_axiator = [](int, int, int approx_order) {
526 return 2 * (approx_order - 1);
527 };
528 feAxiatorRhs->getRuleHook = integration_rule_axiator;
529 feAxiatorLhs->getRuleHook = integration_rule_axiator;
530
531 feThermalRhs = boost::make_shared<DomainEle>(mField);
532 feThermalLhs = boost::make_shared<DomainEle>(mField);
533 auto integration_rule_thermal = [](int, int, int approx_order) {
534 return 2 * approx_order;
535 };
536 feThermalRhs->getRuleHook = integration_rule_thermal;
537 feThermalLhs->getRuleHook = integration_rule_thermal;
538
539 auto add_domain_base_ops = [&](auto &pipeline) {
541
542 auto det_ptr = boost::make_shared<VectorDouble>();
543 auto jac_ptr = boost::make_shared<MatrixDouble>();
544 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
545 pipeline.push_back(new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
546 pipeline.push_back(
547 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
548 pipeline.push_back(
549 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
550
551 if (SPACE_DIM == 2) {
552 pipeline.push_back(new OpMakeHdivFromHcurl());
553 pipeline.push_back(new OpSetContravariantPiolaTransformOnFace2D(jac_ptr));
554 pipeline.push_back(new OpSetInvJacHcurlFace(inv_jac_ptr));
555 pipeline.push_back(new OpSetInvJacL2ForFace(inv_jac_ptr));
556 }
557
558 pipeline.push_back(new OpCalculateScalarFieldValuesDot(
559 "TAU", commonPlasticDataPtr->getPlasticTauDotPtr()));
561 "EP", commonPlasticDataPtr->getPlasticStrainPtr()));
563 "EP", commonPlasticDataPtr->getPlasticStrainDotPtr()));
565 "U", commonPlasticDataPtr->mGradPtr));
566 pipeline.push_back(new OpCalculateScalarFieldValues(
567 "TAU", commonPlasticDataPtr->getPlasticTauPtr()));
568
569 pipeline.push_back(new OpCalculateScalarFieldValues(
570 "T", commonPlasticDataPtr->getTempValPtr()));
571 pipeline.push_back(new OpCalculateScalarFieldValuesDot(
572 "T", commonPlasticDataPtr->getTempValDotPtr()));
574 "FLUX", commonPlasticDataPtr->getTempDivFluxPtr()));
575 pipeline.push_back(new OpCalculateHVecVectorField<3>(
576 "FLUX", commonPlasticDataPtr->getTempFluxValPtr()));
577
579 };
580
581 auto add_domain_stress_ops = [&](auto &pipeline, auto m_D_ptr) {
583
584 pipeline.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
585 "U", commonPlasticDataPtr->mGradPtr, commonPlasticDataPtr->mStrainPtr));
586 pipeline.push_back(new PlasticThermalOps::OpPlasticStressThermal(
587 "U", commonPlasticDataPtr, m_D_ptr, 1));
588
589 if (m_D_ptr != commonPlasticDataPtr->mDPtr_Axiator)
590 pipeline.push_back(
592
594 };
595
596 auto add_domain_ops_lhs_mechanical = [&](auto &pipeline, auto m_D_ptr) {
598 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
599
600 pipeline.push_back(new OpKCauchy("U", "U", m_D_ptr));
602 "U", "T", commonPlasticDataPtr, m_D_ptr));
603 pipeline.push_back(new OpCalculatePlasticInternalForceLhs_dEP(
604 "U", "EP", commonPlasticDataPtr, m_D_ptr));
605
606 pipeline.push_back(new OpUnSetBc("U"));
608 };
609
610 auto add_domain_ops_rhs_mechanical = [&](auto &pipeline) {
612 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
613
615 const std::string block_name = "BODY_FORCE";
616 if (it->getName().compare(0, block_name.size(), block_name) == 0) {
617 std::vector<double> attr;
618 CHKERR it->getAttributes(attr);
619 if (attr.size() == 3) {
620 bodyForces.push_back(
621 FTensor::Tensor1<double, 3>{attr[0], attr[1], attr[2]});
622 } else {
623 SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
624 "Should be three atributes in BODYFORCE blockset, but is %d",
625 attr.size());
626 }
627 }
628 }
629
630 auto get_body_force = [this](const double, const double, const double) {
631 auto *pipeline_mng = mField.getInterface<PipelineManager>();
634 t_source(i) = 0;
635 auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
636 const auto time = fe_domain_rhs->ts_t;
637 // hardcoded gravity load
638 for (auto &t_b : bodyForces)
639 t_source(i) += (scale * t_b(i)) * time;
640 return t_source;
641 };
642
643 pipeline.push_back(new OpBodyForce("U", get_body_force));
644
645 // Calculate internal forece
646 pipeline.push_back(
647 new OpInternalForceCauchy("U", commonPlasticDataPtr->mStressPtr));
648
649 pipeline.push_back(new OpUnSetBc("U"));
651 };
652
653 auto add_domain_ops_lhs_constrain = [&](auto &pipeline, auto m_D_ptr) {
655 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
656
657 pipeline.push_back(new OpCalculatePlasticFlowLhs_dU(
658 "EP", "U", commonPlasticDataPtr, m_D_ptr));
659 pipeline.push_back(new OpCalculateContrainsLhs_dU(
660 "TAU", "U", commonPlasticDataPtr, m_D_ptr));
661
662 pipeline.push_back(new OpCalculatePlasticFlowLhs_dEP(
663 "EP", "EP", commonPlasticDataPtr, m_D_ptr));
664 pipeline.push_back(
666 pipeline.push_back(new OpCalculateContrainsLhs_dEP(
667 "TAU", "EP", commonPlasticDataPtr, m_D_ptr));
668 pipeline.push_back(
670
672 "TAU", "T", commonPlasticDataPtr));
674 "T", "EP", commonPlasticDataPtr));
675
676 pipeline.push_back(new OpUnSetBc("U"));
678 };
679
680 auto add_domain_ops_rhs_constrain = [&](auto &pipeline) {
682 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
683
684 pipeline.push_back(
686 pipeline.push_back(
688
689 pipeline.push_back(new PlasticThermalOps::OpPlasticHeatProduction(
691
693 };
694
695 auto add_boundary_ops_lhs_mechanical = [&](auto &pipeline) {
697 auto &bc_map = mField.getInterface<BcManager>()->getBcMapByBlockName();
698 for (auto bc : bc_map) {
699 if (bc_mng->checkBlock(bc, "FIX_")) {
700 pipeline.push_back(
701 new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
702 pipeline.push_back(new OpBoundaryMass(
703 "U", "U", [](double, double, double) { return 1.; },
704 bc.second->getBcEntsPtr()));
705 pipeline.push_back(new OpUnSetBc("U"));
706 }
707 }
709 };
710
711 auto add_boundary_ops_rhs_mechanical = [&](auto &pipeline) {
713
714 auto get_time = [&](double, double, double) {
715 auto *pipeline_mng = mField.getInterface<PipelineManager>();
716 auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
717 return fe_domain_rhs->ts_t;
718 };
719
720 auto get_time_scaled = [&](double, double, double) {
721 auto *pipeline_mng = mField.getInterface<PipelineManager>();
722 auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
723 return fe_domain_rhs->ts_t * scale;
724 };
725
726 auto get_minus_time = [&](double, double, double) {
727 auto *pipeline_mng = mField.getInterface<PipelineManager>();
728 auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
729 return -fe_domain_rhs->ts_t;
730 };
731
732 auto time_scaled = [&](double, double, double) {
733 auto *pipeline_mng = mField.getInterface<PipelineManager>();
734 auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
736 return amplitude_cycle * sin((2 * fe_domain_rhs->ts_t - phase_shift) *
739 } else {
740 return fe_domain_rhs->ts_t;
741 }
742 };
743
744 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
745
747 if (it->getName().compare(0, 5, "FORCE") == 0) {
748 Range force_edges;
749 std::vector<double> attr_vec;
750 CHKERR it->getMeshsetIdEntitiesByDimension(
751 mField.get_moab(), SPACE_DIM - 1, force_edges, true);
752 it->getAttributes(attr_vec);
753 if (attr_vec.size() < SPACE_DIM)
754 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
755 "Wrong size of boundary attributes vector. Set right block "
756 "size attributes.");
757 auto force_vec_ptr = boost::make_shared<MatrixDouble>(SPACE_DIM, 1);
758 std::copy(&attr_vec[0], &attr_vec[SPACE_DIM],
759 force_vec_ptr->data().begin());
760 pipeline.push_back(
761 new OpBoundaryVec("U", force_vec_ptr, get_time_scaled,
762 boost::make_shared<Range>(force_edges)));
763 }
764 }
765
766 pipeline.push_back(new OpUnSetBc("U"));
767
768 auto u_mat_ptr = boost::make_shared<MatrixDouble>();
769 pipeline.push_back(
770 new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_mat_ptr));
771
772 for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName()) {
773 if (bc_mng->checkBlock(bc, "FIX_")) {
774 pipeline.push_back(
775 new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
776 auto attr_vec = boost::make_shared<MatrixDouble>(SPACE_DIM, 1);
777 attr_vec->clear();
778 if (bc.second->bcAttributes.size() < SPACE_DIM)
779 SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
780 "Wrong size of boundary attributes vector. Set right block "
781 "size attributes. Size of attributes %d",
782 bc.second->bcAttributes.size());
783 std::copy(&bc.second->bcAttributes[0],
784 &bc.second->bcAttributes[SPACE_DIM],
785 attr_vec->data().begin());
786
787 pipeline.push_back(new OpBoundaryVec("U", attr_vec, time_scaled,
788 bc.second->getBcEntsPtr()));
789 pipeline.push_back(new OpBoundaryInternal(
790 "U", u_mat_ptr, [](double, double, double) { return 1.; },
791 bc.second->getBcEntsPtr()));
792
793 pipeline.push_back(new OpUnSetBc("U"));
794 }
795 }
796
798 };
799
800 auto add_domain_ops_lhs_thermal = [&](auto &pipeline, auto &fe) {
802 auto resistance = [](const double, const double, const double) {
803 return (1. / heat_conductivity);
804 };
805 auto capacity = [&](const double, const double, const double) {
806 return -heat_capacity * fe.ts_a;
807 };
808 auto unity = []() { return 1; };
809 pipeline.push_back(new OpHdivHdiv("FLUX", "FLUX", resistance));
810 pipeline.push_back(new OpHdivT("FLUX", "T", unity, true));
811 pipeline.push_back(new OpCapacity("T", "T", capacity));
813 };
814
815 auto add_domain_ops_rhs_thermal = [&](auto &pipeline) {
817 auto resistance = [](const double, const double, const double) {
818 return (1. / heat_conductivity);
819 };
820 auto capacity = [&](const double, const double, const double) {
821 return -heat_capacity;
822 };
823 auto unity = [](const double, const double, const double) { return 1; };
824 pipeline.push_back(new OpHdivFlux(
825 "FLUX", commonPlasticDataPtr->getTempFluxValPtr(), resistance));
826 pipeline.push_back(
827 new OpHDivTemp("FLUX", commonPlasticDataPtr->getTempValPtr(), unity));
828 pipeline.push_back(new OpBaseDivFlux(
829 "T", commonPlasticDataPtr->getTempDivFluxPtr(), unity));
830 // auto source = [&](const double x, const double y, const double z) {
831 // return 1;
832 // };
833 // pipeline.push_back(new OpHeatSource("T", source));
834
835 pipeline.push_back(new OpBaseDotT(
836 "T", commonPlasticDataPtr->getTempValDotPtr(), capacity));
838 };
839
840 auto add_boundary_ops_rhs_thermal = [&](auto &pipeline, auto &fe) {
842
843 if (SPACE_DIM == 2) {
844 pipeline.push_back(new OpSetContravariantPiolaTransformOnEdge2D());
845 } else if (SPACE_DIM == 3) {
846 pipeline.push_back(new OpHOSetCovariantPiolaTransformOnFace3D(HDIV));
847 }
848
850 const std::string block_name = "TEMPERATURE";
851 if (it->getName().compare(0, block_name.size(), block_name) == 0) {
852 Range temp_edges;
853 std::vector<double> attr_vec;
854 CHKERR it->getMeshsetIdEntitiesByDimension(
855 mField.get_moab(), SPACE_DIM - 1, temp_edges, true);
856 it->getAttributes(attr_vec);
857 if (attr_vec.size() != 1)
858 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
859 "Should be one attribute");
860 MOFEM_LOG("EXAMPLE", Sev::inform)
861 << "Set temerature " << attr_vec[0] << " on ents:\n"
862 << temp_edges;
863 bcTemperatureFunctions.push_back(BcTempFun(attr_vec[0], fe));
864 pipeline.push_back(
865 new OpTemperatureBC("FLUX", bcTemperatureFunctions.back(),
866 boost::make_shared<Range>(temp_edges)));
867 }
868 }
870 };
871
872 // Mechanics
873 CHKERR add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
874 CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainLhsPipeline(),
875 commonPlasticDataPtr->mDPtr_Deviator);
876 CHKERR add_domain_ops_lhs_mechanical(pipeline_mng->getOpDomainLhsPipeline(),
877 commonPlasticDataPtr->mDPtr_Deviator);
878 CHKERR add_domain_ops_lhs_constrain(pipeline_mng->getOpDomainLhsPipeline(),
879 commonPlasticDataPtr->mDPtr_Deviator);
880 CHKERR add_boundary_ops_lhs_mechanical(
881 pipeline_mng->getOpBoundaryLhsPipeline());
882
883 // Axiator
884 CHKERR add_domain_base_ops(feAxiatorLhs->getOpPtrVector());
885 CHKERR add_domain_ops_lhs_mechanical(feAxiatorLhs->getOpPtrVector(),
886 commonPlasticDataPtr->mDPtr_Axiator);
887 // Temperature
888 CHKERR add_domain_base_ops(feThermalLhs->getOpPtrVector());
889 CHKERR add_domain_ops_lhs_thermal(feThermalLhs->getOpPtrVector(),
890 *feThermalLhs);
891
892 // Mechanics
893 CHKERR add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
894 CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainRhsPipeline(),
895 commonPlasticDataPtr->mDPtr_Deviator);
896 CHKERR add_domain_ops_rhs_mechanical(pipeline_mng->getOpDomainRhsPipeline());
897 CHKERR add_domain_ops_rhs_constrain(pipeline_mng->getOpDomainRhsPipeline());
898 CHKERR add_boundary_ops_rhs_mechanical(
899 pipeline_mng->getOpBoundaryRhsPipeline());
900
901 // Axiator
902 CHKERR add_domain_base_ops(feAxiatorRhs->getOpPtrVector());
903 CHKERR add_domain_stress_ops(feAxiatorRhs->getOpPtrVector(),
904 commonPlasticDataPtr->mDPtr_Axiator);
905 CHKERR add_domain_ops_rhs_mechanical(feAxiatorRhs->getOpPtrVector());
906
907 // Temperature
908 CHKERR add_domain_base_ops(feThermalRhs->getOpPtrVector());
909 CHKERR add_domain_ops_rhs_thermal(feThermalRhs->getOpPtrVector());
910 CHKERR
911 add_boundary_ops_rhs_thermal(pipeline_mng->getOpBoundaryRhsPipeline(),
912 *pipeline_mng->getBoundaryRhsFE());
913
914 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule_deviator);
915 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule_deviator);
916
917 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule_bc);
918 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule_bc);
919
920 auto create_reaction_pipeline = [&](auto &pipeline) {
922
923 if (reactionMarker) {
924
925 auto det_ptr = boost::make_shared<VectorDouble>();
926 auto jac_ptr = boost::make_shared<MatrixDouble>();
927 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
928 pipeline.push_back(new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
929 pipeline.push_back(
930 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
931 pipeline.push_back(
932 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
933
934 pipeline.push_back(
936 "U", commonPlasticDataPtr->mGradPtr));
938 "EP", commonPlasticDataPtr->getPlasticStrainPtr()));
939 pipeline.push_back(
941 commonPlasticDataPtr->mStrainPtr));
942 pipeline.push_back(new PlasticThermalOps::OpPlasticStressThermal(
944 pipeline.push_back(new OpSetBc("U", false, reactionMarker));
945
946 // Calculate internal forece
947 pipeline.push_back(
948 new OpInternalForceCauchy("U", commonPlasticDataPtr->mStressPtr));
949
950 pipeline.push_back(new OpUnSetBc("U"));
951 }
952
954 };
955
956 reactionFe = boost::make_shared<DomainEle>(mField);
957 reactionFe->getRuleHook = integration_rule_deviator;
958
959 CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
960
962}
963//! [Push operators to pipeline]
964
965//! [Solve]
968
971 ISManager *is_manager = mField.getInterface<ISManager>();
972
973 auto set_section_monitor = [&](auto solver) {
975 SNES snes;
976 CHKERR TSGetSNES(solver, &snes);
977 PetscViewerAndFormat *vf;
978 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
979 PETSC_VIEWER_DEFAULT, &vf);
980 CHKERR SNESMonitorSet(
981 snes,
982 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
983 vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
985 };
986
987 auto create_post_process_element = [&]() {
989 postProcFe = boost::make_shared<PostProcEle>(mField);
990
991 auto det_ptr = boost::make_shared<VectorDouble>();
992 auto jac_ptr = boost::make_shared<MatrixDouble>();
993 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
994 postProcFe->getOpPtrVector().push_back(
995 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
996 postProcFe->getOpPtrVector().push_back(
997 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
998 postProcFe->getOpPtrVector().push_back(
999 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
1000
1001 postProcFe->getOpPtrVector().push_back(
1003 "U", commonPlasticDataPtr->mGradPtr));
1004 postProcFe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
1005 "TAU", commonPlasticDataPtr->getPlasticTauPtr()));
1006 postProcFe->getOpPtrVector().push_back(
1008 "EP", commonPlasticDataPtr->getPlasticStrainPtr()));
1009 postProcFe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
1010 "T", commonPlasticDataPtr->getTempValPtr()));
1011
1012 postProcFe->getOpPtrVector().push_back(new OpSymmetrizeTensor<SPACE_DIM>(
1013 "U", commonPlasticDataPtr->mGradPtr, commonPlasticDataPtr->mStrainPtr));
1014 postProcFe->getOpPtrVector().push_back(
1017 postProcFe->getOpPtrVector().push_back(
1019
1020 auto u_ptr = boost::make_shared<MatrixDouble>();
1021 postProcFe->getOpPtrVector().push_back(
1023
1025
1026 postProcFe->getOpPtrVector().push_back(
1027
1028 new OpPPMap(
1029
1030 postProcFe->getPostProcMesh(), postProcFe->getMapGaussPts(),
1031
1032 {{"PLASTIC_SURFACE", commonPlasticDataPtr->getPlasticSurfacePtr()},
1033 {"PLASTIC_MULTIPLIER", commonPlasticDataPtr->getPlasticTauPtr()},
1034 {"T", commonPlasticDataPtr->getTempValPtr()}},
1035
1036 {{"U", u_ptr}},
1037
1038 {},
1039
1040 {{"STRAIN", commonPlasticDataPtr->mStrainPtr},
1041 {"STRESS", commonPlasticDataPtr->mStressPtr}}
1042
1043 )
1044
1045 );
1046
1048 };
1049
1050 auto scatter_create = [&](auto D, auto coeff) {
1052 CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
1053 ROW, "U", coeff, coeff, is);
1054 int loc_size;
1055 CHKERR ISGetLocalSize(is, &loc_size);
1056 Vec v;
1057 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
1058 VecScatter scatter;
1059 CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
1060 return std::make_tuple(SmartPetscObj<Vec>(v),
1061 SmartPetscObj<VecScatter>(scatter));
1062 };
1063
1064 auto set_time_monitor = [&](auto dm, auto solver) {
1066 boost::shared_ptr<Monitor> monitor_ptr(new Monitor(
1067 dm, postProcFe, reactionFe, uXScatter, uYScatter, uZScatter));
1068 boost::shared_ptr<ForcesAndSourcesCore> null;
1069 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
1070 monitor_ptr, null, null);
1072 };
1073
1074 auto dm = simple->getDM();
1075 auto D = smartCreateDMVector(dm);
1076
1077 boost::shared_ptr<FEMethod> null;
1078 CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), feAxiatorLhs,
1079 null, null);
1080 CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(), feAxiatorRhs,
1081 null, null);
1082 CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), feThermalLhs,
1083 null, null);
1084 CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(), feThermalRhs,
1085 null, null);
1086
1087 CHKERR create_post_process_element();
1088 uXScatter = scatter_create(D, 0);
1089 uYScatter = scatter_create(D, 1);
1090 if (SPACE_DIM == 3)
1091 uZScatter = scatter_create(D, 2);
1092
1093 auto solver = pipeline_mng->createTSIM();
1094
1095 CHKERR TSSetSolution(solver, D);
1096 CHKERR set_section_monitor(solver);
1097 CHKERR set_time_monitor(dm, solver);
1098 CHKERR TSSetSolution(solver, D);
1099 CHKERR TSSetFromOptions(solver);
1100 CHKERR TSSetUp(solver);
1101 CHKERR TSSolve(solver, NULL);
1102
1103 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1104 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1105 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
1106
1108}
1109//! [Solve]
1110
1111static char help[] = "...\n\n";
1112
1113int main(int argc, char *argv[]) {
1114
1115 // Initialisation of MoFEM/PETSc and MOAB data structures
1116 const char param_file[] = "param_file.petsc";
1117 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
1118
1119 // Add logging channel for example
1120 auto core_log = logging::core::get();
1121 core_log->add_sink(
1122 LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
1123 LogManager::setLog("EXAMPLE");
1124 MOFEM_LOG_TAG("EXAMPLE", "example");
1125
1126 try {
1127
1128 //! [Register MoFEM discrete manager in PETSc]
1129 DMType dm_name = "DMMOFEM";
1130 CHKERR DMRegister_MoFEM(dm_name);
1131 //! [Register MoFEM discrete manager in PETSc
1132
1133 //! [Create MoAB]
1134 moab::Core mb_instance; ///< mesh database
1135 moab::Interface &moab = mb_instance; ///< mesh database interface
1136 //! [Create MoAB]
1137
1138 //! [Create MoFEM]
1139 MoFEM::Core core(moab); ///< finite element database
1140 MoFEM::Interface &m_field = core; ///< finite element database insterface
1141 //! [Create MoFEM]
1142
1143 //! [Load mesh]
1144 Simple *simple = m_field.getInterface<Simple>();
1145 CHKERR simple->getOptions();
1146 CHKERR simple->loadFile();
1147 //! [Load mesh]
1148
1149 //! [Example]
1150 Example ex(m_field);
1151 CHKERR ex.runProblem();
1152 //! [Example]
1153 }
1155
1157}
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:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
#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
constexpr double shear_modulus_G
Definition: elastic.cpp:46
constexpr double bulk_modulus_K
Definition: elastic.cpp:45
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
constexpr auto t_kd
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1, 1 > OpBaseTimesScalarField
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:747
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:470
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
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:800
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
#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
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
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:56
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
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:1003
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
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:1086
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
const double D
diffusivity
double A
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: plastic.cpp:70
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: plastic.cpp:68
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: plastic.cpp:56
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: plastic.cpp:77
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Body force]
Definition: plastic.cpp:54
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
Definition: plastic.cpp:79
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:75
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
BcTempFun(double v, FEMethod &fe)
double operator()(const double, const double, const double)
[Example]
Definition: plastic.cpp:111
boost::shared_ptr< DomainEle > feAxiatorRhs
Definition: plastic.cpp:130
std::vector< BcTempFun > bcTemperatureFunctions
boost::shared_ptr< PlasticThermalOps::CommonData > commonPlasticDataPtr
boost::shared_ptr< DomainEle > reactionFe
Definition: plastic.cpp:129
boost::shared_ptr< std::vector< unsigned char > > reactionMarker
Definition: plastic.cpp:138
boost::shared_ptr< PostProcEle > postProcFe
Definition: plastic.cpp:128
boost::shared_ptr< DomainEle > feAxiatorLhs
Definition: plastic.cpp:131
MoFEMErrorCode tsSolve()
FieldApproximationBase base
Definition: plot_base.cpp:68
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: plastic.cpp:134
boost::shared_ptr< DomainEle > feThermalRhs
boost::shared_ptr< DomainEle > feThermalLhs
MoFEMErrorCode createCommonData()
std::vector< FTensor::Tensor1< double, 3 > > bodyForces
Definition: plastic.cpp:140
Example(MoFEM::Interface &m_field)
MoFEMErrorCode OPs()
MoFEMErrorCode runProblem()
MoFEM::Interface & mField
Definition: plastic.cpp:118
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:135
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: plastic.cpp:137
MoFEMErrorCode setupProblem()
MoFEMErrorCode bC()
boost::shared_ptr< PlasticOps::CommonData > commonPlasticDataPtr
Definition: plastic.cpp:126
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:133
boost::shared_ptr< HenckyOps::CommonData > commonHenckyDataPtr
Definition: plastic.cpp:127
Simple interface for fast problem set-up.
Definition: BcManager.hpp:18
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Definition: BcManager.hpp:120
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.
Data on single entity (This is passed as argument to DataOperator::doWork)
structure for User Loop Methods on finite elements
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
Get vector field for H-div approximation.
Calculate divergence of vector field.
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.
transform Hcurl base fluxes from reference element to physical triangle
Make Hdiv space from Hcurl space in 2d.
Post post-proc data at points from hash maps.
Definition: PostProc.hpp:166
Scale base functions by inverses of measure of element.
Set indices on entities on finite element.
Set inverse jacobian to base functions.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
Definition: PostProc.hpp:120
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCore VolEle
double young_modulus
double heat_conductivity
int main(int argc, char *argv[])
double omega_inf
double amplitude_cycle
EntitiesFieldData::EntData EntData
double Qinf
static char help[]
[Solve]
double amplitude_shift
double phase_shift
double rho
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
#define EXECUTABLE_DIMENSION
constexpr int SPACE_DIM
double visH
double poisson_ratio
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< SPACE_DIM > OpHdivT
Integrate Lhs div of base of flux time base of temperature (FLUX x T) and transpose of it,...
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 3, 3, 1 > OpHdivFlux
Integrating Rhs flux base (1/k) flux (FLUX)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpBaseDotT
Integrate Rhs base of temerature time heat capacity times heat rate (T)
int number_of_cycles_in_total_time
double fraction_of_dissipation
double scale
long double hardening(long double tau, double temp)
double H
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpHeatSource
double omega_0
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpCapacity
Integrate Lhs base of temerature times (heat capacity) times base of temperature (T x T)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 1, 2 > OpHDivTemp
Integrate Rhs div flux base times temperature (T)
double heat_capacity
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, 3 > OpHdivHdiv
Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Body force]
long double hardening_dtau(long double tau, double temp)
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
int order
double ref_temp
double b_iso
double coeff_expansion
double sigmaY
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpNormalMixVecTimesScalar< SPACE_DIM > OpTemperatureBC
long double hardening_dtemp(long double tau, double temp)
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
double cn
OpBaseDotT OpBaseDivFlux
Integrate Rhs base of temerature times divergenc of flux (T)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
[Body force]
constexpr EntityType boundary_ent
double omega_h