v0.13.2
Loading...
Searching...
No Matches
thermo_elastic.cpp
Go to the documentation of this file.
1/**
2 * \file thermo_elastic.cpp
3 * \example thermo_elastic.cpp
4 *
5 * Thermo plasticity
6 *
7 */
8
9#ifndef EXECUTABLE_DIMENSION
10#define EXECUTABLE_DIMENSION 3
11#endif
12
13#include <MoFEM.hpp>
14
15using namespace MoFEM;
16
17template <int DIM> struct ElementsAndOps {};
18
19template <> struct ElementsAndOps<2> {
22};
23
24template <> struct ElementsAndOps<3> {
27};
28
29constexpr int SPACE_DIM =
30 EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
31
38
41
42//! [Linear elastic problem]
44 GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM,
45 0>; //< Elastic stiffness matrix
48 GAUSS>::OpGradTimesSymTensor<1, SPACE_DIM,
49 SPACE_DIM>; //< Elastic internal forces
50//! [Linear elastic problem]
51
52//! [Thermal problem]
53/**
54 * @brief Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)
55 *
56 */
59
60/**
61 * @brief Integrate Lhs div of base of flux time base of temperature (FLUX x T)
62 * and transpose of it, i.e. (T x FLAX)
63 *
64 */
66 GAUSS>::OpMixDivTimesScalar<SPACE_DIM>;
67
68/**
69 * @brief Integrate Lhs base of temperature times (heat capacity) times base of
70 * temperature (T x T)
71 *
72 */
75
76/**
77 * @brief Integrating Rhs flux base (1/k) flux (FLUX)
78 */
80 GAUSS>::OpBaseTimesVector<3, SPACE_DIM, 1>;
81
82/**
83 * @brief Integrate Rhs div flux base times temperature (T)
84 *
85 */
87 GAUSS>::OpMixDivTimesU<3, 1, 2>;
88
89/**
90 * @brief Integrate Rhs base of temperature time heat capacity times heat rate
91 * (T)
92 *
93 */
95 GAUSS>::OpBaseTimesScalar<1>;
96
97/**
98 * @brief Integrate Rhs base of temperature times divergence of flux (T)
99 *
100 */
102
103//! [Thermal problem]
104
105//! [Body and heat source]
112//! [Body and heat source]
113
114//! [Natural boundary conditions]
120//! [Natural boundary conditions]
121
122//! [Essential boundary conditions (Least square approach)]
129//! [Essential boundary conditions (Least square approach)]
130
131double young_modulus = 1;
132double poisson_ratio = 0.25;
134double ref_temp = 0.0;
135
137 1; // Force / (time temperature ) or Power /
138 // (length temperature) // Time unit is hour. force unit kN
139double heat_capacity = 1; // length^2/(time^2 temperature) // length is
140 // millimeter time is hour
141
142int order = 2; //< default approximation order
143
144#include <ThermoElasticOps.hpp> //< additional coupling opearyors
145using namespace ThermoElasticOps; //< name space of coupling operators
146
148
150
152
153private:
155
156 MoFEMErrorCode setupProblem(); ///< add fields
157 MoFEMErrorCode createCommonData(); //< read global data from command line
158 MoFEMErrorCode bC(); //< read boundary conditions
159 MoFEMErrorCode OPs(); //< add operators to pipeline
160 MoFEMErrorCode tsSolve(); //< time solver
161
162 boost::shared_ptr<MatrixDouble> getMatDPtr();
163};
164
165//! [Run problem]
170 CHKERR bC();
171 CHKERR OPs();
172 CHKERR tsSolve();
174}
175//! [Run problem]
176
177//! [Set up problem]
181 // Add field
183 // Mechanical fields
184 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
185 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
186 // Temperature
187 const auto flux_space = (SPACE_DIM == 2) ? HCURL : HDIV;
188 CHKERR simple->addDomainField("T", L2, AINSWORTH_LEGENDRE_BASE, 1);
189 CHKERR simple->addDomainField("FLUX", flux_space, DEMKOWICZ_JACOBI_BASE, 1);
190 CHKERR simple->addBoundaryField("FLUX", flux_space, DEMKOWICZ_JACOBI_BASE, 1);
191
192 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
193 CHKERR simple->setFieldOrder("U", order);
194 CHKERR simple->setFieldOrder("FLUX", order);
195 CHKERR simple->setFieldOrder("T", order - 1);
196 CHKERR simple->setUp();
198}
199//! [Set up problem]
200
201//! [Create common data]
204
205 auto get_command_line_parameters = [&]() {
207 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
208 &young_modulus, PETSC_NULL);
209 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
210 &poisson_ratio, PETSC_NULL);
211 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-coeff_expansion",
212 &coeff_expansion, PETSC_NULL);
213 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-ref_temp", &ref_temp,
214 PETSC_NULL);
215
216 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-capacity", &heat_capacity,
217 PETSC_NULL);
218 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-conductivity",
219 &heat_conductivity, PETSC_NULL);
220
221 MOFEM_LOG("ThermoElastic", Sev::inform)
222 << "Young modulus " << young_modulus;
223 MOFEM_LOG("ThermoElastic", Sev::inform)
224 << "Poisson ratio " << poisson_ratio;
225 MOFEM_LOG("ThermoElastic", Sev::inform)
226 << "Coeff_expansion " << coeff_expansion;
227 MOFEM_LOG("ThermoElastic", Sev::inform)
228 << "Reference_temperature " << ref_temp;
229
231 };
232
233 CHKERR get_command_line_parameters();
234
236}
237//! [Create common data]
238
239//! [Boundary condition]
242
244 auto bc_mng = mField.getInterface<BcManager>();
245
246 CHKERR bc_mng->removeBlockDOFsOnEntities<DisplacementCubitBcData>(
247 simple->getProblemName(), "U");
248 CHKERR bc_mng->pushMarkDOFsOnEntities<HeatFluxCubitBcData>(
249 simple->getProblemName(), "FLUX");
251}
252//! [Boundary condition]
253
254//! [Push operators to pipeline]
257 auto pipeline_mng = mField.getInterface<PipelineManager>();
259 auto bc_mng = mField.getInterface<BcManager>();
260
261 auto boundary_marker =
262 bc_mng->getMergedBlocksMarker(vector<string>{"HEATFLUX"});
263 auto mDPtr = getMatDPtr();
264
265 auto integration_rule = [](int, int, int approx_order) {
266 return 2 * approx_order;
267 };
268 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
269 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
270 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
271 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
272
273 auto time_scale = boost::make_shared<TimeScale>();
274
275 auto add_domain_ops = [&](auto &pipeline) {
277
278 auto det_ptr = boost::make_shared<VectorDouble>();
279 auto jac_ptr = boost::make_shared<MatrixDouble>();
280 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
281
282 pipeline.push_back(new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
283 pipeline.push_back(
284 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
285 pipeline.push_back(
286 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
287
288 if constexpr (SPACE_DIM == 2) {
289 pipeline.push_back(new OpMakeHdivFromHcurl());
290 pipeline.push_back(new OpSetContravariantPiolaTransformOnFace2D(jac_ptr));
291 pipeline.push_back(new OpSetInvJacHcurlFace(inv_jac_ptr));
292 pipeline.push_back(new OpSetInvJacL2ForFace(inv_jac_ptr));
293 } else {
294 pipeline.push_back(
295 new OpSetHOContravariantPiolaTransform(HDIV, det_ptr, jac_ptr));
296 pipeline.push_back(new OpSetHOInvJacVectorBase(HDIV, inv_jac_ptr));
297 pipeline.push_back(new OpSetInvJacL2ForFace(inv_jac_ptr));
298 }
299
301 };
302
303 auto add_domain_rhs_ops = [&](auto &pipeline) {
305
306 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
307 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
308 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
309
310 auto vec_temp_ptr = boost::make_shared<VectorDouble>();
311 auto vec_temp_dot_ptr = boost::make_shared<VectorDouble>();
312 auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
313 auto vec_temp_div_ptr = boost::make_shared<VectorDouble>();
314
315 pipeline.push_back(new OpCalculateScalarFieldValues("T", vec_temp_ptr));
316 pipeline.push_back(
317 new OpCalculateScalarFieldValuesDot("T", vec_temp_dot_ptr));
319 "FLUX", vec_temp_div_ptr));
320 pipeline.push_back(
321 new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
322
324 "U", mat_grad_ptr));
325 pipeline.push_back(
326 new OpSymmetrizeTensor<SPACE_DIM>("U", mat_grad_ptr, mat_strain_ptr));
327 pipeline.push_back(new OpStressThermal("U", mat_strain_ptr, vec_temp_ptr,
328 mDPtr, mat_stress_ptr));
329
330 pipeline.push_back(new OpSetBc("FLUX", true, boundary_marker));
331
332 pipeline.push_back(new OpInternalForceCauchy(
333 "U", mat_stress_ptr,
334 [](double, double, double) constexpr { return 1; }));
335
336 auto resistance = [](const double, const double, const double) {
337 return (1. / heat_conductivity);
338 };
339 auto capacity = [&](const double, const double, const double) {
340 return heat_capacity;
341 };
342 auto unity = [](const double, const double, const double) { return -1.; };
343 pipeline.push_back(new OpHdivFlux("FLUX", mat_flux_ptr, resistance));
344 pipeline.push_back(new OpHDivTemp("FLUX", vec_temp_ptr, unity));
345 pipeline.push_back(new OpBaseDivFlux("T", vec_temp_div_ptr, unity));
346 pipeline.push_back(new OpBaseDotT("T", vec_temp_dot_ptr, capacity));
347
349 pipeline, mField, "T", {time_scale}, "HEAT_SOURCE", Sev::inform);
351 pipeline_mng->getOpDomainRhsPipeline(), mField, "U", {time_scale},
352 "BODY_FORCE", Sev::inform);
353
354 pipeline.push_back(new OpUnSetBc("FLUX"));
356 };
357
358 auto add_domain_lhs_ops = [&](auto &pipeline) {
360 pipeline.push_back(new OpSetBc("FLUX", true, boundary_marker));
361
362 pipeline.push_back(new OpKCauchy("U", "U", mDPtr));
363 pipeline.push_back(
365
366 auto resistance = [](const double, const double, const double) {
367 return (1. / heat_conductivity);
368 };
369 auto capacity = [](const double, const double, const double) {
370 return heat_capacity;
371 };
372 pipeline.push_back(new OpHdivHdiv("FLUX", "FLUX", resistance));
373 pipeline.push_back(new OpHdivT(
374 "FLUX", "T", []() { return -1; }, true));
375
376 auto op_capacity = new OpCapacity("T", "T", capacity);
377 op_capacity->feScalingFun = [](const FEMethod *fe_ptr) {
378 return fe_ptr->ts_a;
379 };
380 pipeline.push_back(op_capacity);
381
382 pipeline.push_back(new OpUnSetBc("FLUX"));
384 };
385
386 auto add_boundary_rhs_ops = [&](auto &pipeline) {
388
389 if constexpr (SPACE_DIM == 2) {
390 pipeline.push_back(new OpSetContravariantPiolaTransformOnEdge2D());
391 } else {
392 pipeline.push_back(new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
393 }
394
395 pipeline.push_back(new OpSetBc("FLUX", true, boundary_marker));
396
398 pipeline_mng->getOpBoundaryRhsPipeline(), mField, "U", {time_scale},
399 "FORCE", Sev::inform);
400
402 pipeline_mng->getOpBoundaryRhsPipeline(), mField, "FLUX", {time_scale},
403 "TEMPERATURE", Sev::inform);
404
405 pipeline.push_back(new OpUnSetBc("FLUX"));
406
407 auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
408 pipeline.push_back(
409 new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
412 mField, pipeline, simple->getProblemName(), "FLUX", mat_flux_ptr,
413 {time_scale});
414
416 };
417
418 auto add_boundary_lhs_ops = [&](auto &pipeline) {
420
421 if constexpr (SPACE_DIM == 2) {
422 pipeline.push_back(new OpSetContravariantPiolaTransformOnEdge2D());
423 } else {
424 pipeline.push_back(new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
425 }
426
429 mField, pipeline, simple->getProblemName(), "FLUX");
430
432 };
433
434 auto get_bc_hook_rhs = [&]() {
436 mField, pipeline_mng->getDomainRhsFE(), {time_scale});
437 return hook;
438 };
439 auto get_bc_hook_lhs = [&]() {
441 mField, pipeline_mng->getDomainLhsFE(), {time_scale});
442 return hook;
443 };
444
445 pipeline_mng->getDomainRhsFE()->preProcessHook = get_bc_hook_rhs();
446 pipeline_mng->getDomainLhsFE()->preProcessHook = get_bc_hook_lhs();
447
448 CHKERR add_domain_ops(pipeline_mng->getOpDomainRhsPipeline());
449 CHKERR add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
450
451 CHKERR add_domain_ops(pipeline_mng->getOpDomainLhsPipeline());
452 CHKERR add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
453
454 CHKERR add_boundary_rhs_ops(pipeline_mng->getOpBoundaryRhsPipeline());
455 CHKERR add_boundary_lhs_ops(pipeline_mng->getOpBoundaryLhsPipeline());
456
458}
459//! [Push operators to pipeline]
460
461//! [Solve]
464
467 ISManager *is_manager = mField.getInterface<ISManager>();
468
469 auto dm = simple->getDM();
470 auto solver = pipeline_mng->createTSIM();
471 auto snes_ctx_ptr = getDMSnesCtx(dm);
472
473 auto set_section_monitor = [&](auto solver) {
475 SNES snes;
476 CHKERR TSGetSNES(solver, &snes);
477 CHKERR SNESMonitorSet(snes,
478 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
480 (void *)(snes_ctx_ptr.get()), nullptr);
482 };
483
484 auto create_post_process_element = [&]() {
485 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
486
487 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
488 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
489 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
490
491 auto vec_temp_ptr = boost::make_shared<VectorDouble>();
492 auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
493
494 auto det_ptr = boost::make_shared<VectorDouble>();
495 auto jac_ptr = boost::make_shared<MatrixDouble>();
496 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
497
498 post_proc_fe->getOpPtrVector().push_back(
499 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
500 post_proc_fe->getOpPtrVector().push_back(
501 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
502 post_proc_fe->getOpPtrVector().push_back(
503 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
504
505 if constexpr (SPACE_DIM == 2) {
506 post_proc_fe->getOpPtrVector().push_back(new OpMakeHdivFromHcurl());
507 post_proc_fe->getOpPtrVector().push_back(
509 post_proc_fe->getOpPtrVector().push_back(
510 new OpSetInvJacHcurlFace(inv_jac_ptr));
511 post_proc_fe->getOpPtrVector().push_back(
512 new OpSetInvJacL2ForFace(inv_jac_ptr));
513 } else {
514 post_proc_fe->getOpPtrVector().push_back(
515 new OpSetHOContravariantPiolaTransform(HDIV, det_ptr, jac_ptr));
516 post_proc_fe->getOpPtrVector().push_back(
517 new OpSetHOInvJacVectorBase(HDIV, inv_jac_ptr));
518 }
519
520 post_proc_fe->getOpPtrVector().push_back(
521 new OpCalculateScalarFieldValues("T", vec_temp_ptr));
522 post_proc_fe->getOpPtrVector().push_back(
523 new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
524
525 auto u_ptr = boost::make_shared<MatrixDouble>();
526 post_proc_fe->getOpPtrVector().push_back(
528 post_proc_fe->getOpPtrVector().push_back(
530 mat_grad_ptr));
531 post_proc_fe->getOpPtrVector().push_back(
532 new OpSymmetrizeTensor<SPACE_DIM>("U", mat_grad_ptr, mat_strain_ptr));
533 post_proc_fe->getOpPtrVector().push_back(new OpStressThermal(
534 "U", mat_strain_ptr, vec_temp_ptr, getMatDPtr(), mat_stress_ptr));
535
537
538 post_proc_fe->getOpPtrVector().push_back(
539
540 new OpPPMap(
541
542 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
543
544 {{"T", vec_temp_ptr}},
545
546 {{"U", u_ptr}, {"FLUX", mat_flux_ptr}},
547
548 {},
549
550 {{"STRAIN", mat_strain_ptr}, {"STRESS", mat_stress_ptr}}
551
552 )
553
554 );
555
556 return post_proc_fe;
557 };
558
559 auto monitor_ptr = boost::make_shared<FEMethod>();
560 auto post_proc_fe = create_post_process_element();
561
562 auto set_time_monitor = [&](auto dm, auto solver) {
564 monitor_ptr->preProcessHook = [&]() {
566 CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
567 CHKERR post_proc_fe->writeFile(
568 "out_" + boost::lexical_cast<std::string>(monitor_ptr->ts_step) +
569 ".h5m");
571 };
572 auto null = boost::shared_ptr<FEMethod>();
573 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(), null,
574 monitor_ptr, null);
576 };
577
578 auto set_fieldsplit_preconditioner = [&](auto solver) {
580
581 SNES snes;
582 CHKERR TSGetSNES(solver, &snes);
583 KSP ksp;
584 CHKERR SNESGetKSP(snes, &ksp);
585 PC pc;
586 CHKERR KSPGetPC(ksp, &pc);
587 PetscBool is_pcfs = PETSC_FALSE;
588 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
589
590 // Setup fieldsplit (block) solver - optional: yes/no
591 if (is_pcfs == PETSC_TRUE) {
592 auto bc_mng = mField.getInterface<BcManager>();
593 auto name_prb = simple->getProblemName();
594 auto is_all_bc = bc_mng->getBlockIS(name_prb, "HEATFLUX", "FLUX", 0, 0);
595 int is_all_bc_size;
596 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
597 MOFEM_LOG("ThermoElastic", Sev::inform)
598 << "Field split block size " << is_all_bc_size;
599 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
600 is_all_bc); // boundary block
601 }
602
604 };
605
606 auto D = createDMVector(dm);
607 CHKERR TSSetSolution(solver, D);
608 CHKERR TSSetFromOptions(solver);
609 CHKERR set_section_monitor(solver);
610 CHKERR set_time_monitor(dm, solver);
611 CHKERR set_fieldsplit_preconditioner(solver);
612 CHKERR TSSetUp(solver);
613 CHKERR TSSolve(solver, NULL);
614
616}
617//! [Solve]
618
619boost::shared_ptr<MatrixDouble> ThermoElasticProblem::getMatDPtr() {
620 auto set_matrial_stiffness = [&](auto mDPtr) {
626 const double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
627 const double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
628
629 // Plane stress or when 1, plane strain or 3d
630 double A = (SPACE_DIM == 2)
631 ? 2 * shear_modulus_G /
632 (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
633 : 1;
634 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mDPtr);
635 t_D(i, j, k, l) = 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
636 A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
637 t_kd(i, j) * t_kd(k, l);
638 return mDPtr;
639 };
640
641 constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
642 return set_matrial_stiffness(
643 boost::make_shared<MatrixDouble>(size_symm * size_symm, 1));
644}
645
646static char help[] = "...\n\n";
647
648int main(int argc, char *argv[]) {
649
650 // Initialisation of MoFEM/PETSc and MOAB data structures
651 const char param_file[] = "param_file.petsc";
653
654 // Add logging channel for example
655 auto core_log = logging::core::get();
656 core_log->add_sink(
658 LogManager::setLog("ThermoElastic");
659 MOFEM_LOG_TAG("ThermoElastic", "ThermoElastic");
660
661 try {
662
663 //! [Register MoFEM discrete manager in PETSc]
664 DMType dm_name = "DMMOFEM";
665 CHKERR DMRegister_MoFEM(dm_name);
666 //! [Register MoFEM discrete manager in PETSc
667
668 //! [Create MoAB]
669 moab::Core mb_instance; ///< mesh database
670 moab::Interface &moab = mb_instance; ///< mesh database interface
671 //! [Create MoAB]
672
673 //! [Create MoFEM]
674 MoFEM::Core core(moab); ///< finite element database
675 MoFEM::Interface &m_field = core; ///< finite element database interface
676 //! [Create MoFEM]
677
678 //! [Load mesh]
679 Simple *simple = m_field.getInterface<Simple>();
680 CHKERR simple->getOptions();
681 CHKERR simple->loadFile();
682 //! [Load mesh]
683
684 //! [ThermoElasticProblem]
685 ThermoElasticProblem ex(m_field);
686 CHKERR ex.runProblem();
687 //! [ThermoElasticProblem]
688 }
690
692}
std::string param_file
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
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
#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
#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
@ HCURL
field with continuous tangents
Definition: definitions.h:86
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#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:79
constexpr double bulk_modulus_K
Definition: elastic.cpp:78
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
auto integration_rule
constexpr auto t_kd
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:572
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
FTensor::Index< 'i', SPACE_DIM > i
double D
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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: DMMoFEM.cpp:1042
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:226
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:1031
constexpr AssemblyType A
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< G >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: plastic.cpp:54
constexpr auto size_symm
Definition: plastic.cpp:33
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
Simple interface for fast problem set-up.
Definition: BcManager.hpp:24
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
Class (Function) to enforce essential constrains.
Definition: Essential.hpp:39
structure for User Loop Methods on finite elements
Definition of the heat flux bc data structure.
Definition: BCData.hpp:422
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
Assembly methods.
Definition: Natural.hpp:65
Get vector field for H-div approximation.
Calculate divergence of vector field.
Get value at integration points for scalar field.
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
transform Hdiv 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.
Set indices on entities on finite element.
Apply contravariant (Piola) transfer to Hdiv space for HO geometry.
Set inverse jacobian to base functions.
transform local reference derivatives of shape function to global derivatives if higher order geometr...
PipelineManager interface.
MoFEM::VolumeElementForcesAndSourcesCore VolEle
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode runProblem()
[Run problem]
MoFEM::Interface & mField
MoFEMErrorCode createCommonData()
[Set up problem]
ThermoElasticProblem(MoFEM::Interface &m_field)
boost::shared_ptr< MatrixDouble > getMatDPtr()
[Solve]
MoFEMErrorCode setupProblem()
add fields
MoFEMErrorCode bC()
[Create common data]
MoFEMErrorCode OPs()
[Boundary condition]
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
double young_modulus
[Essential boundary conditions (Least square approach)]
double heat_conductivity
static char help[]
#define EXECUTABLE_DIMENSION
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBaseDotT
Integrate Rhs base of temperature time heat capacity times heat rate (T)
constexpr int SPACE_DIM
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, SPACE_DIM, 1 > OpHdivFlux
Integrating Rhs flux base (1/k) flux (FLUX)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpCapacity
Integrate Lhs base of temperature times (heat capacity) times base of temperature (T x T)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, SPACE_DIM > OpHdivHdiv
[Linear elastic problem]
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 >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Linear elastic problem]
int order
double ref_temp
double coeff_expansion
OpBaseDotT OpBaseDivFlux
Integrate Rhs base of temperature times divergence of flux (T)