v0.14.0
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
17constexpr int SPACE_DIM =
18 EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
19
27
30
31//! [Linear elastic problem]
33 GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM,
34 0>; //< Elastic stiffness matrix
37 GAUSS>::OpGradTimesSymTensor<1, SPACE_DIM,
38 SPACE_DIM>; //< Elastic internal forces
39//! [Linear elastic problem]
40
41//! [Thermal problem]
42/**
43 * @brief Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)
44 *
45 */
48
49/**
50 * @brief Integrate Lhs div of base of flux time base of temperature (FLUX x T)
51 * and transpose of it, i.e. (T x FLAX)
52 *
53 */
55 GAUSS>::OpMixDivTimesScalar<SPACE_DIM>;
56
57/**
58 * @brief Integrate Lhs base of temperature times (heat capacity) times base of
59 * temperature (T x T)
60 *
61 */
64
65/**
66 * @brief Integrating Rhs flux base (1/k) flux (FLUX)
67 */
69 GAUSS>::OpBaseTimesVector<3, SPACE_DIM, 1>;
70
71/**
72 * @brief Integrate Rhs div flux base times temperature (T)
73 *
74 */
76 GAUSS>::OpMixDivTimesU<3, 1, SPACE_DIM>;
77
78/**
79 * @brief Integrate Rhs base of temperature time heat capacity times heat rate
80 * (T)
81 *
82 */
84 GAUSS>::OpBaseTimesScalar<1>;
85
86/**
87 * @brief Integrate Rhs base of temperature times divergence of flux (T)
88 *
89 */
91
92//! [Thermal problem]
93
94//! [Body and heat source]
103//! [Body and heat source]
104
105//! [Natural boundary conditions]
111//! [Natural boundary conditions]
112
113//! [Essential boundary conditions (Least square approach)]
116 GAUSS>::OpEssentialRhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
119 GAUSS>::OpEssentialLhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
120//! [Essential boundary conditions (Least square approach)]
121
124double ref_temp = 0.0;
125
128 1; // Force / (time temperature ) or Power /
129 // (length temperature) // Time unit is hour. force unit kN
130double default_heat_capacity = 1; // length^2/(time^2 temperature) // length is
131 // millimeter time is hour
132
133int order = 2; //< default approximation order
134
135#include <ThermoElasticOps.hpp> //< additional coupling opearyors
136using namespace ThermoElasticOps; //< name space of coupling operators
137
142
143auto save_range = [](moab::Interface &moab, const std::string name,
144 const Range r) {
146 auto out_meshset = get_temp_meshset_ptr(moab);
147 CHKERR moab.add_entities(*out_meshset, r);
148 CHKERR moab.write_file(name.c_str(), "VTK", "", out_meshset->get_ptr(), 1);
150};
151
153
155
157
158private:
160
161 PetscBool doEvalField;
162 std::array<double, SPACE_DIM> fieldEvalCoords;
163 boost::shared_ptr<FieldEvaluatorInterface::SetPtsData> fieldEvalData;
164
165 boost::shared_ptr<VectorDouble> scalarFieldPtr;
166 boost::shared_ptr<MatrixDouble> vectorFieldPtr;
167 boost::shared_ptr<MatrixDouble> tensorFieldPtr;
168
169 MoFEMErrorCode setupProblem(); ///< add fields
170 MoFEMErrorCode createCommonData(); //< read global data from command line
171 MoFEMErrorCode bC(); //< read boundary conditions
172 MoFEMErrorCode OPs(); //< add operators to pipeline
173 MoFEMErrorCode tsSolve(); //< time solver
174
176 : public boost::enable_shared_from_this<BlockedParameters> {
181
182 inline auto getDPtr() {
183 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &mD);
184 }
185
186 inline auto getCoeffExpansionPtr() {
187 return boost::shared_ptr<double>(shared_from_this(), &coeffExpansion);
188 }
189
191 return boost::shared_ptr<double>(shared_from_this(), &heatConductivity);
192 }
193
194 inline auto getHeatCapacityPtr() {
195 return boost::shared_ptr<double>(shared_from_this(), &heatCapacity);
196 }
197 };
198
200 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
201 std::string block_elastic_name, std::string block_thermal_name,
202 boost::shared_ptr<BlockedParameters> blockedParamsPtr, Sev sev);
203};
204
206 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
207 std::string block_elastic_name, std::string block_thermal_name,
208 boost::shared_ptr<BlockedParameters> blockedParamsPtr, Sev sev) {
210
211 struct OpMatElasticBlocks : public DomainEleOp {
212 OpMatElasticBlocks(boost::shared_ptr<MatrixDouble> m, double bulk_modulus_K,
213 double shear_modulus_G, MoFEM::Interface &m_field,
214 Sev sev,
215 std::vector<const CubitMeshSets *> meshset_vec_ptr)
217 bulkModulusKDefault(bulk_modulus_K),
218 shearModulusGDefault(shear_modulus_G) {
219 CHK_THROW_MESSAGE(extractElasticBlockData(m_field, meshset_vec_ptr, sev),
220 "Can not get data from block");
221 }
222
223 MoFEMErrorCode doWork(int side, EntityType type,
226
227 for (auto &b : blockData) {
228
229 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
230 CHKERR getMatDPtr(matDPtr, b.bulkModulusK, b.shearModulusG);
232 }
233 }
234
235 CHKERR getMatDPtr(matDPtr, bulkModulusKDefault, shearModulusGDefault);
237 }
238
239 private:
240 boost::shared_ptr<MatrixDouble> matDPtr;
241
242 struct BlockData {
243 double bulkModulusK;
244 double shearModulusG;
245 Range blockEnts;
246 };
247
248 double bulkModulusKDefault;
249 double shearModulusGDefault;
250 std::vector<BlockData> blockData;
251
253 extractElasticBlockData(MoFEM::Interface &m_field,
254 std::vector<const CubitMeshSets *> meshset_vec_ptr,
255 Sev sev) {
257
258 for (auto m : meshset_vec_ptr) {
259 MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Elastic Block") << *m;
260 std::vector<double> block_data;
261 CHKERR m->getAttributes(block_data);
262 if (block_data.size() < 2) {
263 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
264 "Expected that block has two attributes");
265 }
266 auto get_block_ents = [&]() {
267 Range ents;
268 CHKERR
269 m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
270 return ents;
271 };
272
273 double young_modulus = block_data[0];
274 double poisson_ratio = block_data[1];
275 double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
276 double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
277
278 MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Elastic Block")
279 << m->getName() << ": E = " << young_modulus
280 << " nu = " << poisson_ratio;
281
282 blockData.push_back(
283 {bulk_modulus_K, shear_modulus_G, get_block_ents()});
284 }
285 MOFEM_LOG_CHANNEL("WORLD");
287 }
288
289 MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
290 double bulk_modulus_K, double shear_modulus_G) {
292 //! [Calculate elasticity tensor]
293 auto set_material_stiffness = [&]() {
299 double A = (SPACE_DIM == 2)
300 ? 2 * shear_modulus_G /
301 (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
302 : 1;
303 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
304 t_D(i, j, k, l) =
305 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
306 A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
307 t_kd(k, l);
308 };
309 //! [Calculate elasticity tensor]
310 constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
311 mat_D_ptr->resize(size_symm * size_symm, 1);
312 set_material_stiffness();
314 }
315 };
316
317 double default_bulk_modulus_K =
319 double default_shear_modulus_G =
321
322 pipeline.push_back(new OpMatElasticBlocks(
323 blockedParamsPtr->getDPtr(), default_bulk_modulus_K,
324 default_bulk_modulus_K, mField, sev,
325
326 // Get blockset using regular expression
327 mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
328
329 (boost::format("%s(.*)") % block_elastic_name).str()
330
331 ))
332
333 ));
334
335 struct OpMatThermalBlocks : public DomainEleOp {
336 OpMatThermalBlocks(boost::shared_ptr<double> expansion_ptr,
337 boost::shared_ptr<double> conductivity_ptr,
338 boost::shared_ptr<double> capacity_ptr,
339 MoFEM::Interface &m_field, Sev sev,
340 std::vector<const CubitMeshSets *> meshset_vec_ptr)
342 expansionPtr(expansion_ptr), conductivityPtr(conductivity_ptr),
343 capacityPtr(capacity_ptr) {
344 CHK_THROW_MESSAGE(extractThermallockData(m_field, meshset_vec_ptr, sev),
345 "Can not get data from block");
346 }
347
348 MoFEMErrorCode doWork(int side, EntityType type,
351
352 for (auto &b : blockData) {
353
354 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
355 *expansionPtr = b.expansion;
356 *conductivityPtr = b.conductivity;
357 *capacityPtr = b.capacity;
359 }
360 }
361
362 *expansionPtr = default_coeff_expansion;
363 *conductivityPtr = default_heat_conductivity;
364 *capacityPtr = default_heat_capacity;
365
367 }
368
369 private:
370 struct BlockData {
371 double expansion;
372 double conductivity;
373 double capacity;
374 Range blockEnts;
375 };
376
377 std::vector<BlockData> blockData;
378
380 extractThermallockData(MoFEM::Interface &m_field,
381 std::vector<const CubitMeshSets *> meshset_vec_ptr,
382 Sev sev) {
384
385 for (auto m : meshset_vec_ptr) {
386 MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Thermal Block") << *m;
387 std::vector<double> block_data;
388 CHKERR m->getAttributes(block_data);
389 if (block_data.size() < 3) {
390 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
391 "Expected that block has two attributes");
392 }
393 auto get_block_ents = [&]() {
394 Range ents;
395 CHKERR
396 m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
397 return ents;
398 };
399
400 MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Thermal Block")
401 << m->getName() << ": expansion = " << block_data[0]
402 << " conductivity = " << block_data[1] << " capacity "
403 << block_data[2];
404
405 blockData.push_back(
406 {block_data[0], block_data[1], block_data[2], get_block_ents()});
407 }
408 MOFEM_LOG_CHANNEL("WORLD");
410 }
411
412 boost::shared_ptr<double> expansionPtr;
413 boost::shared_ptr<double> conductivityPtr;
414 boost::shared_ptr<double> capacityPtr;
415 };
416
417 pipeline.push_back(new OpMatThermalBlocks(
418 blockedParamsPtr->getCoeffExpansionPtr(),
419 blockedParamsPtr->getHeatConductivityPtr(),
420 blockedParamsPtr->getHeatCapacityPtr(), mField, sev,
421
422 // Get blockset using regular expression
424
425 (boost::format("%s(.*)") % block_thermal_name).str()
426
427 ))
428
429 ));
430
432}
433
434//! [Run problem]
439 CHKERR bC();
440 CHKERR OPs();
441 CHKERR tsSolve();
443}
444//! [Run problem]
445
446//! [Set up problem]
450 // Add field
452 // Mechanical fields
453 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
454 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
455 // Temperature
456 constexpr auto flux_space = (SPACE_DIM == 2) ? HCURL : HDIV;
457 CHKERR simple->addDomainField("T", L2, AINSWORTH_LEGENDRE_BASE, 1);
458 CHKERR simple->addDomainField("FLUX", flux_space, DEMKOWICZ_JACOBI_BASE, 1);
459 CHKERR simple->addBoundaryField("FLUX", flux_space, DEMKOWICZ_JACOBI_BASE, 1);
460
461 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
462 CHKERR simple->setFieldOrder("U", order + 1);
463 CHKERR simple->setFieldOrder("FLUX", order + 1);
464 CHKERR simple->setFieldOrder("T", order);
465 CHKERR simple->setUp();
466
467 int coords_dim = SPACE_DIM;
468 CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
469 fieldEvalCoords.data(), &coords_dim,
470 &doEvalField);
471
472 scalarFieldPtr = boost::make_shared<VectorDouble>();
473 vectorFieldPtr = boost::make_shared<MatrixDouble>();
474 tensorFieldPtr = boost::make_shared<MatrixDouble>();
475
476 if (doEvalField) {
478 mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
479
480 if constexpr (SPACE_DIM == 3) {
482 fieldEvalData, simple->getDomainFEName());
483 } else {
485 fieldEvalData, simple->getDomainFEName());
486 }
487
488 fieldEvalData->setEvalPoints(fieldEvalCoords.data(), 1);
489 auto no_rule = [](int, int, int) { return -1; };
490
491 auto field_eval_fe_ptr = fieldEvalData->feMethodPtr.lock();
492 field_eval_fe_ptr->getRuleHook = no_rule;
493
494 field_eval_fe_ptr->getOpPtrVector().push_back(
496 field_eval_fe_ptr->getOpPtrVector().push_back(
498 field_eval_fe_ptr->getOpPtrVector().push_back(
500 "U", tensorFieldPtr));
501 }
502
504}
505//! [Set up problem]
506
507//! [Create common data]
510
511 auto get_command_line_parameters = [&]() {
513 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
514 &default_young_modulus, PETSC_NULL);
515 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
516 &default_poisson_ratio, PETSC_NULL);
517 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-coeff_expansion",
518 &default_coeff_expansion, PETSC_NULL);
519 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-ref_temp", &ref_temp,
520 PETSC_NULL);
521
522 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-capacity",
523 &default_heat_capacity, PETSC_NULL);
524 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-conductivity",
525 &default_heat_conductivity, PETSC_NULL);
526
527 MOFEM_LOG("ThermoElastic", Sev::inform)
528 << "Young modulus " << default_young_modulus;
529 MOFEM_LOG("ThermoElastic", Sev::inform)
530 << "Poisson ratio " << default_poisson_ratio;
531 MOFEM_LOG("ThermoElastic", Sev::inform)
532 << "Coeff of expansion " << default_coeff_expansion;
533 MOFEM_LOG("ThermoElastic", Sev::inform)
534 << "Capacity " << default_heat_capacity;
535 MOFEM_LOG("ThermoElastic", Sev::inform)
536 << "Heat conductivity " << default_heat_conductivity;
537
538 MOFEM_LOG("ThermoElastic", Sev::inform)
539 << "Reference_temperature " << ref_temp;
540
542 };
543
544 CHKERR get_command_line_parameters();
545
547}
548//! [Create common data]
549
550//! [Boundary condition]
553
554 MOFEM_LOG("SYNC", Sev::noisy) << "bC";
556
558 auto bc_mng = mField.getInterface<BcManager>();
559
560 CHKERR bc_mng->removeBlockDOFsOnEntities<DisplacementCubitBcData>(
561 simple->getProblemName(), "U");
562 CHKERR bc_mng->pushMarkDOFsOnEntities<HeatFluxCubitBcData>(
563 simple->getProblemName(), "FLUX", false);
564
565 auto get_skin = [&]() {
566 Range body_ents;
567 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
568 Skinner skin(&mField.get_moab());
569 Range skin_ents;
570 CHKERR skin.find_skin(0, body_ents, false, skin_ents);
571 return skin_ents;
572 };
573
574 auto filter_flux_blocks = [&](auto skin) {
575 auto remove_cubit_blocks = [&](auto c) {
577 for (auto m :
578
580
581 ) {
582 Range ents;
583 CHKERR mField.get_moab().get_entities_by_dimension(
584 m->getMeshset(), SPACE_DIM - 1, ents, true);
585 skin = subtract(skin, ents);
586 }
588 };
589
590 auto remove_named_blocks = [&](auto n) {
593 std::regex(
594
595 (boost::format("%s(.*)") % n).str()
596
597 ))
598
599 ) {
600 Range ents;
601 CHKERR mField.get_moab().get_entities_by_dimension(
602 m->getMeshset(), SPACE_DIM - 1, ents, true);
603 skin = subtract(skin, ents);
604 }
606 };
607
608 CHK_THROW_MESSAGE(remove_cubit_blocks(NODESET | TEMPERATURESET),
609 "remove_cubit_blocks");
610 CHK_THROW_MESSAGE(remove_cubit_blocks(SIDESET | HEATFLUXSET),
611 "remove_cubit_blocks");
612 CHK_THROW_MESSAGE(remove_named_blocks("TEMPERATURE"),
613 "remove_named_blocks");
614 CHK_THROW_MESSAGE(remove_named_blocks("HEATFLUX"), "remove_named_blocks");
615
616 return skin;
617 };
618
619 auto filter_true_skin = [&](auto skin) {
620 Range boundary_ents;
621 ParallelComm *pcomm =
622 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
623 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
624 PSTATUS_NOT, -1, &boundary_ents);
625 return boundary_ents;
626 };
627
628 auto remove_flux_ents = filter_true_skin(filter_flux_blocks(get_skin()));
629
630 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
631 remove_flux_ents);
632
633 MOFEM_LOG("SYNC", Sev::noisy) << remove_flux_ents << endl;
635
636#ifndef NDEBUG
637
640 (boost::format("flux_remove_%d.vtk") % mField.get_comm_rank()).str(),
641 remove_flux_ents);
642
643#endif
644
645 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
646 simple->getProblemName(), "FLUX", remove_flux_ents);
647
649}
650//! [Boundary condition]
651
652//! [Push operators to pipeline]
655
656 MOFEM_LOG("SYNC", Sev::noisy) << "OPs";
658
659 auto pipeline_mng = mField.getInterface<PipelineManager>();
661 auto bc_mng = mField.getInterface<BcManager>();
662
663 auto boundary_marker =
664 bc_mng->getMergedBlocksMarker(vector<string>{"HEATFLUX"});
665
666 auto integration_rule = [](int, int, int approx_order) {
667 return 2 * approx_order;
668 };
669 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
670 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
671 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
672 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
673
674 auto block_params = boost::make_shared<BlockedParameters>();
675 auto mDPtr = block_params->getDPtr();
676 auto coeff_expansion_ptr = block_params->getCoeffExpansionPtr();
677 auto heat_conductivity_ptr = block_params->getHeatConductivityPtr();
678 auto heat_capacity_ptr = block_params->getHeatCapacityPtr();
679
680 // Default time scaling of BCs and sources, from command line arguments
681 auto time_scale = boost::make_shared<TimeScale>();
682
683 // Files which define scaling for separate variables, if provided
684 auto time_bodyforce_scaling =
685 boost::make_shared<TimeScale>("bodyforce_scale.txt");
686 auto time_heatsource_scaling =
687 boost::make_shared<TimeScale>("heatsource_scale.txt");
688 auto time_temperature_scaling =
689 boost::make_shared<TimeScale>("temperature_bc_scale.txt");
690 auto time_displacement_scaling =
691 boost::make_shared<TimeScale>("displacement_bc_scale.txt");
692 auto time_flux_scaling = boost::make_shared<TimeScale>("flux_bc_scale.txt");
693 auto time_force_scaling =
694 boost::make_shared<TimeScale>("force_bc_scale.txt");
695
696 auto add_domain_rhs_ops = [&](auto &pipeline) {
698 CHKERR addMatBlockOps(pipeline, "MAT_ELASTIC", "MAT_THERMAL", block_params,
699 Sev::inform);
701
702 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
703 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
704 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
705
706 auto vec_temp_ptr = boost::make_shared<VectorDouble>();
707 auto vec_temp_dot_ptr = boost::make_shared<VectorDouble>();
708 auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
709 auto vec_temp_div_ptr = boost::make_shared<VectorDouble>();
710
711 pipeline.push_back(new OpCalculateScalarFieldValues("T", vec_temp_ptr));
712 pipeline.push_back(
713 new OpCalculateScalarFieldValuesDot("T", vec_temp_dot_ptr));
715 "FLUX", vec_temp_div_ptr));
716 pipeline.push_back(
717 new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
718
720 "U", mat_grad_ptr));
721 pipeline.push_back(
722 new OpSymmetrizeTensor<SPACE_DIM>(mat_grad_ptr, mat_strain_ptr));
723 pipeline.push_back(new OpStressThermal(mat_strain_ptr, vec_temp_ptr, mDPtr,
724 coeff_expansion_ptr,
725 mat_stress_ptr));
726
727 pipeline.push_back(new OpSetBc("FLUX", true, boundary_marker));
728
729 pipeline.push_back(new OpInternalForceCauchy(
730 "U", mat_stress_ptr,
731 [](double, double, double) constexpr { return 1; }));
732
733 auto resistance = [heat_conductivity_ptr](const double, const double,
734 const double) {
735 return (1. / (*heat_conductivity_ptr));
736 };
737 // negative value is consequence of symmetric system
738 auto capacity = [heat_capacity_ptr](const double, const double,
739 const double) {
740 return -(*heat_capacity_ptr);
741 };
742 auto unity = [](const double, const double, const double) constexpr {
743 return -1.;
744 };
745 pipeline.push_back(new OpHdivFlux("FLUX", mat_flux_ptr, resistance));
746 pipeline.push_back(new OpHDivTemp("FLUX", vec_temp_ptr, unity));
747 pipeline.push_back(new OpBaseDivFlux("T", vec_temp_div_ptr, unity));
748 pipeline.push_back(new OpBaseDotT("T", vec_temp_dot_ptr, capacity));
749
750 pipeline.push_back(new OpUnSetBc("FLUX"));
751
752 // Set source terms
754 pipeline, mField, "T", {time_scale, time_heatsource_scaling},
755 "HEAT_SOURCE", Sev::inform);
757 pipeline, mField, "U", {time_scale, time_bodyforce_scaling},
758 "BODY_FORCE", Sev::inform);
760 pipeline, mField, "T", vec_temp_ptr, "SETTEMP", Sev::inform);
761
763 };
764
765 auto add_domain_lhs_ops = [&](auto &pipeline) {
767 CHKERR addMatBlockOps(pipeline, "MAT_ELASTIC", "MAT_THERMAL", block_params,
768 Sev::verbose);
770
771 pipeline.push_back(new OpSetBc("FLUX", true, boundary_marker));
772
773 pipeline.push_back(new OpKCauchy("U", "U", mDPtr));
774 pipeline.push_back(new ThermoElasticOps::OpKCauchyThermoElasticity(
775 "U", "T", mDPtr, coeff_expansion_ptr));
776
777 auto resistance = [heat_conductivity_ptr](const double, const double,
778 const double) {
779 return (1. / (*heat_conductivity_ptr));
780 };
781 auto capacity = [heat_capacity_ptr](const double, const double,
782 const double) {
783 return -(*heat_capacity_ptr);
784 };
785 pipeline.push_back(new OpHdivHdiv("FLUX", "FLUX", resistance));
786 pipeline.push_back(new OpHdivT(
787 "FLUX", "T", []() constexpr { return -1; }, true));
788
789 auto op_capacity = new OpCapacity("T", "T", capacity);
790 op_capacity->feScalingFun = [](const FEMethod *fe_ptr) {
791 return fe_ptr->ts_a;
792 };
793 pipeline.push_back(op_capacity);
794
795 pipeline.push_back(new OpUnSetBc("FLUX"));
796
797 auto vec_temp_ptr = boost::make_shared<VectorDouble>();
798 pipeline.push_back(new OpCalculateScalarFieldValues("T", vec_temp_ptr));
800 pipeline, mField, "T", vec_temp_ptr, "SETTEMP", Sev::verbose);
801
803 };
804
805 auto add_boundary_rhs_ops = [&](auto &pipeline) {
807
809
810 pipeline.push_back(new OpSetBc("FLUX", true, boundary_marker));
811
812 // Set BCs using the least squares imposition
814 pipeline_mng->getOpBoundaryRhsPipeline(), mField, "U",
815 {time_scale, time_force_scaling}, "FORCE", Sev::inform);
816
818 pipeline_mng->getOpBoundaryRhsPipeline(), mField, "FLUX",
819 {time_scale, time_temperature_scaling}, "TEMPERATURE", Sev::inform);
820
821 pipeline.push_back(new OpUnSetBc("FLUX"));
822
823 auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
824 pipeline.push_back(
825 new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
828 mField, pipeline, simple->getProblemName(), "FLUX", mat_flux_ptr,
829 {time_scale, time_flux_scaling});
830
832 };
833
834 auto add_boundary_lhs_ops = [&](auto &pipeline) {
836
838
841 mField, pipeline, simple->getProblemName(), "FLUX");
842
844 };
845
846 // Set BCs by eliminating degrees of freedom
847 auto get_bc_hook_rhs = [&]() {
849 mField, pipeline_mng->getDomainRhsFE(),
850 {time_scale, time_displacement_scaling});
851 return hook;
852 };
853 auto get_bc_hook_lhs = [&]() {
855 mField, pipeline_mng->getDomainLhsFE(),
856 {time_scale, time_displacement_scaling});
857 return hook;
858 };
859
860 pipeline_mng->getDomainRhsFE()->preProcessHook = get_bc_hook_rhs();
861 pipeline_mng->getDomainLhsFE()->preProcessHook = get_bc_hook_lhs();
862
863 CHKERR add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
864 CHKERR add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
865 CHKERR add_boundary_rhs_ops(pipeline_mng->getOpBoundaryRhsPipeline());
866 CHKERR add_boundary_lhs_ops(pipeline_mng->getOpBoundaryLhsPipeline());
867
869}
870//! [Push operators to pipeline]
871
872//! [Solve]
875
878 ISManager *is_manager = mField.getInterface<ISManager>();
879
880 auto dm = simple->getDM();
881 auto solver = pipeline_mng->createTSIM();
882 auto snes_ctx_ptr = getDMSnesCtx(dm);
883
884 auto set_section_monitor = [&](auto solver) {
886 SNES snes;
887 CHKERR TSGetSNES(solver, &snes);
888 CHKERR SNESMonitorSet(snes,
889 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
891 (void *)(snes_ctx_ptr.get()), nullptr);
893 };
894
895 auto create_post_process_element = [&]() {
896 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
897
898 auto block_params = boost::make_shared<BlockedParameters>();
899 auto mDPtr = block_params->getDPtr();
900 auto coeff_expansion_ptr = block_params->getCoeffExpansionPtr();
901
902 CHKERR addMatBlockOps(post_proc_fe->getOpPtrVector(), "MAT_ELASTIC",
903 "MAT_THERMAL", block_params, Sev::verbose);
904
906 post_proc_fe->getOpPtrVector(), {H1, HDIV});
907
908 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
909 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
910 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
911
912 auto vec_temp_ptr = boost::make_shared<VectorDouble>();
913 auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
914
915 post_proc_fe->getOpPtrVector().push_back(
916 new OpCalculateScalarFieldValues("T", vec_temp_ptr));
917 post_proc_fe->getOpPtrVector().push_back(
918 new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
919
920 auto u_ptr = boost::make_shared<MatrixDouble>();
921 post_proc_fe->getOpPtrVector().push_back(
923 post_proc_fe->getOpPtrVector().push_back(
925 mat_grad_ptr));
926 post_proc_fe->getOpPtrVector().push_back(
927 new OpSymmetrizeTensor<SPACE_DIM>(mat_grad_ptr, mat_strain_ptr));
928 post_proc_fe->getOpPtrVector().push_back(
929 new OpStressThermal(mat_strain_ptr, vec_temp_ptr, mDPtr,
930 coeff_expansion_ptr, mat_stress_ptr));
931
933
934 post_proc_fe->getOpPtrVector().push_back(
935
936 new OpPPMap(
937
938 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
939
940 {{"T", vec_temp_ptr}},
941
942 {{"U", u_ptr}, {"FLUX", mat_flux_ptr}},
943
944 {},
945
946 {{"STRAIN", mat_strain_ptr}, {"STRESS", mat_stress_ptr}}
947
948 )
949
950 );
951
952 return post_proc_fe;
953 };
954
955 auto monitor_ptr = boost::make_shared<FEMethod>();
956 auto post_proc_fe = create_post_process_element();
957
958 auto set_time_monitor = [&](auto dm, auto solver) {
960 monitor_ptr->preProcessHook = [&]() {
962
963 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
964 post_proc_fe,
965 monitor_ptr->getCacheWeakPtr());
966 CHKERR post_proc_fe->writeFile(
967 "out_" + boost::lexical_cast<std::string>(monitor_ptr->ts_step) +
968 ".h5m");
969
970 if (doEvalField) {
971 if constexpr (SPACE_DIM == 3) {
972 CHKERR mField.getInterface<FieldEvaluatorInterface>()
973 ->evalFEAtThePoint3D(
974 fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
975 simple->getDomainFEName(), fieldEvalData,
976 mField.get_comm_rank(), mField.get_comm_rank(), nullptr,
977 MF_EXIST, QUIET);
978 } else {
979 CHKERR mField.getInterface<FieldEvaluatorInterface>()
980 ->evalFEAtThePoint2D(
981 fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
982 simple->getDomainFEName(), fieldEvalData,
983 mField.get_comm_rank(), mField.get_comm_rank(), nullptr,
984 MF_EXIST, QUIET);
985 }
986
987 if (scalarFieldPtr->size()) {
988 auto t_temp = getFTensor0FromVec(*scalarFieldPtr);
989 MOFEM_LOG("ThermoElasticSync", Sev::inform)
990 << "Eval point T: " << t_temp;
991 }
992 if (vectorFieldPtr->size1()) {
994 auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vectorFieldPtr);
995 MOFEM_LOG("ThermoElasticSync", Sev::inform)
996 << "Eval point U magnitude: " << sqrt(t_disp(i) * t_disp(i));
997 }
998 if (tensorFieldPtr->size1()) {
1000 auto t_disp_grad =
1001 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*tensorFieldPtr);
1002 MOFEM_LOG("ThermoElasticSync", Sev::inform)
1003 << "Eval point U_GRAD trace: " << t_disp_grad(i, i);
1004 }
1005
1006 MOFEM_LOG_SYNCHRONISE(mField.get_comm());
1007 }
1008
1010 };
1011 auto null = boost::shared_ptr<FEMethod>();
1012 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(), null,
1013 monitor_ptr, null);
1015 };
1016
1017 auto set_fieldsplit_preconditioner = [&](auto solver) {
1019
1020 SNES snes;
1021 CHKERR TSGetSNES(solver, &snes);
1022 KSP ksp;
1023 CHKERR SNESGetKSP(snes, &ksp);
1024 PC pc;
1025 CHKERR KSPGetPC(ksp, &pc);
1026 PetscBool is_pcfs = PETSC_FALSE;
1027 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1028
1029 // Setup fieldsplit (block) solver - optional: yes/no
1030 if (is_pcfs == PETSC_TRUE) {
1031 auto bc_mng = mField.getInterface<BcManager>();
1032 auto is_mng = mField.getInterface<ISManager>();
1033 auto name_prb = simple->getProblemName();
1034
1035 SmartPetscObj<IS> is_u;
1036 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "U", 0,
1037 SPACE_DIM, is_u);
1038 SmartPetscObj<IS> is_flux;
1039 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "FLUX", 0, 0,
1040 is_flux);
1041 SmartPetscObj<IS> is_T;
1042 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "T", 0, 0,
1043 is_T);
1044 IS is_tmp;
1045 CHKERR ISExpand(is_T, is_flux, &is_tmp);
1046 auto is_TFlux = SmartPetscObj<IS>(is_tmp);
1047
1048 auto is_all_bc = bc_mng->getBlockIS(name_prb, "HEATFLUX", "FLUX", 0, 0);
1049 int is_all_bc_size;
1050 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
1051 MOFEM_LOG("ThermoElastic", Sev::inform)
1052 << "Field split block size " << is_all_bc_size;
1053 if (is_all_bc_size) {
1054 IS is_tmp2;
1055 CHKERR ISDifference(is_TFlux, is_all_bc, &is_tmp2);
1056 is_TFlux = SmartPetscObj<IS>(is_tmp2);
1057 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
1058 is_all_bc); // boundary block
1059 }
1060
1061 CHKERR ISSort(is_u);
1062 CHKERR ISSort(is_TFlux);
1063 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_TFlux);
1064 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
1065 }
1066
1068 };
1069
1070 auto D = createDMVector(dm);
1071 CHKERR TSSetSolution(solver, D);
1072 CHKERR TSSetFromOptions(solver);
1073
1074 CHKERR set_section_monitor(solver);
1075 CHKERR set_fieldsplit_preconditioner(solver);
1076 CHKERR set_time_monitor(dm, solver);
1077
1078 CHKERR TSSetUp(solver);
1079 CHKERR TSSolve(solver, NULL);
1080
1082}
1083//! [Solve]
1084
1085static char help[] = "...\n\n";
1086
1087int main(int argc, char *argv[]) {
1088
1089 // Initialisation of MoFEM/PETSc and MOAB data structures
1090 const char param_file[] = "param_file.petsc";
1091 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
1092
1093 // Add logging channel for example
1094 auto core_log = logging::core::get();
1095 core_log->add_sink(
1097 LogManager::setLog("ThermoElastic");
1098 MOFEM_LOG_TAG("ThermoElastic", "ThermoElastic");
1099
1100 core_log->add_sink(
1101 LogManager::createSink(LogManager::getStrmSync(), "ThermoElasticSync"));
1102 LogManager::setLog("ThermoElasticSync");
1103 MOFEM_LOG_TAG("ThermoElasticSync", "ThermoElasticSync");
1104
1105 try {
1106
1107 //! [Register MoFEM discrete manager in PETSc]
1108 DMType dm_name = "DMMOFEM";
1109 CHKERR DMRegister_MoFEM(dm_name);
1110 //! [Register MoFEM discrete manager in PETSc
1111
1112 //! [Create MoAB]
1113 moab::Core mb_instance; ///< mesh database
1114 moab::Interface &moab = mb_instance; ///< mesh database interface
1115 //! [Create MoAB]
1116
1117 //! [Create MoFEM]
1118 MoFEM::Core core(moab); ///< finite element database
1119 MoFEM::Interface &m_field = core; ///< finite element database interface
1120 //! [Create MoFEM]
1121
1122 //! [Load mesh]
1123 Simple *simple = m_field.getInterface<Simple>();
1124 CHKERR simple->getOptions();
1125 CHKERR simple->loadFile();
1126 //! [Load mesh]
1127
1128 //! [ThermoElasticProblem]
1129 ThermoElasticProblem ex(m_field);
1130 CHKERR ex.runProblem();
1131 //! [ThermoElasticProblem]
1132 }
1134
1136}
std::string param_file
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
Definition: LogManager.hpp:352
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
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
constexpr int SPACE_DIM
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
@ QUIET
Definition: definitions.h:208
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ MF_EXIST
Definition: definitions.h:100
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 CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ NOSPACE
Definition: definitions.h:83
@ 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
@ TEMPERATURESET
Definition: definitions.h:155
@ HEATFLUXSET
Definition: definitions.h:156
@ NODESET
Definition: definitions.h:146
@ SIDESET
Definition: definitions.h:147
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#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
double bulk_modulus_K
double shear_modulus_G
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
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
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
FTensor::Index< 'i', SPACE_DIM > i
PetscBool doEvalField
const double c
speed of light (cm/ns)
double D
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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:232
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1857
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
double young_modulus
Young modulus.
Definition: plastic.cpp:172
constexpr auto size_symm
Definition: plastic.cpp:42
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
Data on single entity (This is passed as argument to DataOperator::doWork)
Essential boundary conditions.
Definition: Essential.hpp:111
Class (Function) to enforce essential constrains.
Definition: Essential.hpp:25
structure for User Loop Methods on finite elements
Field evaluator interface.
@ OPSPACE
operator do Work is execute on space data
Definition of the heat flux bc data structure.
Definition: BCData.hpp:427
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
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
Definition: LogManager.cpp:348
Interface for managing meshsets containing materials and boundary conditions.
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:81
Enforce essential constrains on rhs.
Definition: Essential.hpp:65
Post post-proc data at points from hash maps.
Set indices on entities on finite element.
PipelineManager interface.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
boost::shared_ptr< FieldEvaluatorInterface::SetPtsData > fieldEvalData
MoFEMErrorCode runProblem()
[Run problem]
std::array< double, SPACE_DIM > fieldEvalCoords
MoFEM::Interface & mField
MoFEMErrorCode createCommonData()
[Set up problem]
boost::shared_ptr< MatrixDouble > vectorFieldPtr
ThermoElasticProblem(MoFEM::Interface &m_field)
MoFEMErrorCode addMatBlockOps(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string block_elastic_name, std::string block_thermal_name, boost::shared_ptr< BlockedParameters > blockedParamsPtr, Sev sev)
boost::shared_ptr< VectorDouble > scalarFieldPtr
MoFEMErrorCode setupProblem()
add fields
MoFEMErrorCode bC()
[Create common data]
MoFEMErrorCode OPs()
[Boundary condition]
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
boost::shared_ptr< MatrixDouble > tensorFieldPtr
auto save_range
static char help[]
[Solve]
#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
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)
double default_heat_capacity
double default_young_modulus
[Essential boundary conditions (Least square approach)]
double default_coeff_expansion
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 >::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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 1, SPACE_DIM > OpHDivTemp
Integrate Rhs div flux base times temperature (T)
double default_heat_conductivity
OpBaseDotT OpBaseDivFlux
Integrate Rhs base of temperature times divergence of flux (T)
double default_poisson_ratio