v0.15.0
Loading...
Searching...
No Matches
ThermoElasticProblem Struct Reference
Collaboration diagram for ThermoElasticProblem:
[legend]

Classes

struct  BlockedThermalParameters
 
struct  BlockedThermoElasticParameters
 

Public Member Functions

 ThermoElasticProblem (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 [Run problem]
 

Private Member Functions

MoFEMErrorCode setupProblem ()
 add fields
 
MoFEMErrorCode getCommandLineParameters ()
 [Run problem]
 
MoFEMErrorCode bC ()
 [Set up problem]
 
MoFEMErrorCode OPs ()
 [Boundary condition]
 
MoFEMErrorCode tsSolve ()
 [Push operators to pipeline]
 
MoFEMErrorCode addMatThermalBlockOps (boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string block_name, boost::shared_ptr< BlockedThermalParameters > blockedParamsPtr, Sev sev)
 
MoFEMErrorCode addMatThermoElasticBlockOps (boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string block_name, boost::shared_ptr< BlockedThermoElasticParameters > blockedParamsPtr, Sev sev)
 
template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >
MoFEMErrorCode opThermoElasticFactoryDomainRhs (MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HenckyOps::CommonData > elastic_common_ptr, boost::shared_ptr< ThermoElasticProblem::BlockedThermalParameters > thermal_common_ptr, boost::shared_ptr< ThermoElasticProblem::BlockedThermoElasticParameters > thermoelastic_common_ptr, Sev sev)
 
template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >
MoFEMErrorCode opThermoElasticFactoryDomainRhs (MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string elastic_block_name, std::string thermal_block_name, std::string thermoelastic_block_name, Sev sev, double scale=1)
 
template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >
MoFEMErrorCode opThermoElasticFactoryDomainLhs (MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string coupled_field_name, boost::shared_ptr< HenckyOps::CommonData > elastic_common_ptr, boost::shared_ptr< ThermoElasticProblem::BlockedThermalParameters > thermal_common_ptr, boost::shared_ptr< ThermoElasticProblem::BlockedThermoElasticParameters > thermoelastic_common_ptr, Sev sev)
 
template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >
MoFEMErrorCode opThermoElasticFactoryDomainLhs (MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string coupled_field_name, std::string elastic_block_name, std::string thermal_block_name, std::string thermoelastic_block_name, Sev sev, double scale=1)
 

Private Attributes

MoFEM::InterfacemField
 
PetscBool doEvalField
 
std::array< double, 3 > fieldEvalCoords = {0.0, 0.0, 0.0}
 
boost::shared_ptr< FieldEvaluatorInterface::SetPtsDatafieldEvalData
 
boost::shared_ptr< VectorDoubletempFieldPtr
 
boost::shared_ptr< MatrixDoublefluxFieldPtr
 
boost::shared_ptr< MatrixDoubledispFieldPtr
 
boost::shared_ptr< MatrixDoubledispGradPtr
 
boost::shared_ptr< MatrixDoublestrainFieldPtr
 
boost::shared_ptr< MatrixDoublestressFieldPtr
 

Detailed Description

Examples
thermo_elastic.cpp.

Definition at line 130 of file thermo_elastic.cpp.

Constructor & Destructor Documentation

◆ ThermoElasticProblem()

ThermoElasticProblem::ThermoElasticProblem ( MoFEM::Interface & m_field)
inline
Examples
thermo_elastic.cpp.

Definition at line 132 of file thermo_elastic.cpp.

132: mField(m_field) {}
MoFEM::Interface & mField

Member Function Documentation

◆ addMatThermalBlockOps()

MoFEMErrorCode ThermoElasticProblem::addMatThermalBlockOps ( boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > & pipeline,
std::string block_name,
boost::shared_ptr< BlockedThermalParameters > blockedParamsPtr,
Sev sev )
private
Examples
thermo_elastic.cpp.

Definition at line 318 of file thermo_elastic.cpp.

321 {
323
324 struct OpMatThermalBlocks : public DomainEleOp {
325 OpMatThermalBlocks(boost::shared_ptr<double> conductivity_ptr,
326 boost::shared_ptr<double> capacity_ptr,
327 MoFEM::Interface &m_field, Sev sev,
328 std::vector<const CubitMeshSets *> meshset_vec_ptr)
329 : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE),
330 conductivityPtr(conductivity_ptr), capacityPtr(capacity_ptr) {
331 CHK_THROW_MESSAGE(extractThermalBlockData(m_field, meshset_vec_ptr, sev),
332 "Cannot get data from thermal block");
333 }
334
335 MoFEMErrorCode doWork(int side, EntityType type,
338
339 for (auto &b : blockData) {
340
341 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
342 *conductivityPtr = b.conductivity;
343 *capacityPtr = b.capacity;
345 }
346 }
347
348 *conductivityPtr = default_heat_conductivity;
349 *capacityPtr = default_heat_capacity;
350
352 }
353
354 private:
355 struct BlockData {
356 double conductivity;
357 double capacity;
358 Range blockEnts;
359 };
360
361 std::vector<BlockData> blockData;
362
364 extractThermalBlockData(MoFEM::Interface &m_field,
365 std::vector<const CubitMeshSets *> meshset_vec_ptr,
366 Sev sev) {
368
369 for (auto m : meshset_vec_ptr) {
370 MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Thermal Block") << *m;
371 std::vector<double> block_data;
372 CHKERR m->getAttributes(block_data);
373 if (block_data.size() < 2) {
374 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
375 "Expected that block has at least two attributes");
376 }
377 auto get_block_ents = [&]() {
378 Range ents;
379 CHKERR
380 m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
381 return ents;
382 };
383
384 MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Thermal Block")
385 << m->getName() << ": conductivity = " << block_data[0]
386 << " capacity = " << block_data[1];
387
388 blockData.push_back({block_data[0], block_data[1], get_block_ents()});
389 }
390 MOFEM_LOG_CHANNEL("WORLD");
392 }
393
394 boost::shared_ptr<double> conductivityPtr;
395 boost::shared_ptr<double> capacityPtr;
396 };
397
398 pipeline.push_back(new OpMatThermalBlocks(
399 blockedParamsPtr->getHeatConductivityPtr(),
400 blockedParamsPtr->getHeatCapacityPtr(), mField, sev,
401
402 // Get blockset using regular expression
404
405 (boost::format("%s(.*)") % block_name).str()
406
407 ))
408
409 ));
410
412}
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ NOSPACE
Definition definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
SeverityLevel
Severity levels.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
FTensor::Index< 'm', 3 > m
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
double default_heat_capacity
double default_heat_conductivity

◆ addMatThermoElasticBlockOps()

MoFEMErrorCode ThermoElasticProblem::addMatThermoElasticBlockOps ( boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > & pipeline,
std::string block_name,
boost::shared_ptr< BlockedThermoElasticParameters > blockedParamsPtr,
Sev sev )
private
Examples
thermo_elastic.cpp.

Definition at line 414 of file thermo_elastic.cpp.

418 {
420
421 struct OpMatThermoElasticBlocks : public DomainEleOp {
422 OpMatThermoElasticBlocks(boost::shared_ptr<VectorDouble> expansion_ptr,
423 boost::shared_ptr<double> ref_temp_ptr,
424 MoFEM::Interface &m_field, Sev sev,
425 std::vector<const CubitMeshSets *> meshset_vec_ptr)
426 : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE),
427 expansionPtr(expansion_ptr), refTempPtr(ref_temp_ptr) {
429 extractThermoElasticBlockData(m_field, meshset_vec_ptr, sev),
430 "Cannot get data from thermoelastic block");
431 }
432
433 MoFEMErrorCode doWork(int side, EntityType type,
436
437 for (auto &b : blockData) {
438
439 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
440 *expansionPtr = b.expansion;
441 *refTempPtr = b.ref_temp;
443 }
444 }
445
447 *refTempPtr = default_ref_temp;
448
450 }
451
452 private:
453 struct BlockData {
454 double ref_temp;
455 VectorDouble expansion;
456 Range blockEnts;
457 };
458
459 std::vector<BlockData> blockData;
460
461 MoFEMErrorCode extractThermoElasticBlockData(
462 MoFEM::Interface &m_field,
463 std::vector<const CubitMeshSets *> meshset_vec_ptr, Sev sev) {
465
466 for (auto m : meshset_vec_ptr) {
467 MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Thermoelastic Block") << *m;
468 std::vector<double> block_data;
469 CHKERR m->getAttributes(block_data);
470 if (block_data.size() < 2) {
471 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
472 "Expected that block has at least two attributes");
473 }
474 auto get_block_ents = [&]() {
475 Range ents;
476 CHKERR
477 m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
478 return ents;
479 };
480
481 auto get_expansion = [&]() {
482 VectorDouble expansion(SPACE_DIM, block_data[1]);
483 if (block_data.size() > 2) {
484 expansion[1] = block_data[2];
485 }
486 if (SPACE_DIM == 3 && block_data.size() > 3) {
487 expansion[2] = block_data[3];
488 }
489 return expansion;
490 };
491
492 auto coeff_exp_vec = get_expansion();
493
494 MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Thermoelastic Block")
495 << " ref_temp = " << block_data[0]
496 << " expansion = " << coeff_exp_vec;
497
498 blockData.push_back({block_data[0], coeff_exp_vec, get_block_ents()});
499 }
500 MOFEM_LOG_CHANNEL("WORLD");
502 }
503
504 boost::shared_ptr<VectorDouble> expansionPtr;
505 boost::shared_ptr<double> refTempPtr;
506 };
507
508 pipeline.push_back(new OpMatThermoElasticBlocks(
509 blockedParamsPtr->getCoeffExpansionPtr(),
510 blockedParamsPtr->getRefTempPtr(), mField, sev,
511
512 // Get blockset using regular expression
514
515 (boost::format("%s(.*)") % block_name).str()
516
517 ))
518
519 ));
520
522}
constexpr int SPACE_DIM
UBlasVector< double > VectorDouble
Definition Types.hpp:68
double default_ref_temp
double default_coeff_expansion

◆ bC()

MoFEMErrorCode ThermoElasticProblem::bC ( )
private

[Set up problem]

[Boundary condition]

Examples
thermo_elastic.cpp.

Definition at line 723 of file thermo_elastic.cpp.

723 {
725
726 MOFEM_LOG("SYNC", Sev::noisy) << "bC";
728
730 auto bc_mng = mField.getInterface<BcManager>();
731
732 auto get_skin = [&]() {
733 Range body_ents;
734 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
735 Skinner skin(&mField.get_moab());
736 Range skin_ents;
737 CHKERR skin.find_skin(0, body_ents, false, skin_ents);
738 return skin_ents;
739 };
740
741 auto filter_flux_blocks = [&](auto skin, bool temp_bc = false) {
742 auto remove_cubit_blocks = [&](auto c) {
744 for (auto m :
745
747
748 ) {
749 Range ents;
750 CHKERR mField.get_moab().get_entities_by_dimension(
751 m->getMeshset(), SPACE_DIM - 1, ents, true);
752 skin = subtract(skin, ents);
753 }
755 };
756
757 auto remove_named_blocks = [&](auto n) {
760 std::regex(
761
762 (boost::format("%s(.*)") % n).str()
763
764 ))
765
766 ) {
767 Range ents;
768 CHKERR mField.get_moab().get_entities_by_dimension(
769 m->getMeshset(), SPACE_DIM - 1, ents, true);
770 skin = subtract(skin, ents);
771 }
773 };
774 if (!temp_bc) {
775 CHK_THROW_MESSAGE(remove_cubit_blocks(NODESET | TEMPERATURESET),
776 "remove_cubit_blocks");
777 CHK_THROW_MESSAGE(remove_named_blocks("TEMPERATURE"),
778 "remove_named_blocks");
779 }
780 CHK_THROW_MESSAGE(remove_cubit_blocks(SIDESET | HEATFLUXSET),
781 "remove_cubit_blocks");
782 CHK_THROW_MESSAGE(remove_named_blocks("HEATFLUX"), "remove_named_blocks");
783 CHK_THROW_MESSAGE(remove_named_blocks("CONVECTION"), "remove_named_blocks");
784 CHK_THROW_MESSAGE(remove_named_blocks("RADIATION"), "remove_named_blocks");
785 return skin;
786 };
787
788 auto filter_true_skin = [&](auto skin) {
789 Range boundary_ents;
790 ParallelComm *pcomm =
791 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
792 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
793 PSTATUS_NOT, -1, &boundary_ents);
794 return boundary_ents;
795 };
796
797 auto remove_flux_ents = filter_true_skin(filter_flux_blocks(get_skin()));
798 auto remove_temp_bc_ents =
799 filter_true_skin(filter_flux_blocks(get_skin(), true));
800
801 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
802 remove_flux_ents);
803 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
804 remove_temp_bc_ents);
805
806 MOFEM_LOG("SYNC", Sev::noisy) << remove_flux_ents << endl;
808
809 MOFEM_LOG("SYNC", Sev::noisy) << remove_temp_bc_ents << endl;
811
812#ifndef NDEBUG
813
816 (boost::format("flux_remove_%d.vtk") % mField.get_comm_rank()).str(),
817 remove_flux_ents);
818
821 (boost::format("temp_bc_remove_%d.vtk") % mField.get_comm_rank()).str(),
822 remove_temp_bc_ents);
823
824#endif
825
826 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
827 simple->getProblemName(), "FLUX", remove_flux_ents);
828 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
829 simple->getProblemName(), "TBC", remove_temp_bc_ents);
830
831 auto set_init_temp = [](boost::shared_ptr<FieldEntity> field_entity_ptr) {
832 field_entity_ptr->getEntFieldData()[0] = default_init_temp;
833 return 0;
834 };
835 CHKERR mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(set_init_temp,
836 "T");
837
838 CHKERR bc_mng->removeBlockDOFsOnEntities<DisplacementCubitBcData>(
839 simple->getProblemName(), "U");
840 CHKERR bc_mng->pushMarkDOFsOnEntities<HeatFluxCubitBcData>(
841 simple->getProblemName(), "FLUX", false);
842
844}
static auto filter_true_skin(MoFEM::Interface &m_field, Range &&skin)
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
#define MYPCOMM_INDEX
default communicator number PCOMM
@ TEMPERATURESET
@ HEATFLUXSET
@ NODESET
@ SIDESET
#define MOFEM_LOG(channel, severity)
Log.
const double c
speed of light (cm/ns)
const double n
refractive index of diffusive medium
Simple interface for fast problem set-up.
Definition BcManager.hpp:29
Managing BitRefLevels.
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Definition of the displacement bc data structure.
Definition BCData.hpp:72
Basic algebra on fields.
Definition FieldBlas.hpp:21
Definition of the heat flux bc data structure.
Definition BCData.hpp:423
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
constexpr int SPACE_DIM
auto save_range
double default_init_temp

◆ getCommandLineParameters()

MoFEMErrorCode ThermoElasticProblem::getCommandLineParameters ( )
private

[Run problem]

[Get command line parameters]

Examples
thermo_elastic.cpp.

Definition at line 537 of file thermo_elastic.cpp.

537 {
539
540 auto get_command_line_parameters = [&]() {
542
543 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order_temp,
544 PETSC_NULLPTR);
545 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order_temp", &order_temp,
546 PETSC_NULLPTR);
548 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order_flux", &order_flux,
549 PETSC_NULLPTR);
551 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order_disp", &order_disp,
552 PETSC_NULLPTR);
553 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-atom_test", &atom_test,
554 PETSC_NULLPTR);
555 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-save_every", &save_every,
556 PETSC_NULLPTR);
557
558 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-young_modulus",
559 &default_young_modulus, PETSC_NULLPTR);
560 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-poisson_ratio",
561 &default_poisson_ratio, PETSC_NULLPTR);
562 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-coeff_expansion",
563 &default_coeff_expansion, PETSC_NULLPTR);
564 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-ref_temp", &default_ref_temp,
565 PETSC_NULLPTR);
566 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-init_temp",
567 &default_init_temp, PETSC_NULLPTR);
568
569 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-capacity",
570 &default_heat_capacity, PETSC_NULLPTR);
571 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-conductivity",
572 &default_heat_conductivity, PETSC_NULLPTR);
573
574 if constexpr (SPACE_DIM == 2) {
575 do_output_domain = PETSC_TRUE;
576 do_output_skin = PETSC_FALSE;
577 } else {
578 do_output_domain = PETSC_FALSE;
579 do_output_skin = PETSC_TRUE;
580 }
581
582 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-output_domain",
583 &do_output_domain, PETSC_NULLPTR);
584 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-output_skin", &do_output_skin,
585 PETSC_NULLPTR);
586
587 MOFEM_LOG("ThermoElastic", Sev::inform)
588 << "Young's modulus " << default_young_modulus;
589 MOFEM_LOG("ThermoElastic", Sev::inform)
590 << "Poisson's ratio " << default_poisson_ratio;
591 MOFEM_LOG("ThermoElastic", Sev::inform)
592 << "Coeff of expansion " << default_coeff_expansion;
593 MOFEM_LOG("ThermoElastic", Sev::inform)
594 << "Default heat capacity " << default_heat_capacity;
595 MOFEM_LOG("ThermoElastic", Sev::inform)
596 << "Heat conductivity " << default_heat_conductivity;
597
598 MOFEM_LOG("ThermoElastic", Sev::inform)
599 << "Reference temperature " << default_ref_temp;
600 MOFEM_LOG("ThermoElastic", Sev::inform)
601 << "Initial temperature " << default_init_temp;
602
604 };
605
606 CHKERR get_command_line_parameters();
607
609}
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)
int save_every
int order_flux
int atom_test
int order_disp
double default_young_modulus
[Essential boundary conditions (Least square approach)]
PetscBool do_output_skin
int order_temp
PetscBool do_output_domain
double default_poisson_ratio

◆ OPs()

MoFEMErrorCode ThermoElasticProblem::OPs ( )
private

[Boundary condition]

[Push operators to pipeline]

Examples
thermo_elastic.cpp.

Definition at line 848 of file thermo_elastic.cpp.

848 {
850
851 MOFEM_LOG("SYNC", Sev::noisy) << "OPs";
853
854 auto pipeline_mng = mField.getInterface<PipelineManager>();
855 auto bc_mng = mField.getInterface<BcManager>();
856
857 auto boundary_marker =
858 bc_mng->getMergedBlocksMarker(vector<string>{"HEATFLUX"});
859
860 auto integration_rule = [](int, int, int approx_order) {
861 return 2 * approx_order;
862 };
864 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
865 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
866 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
867
868 auto block_thermal_params = boost::make_shared<BlockedThermalParameters>();
869 auto heat_conductivity_ptr = block_thermal_params->getHeatConductivityPtr();
870 auto heat_capacity_ptr = block_thermal_params->getHeatCapacityPtr();
871
872 auto block_thermoelastic_params =
873 boost::make_shared<BlockedThermoElasticParameters>();
874 auto coeff_expansion_ptr = block_thermoelastic_params->getCoeffExpansionPtr();
875 auto ref_temp_ptr = block_thermoelastic_params->getRefTempPtr();
876
877 // Default time scaling of BCs and sources, from command line arguments
878 auto time_scale =
879 boost::make_shared<TimeScale>("", false, [](const double) { return 1; });
880 auto def_time_scale = [time_scale](const double time) {
881 return (!time_scale->argFileScale) ? time : 1;
882 };
883 auto def_file_name = [time_scale](const std::string &&name) {
884 return (!time_scale->argFileScale) ? name : "";
885 };
886
887 // Files which define scaling for separate variables, if provided
888 auto time_bodyforce_scaling = boost::make_shared<TimeScale>(
889 def_file_name("bodyforce_scale.txt"), false, def_time_scale);
890 auto time_heatsource_scaling = boost::make_shared<TimeScale>(
891 def_file_name("heatsource_scale.txt"), false, def_time_scale);
892 auto time_temperature_scaling = boost::make_shared<TimeScale>(
893 def_file_name("temperature_bc_scale.txt"), false, def_time_scale);
894 auto time_displacement_scaling = boost::make_shared<TimeScale>(
895 def_file_name("displacement_bc_scale.txt"), false, def_time_scale);
896 auto time_flux_scaling = boost::make_shared<TimeScale>(
897 def_file_name("flux_bc_scale.txt"), false, def_time_scale);
898 auto time_force_scaling = boost::make_shared<TimeScale>(
899 def_file_name("force_bc_scale.txt"), false, def_time_scale);
900
901 auto add_domain_rhs_ops = [&](auto &pipeline) {
903 CHKERR addMatThermalBlockOps(pipeline, "MAT_THERMAL", block_thermal_params,
904 Sev::inform);
905 CHKERR addMatThermoElasticBlockOps(pipeline, "MAT_THERMOELASTIC",
906 block_thermoelastic_params, Sev::inform);
908
909 auto hencky_common_data_ptr =
911 mField, pipeline, "U", "MAT_ELASTIC", Sev::inform);
912 auto mat_D_ptr = hencky_common_data_ptr->matDPtr;
913 auto mat_grad_ptr = hencky_common_data_ptr->matGradPtr;
914 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
915 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
916
917 auto vec_temp_ptr = boost::make_shared<VectorDouble>();
918 auto vec_temp_dot_ptr = boost::make_shared<VectorDouble>();
919 auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
920 auto vec_temp_div_ptr = boost::make_shared<VectorDouble>();
921
922 pipeline.push_back(new OpCalculateScalarFieldValues("T", vec_temp_ptr));
923 pipeline.push_back(
924 new OpCalculateScalarFieldValuesDot("T", vec_temp_dot_ptr));
926 "FLUX", vec_temp_div_ptr));
927 pipeline.push_back(
928 new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
929
931 mField, pipeline, "U", "MAT_ELASTIC", "MAT_THERMAL",
932 "MAT_THERMOELASTIC", Sev::inform);
933
934 auto resistance = [heat_conductivity_ptr](const double, const double,
935 const double) {
936 return (1. / (*heat_conductivity_ptr));
937 };
938 // negative value is consequence of symmetric system
939 auto capacity = [heat_capacity_ptr](const double, const double,
940 const double) {
941 return -(*heat_capacity_ptr);
942 };
943 auto unity = [](const double, const double, const double) constexpr {
944 return -1.;
945 };
946 pipeline.push_back(
947 new OpHdivFlux("FLUX", mat_flux_ptr, resistance, mat_grad_ptr));
948 pipeline.push_back(new OpHDivTemp("FLUX", vec_temp_ptr, unity));
949 pipeline.push_back(new OpBaseDivFlux("T", vec_temp_div_ptr, unity));
950 pipeline.push_back(new OpBaseDotT("T", vec_temp_dot_ptr, capacity));
951
952 // Set source terms
954 pipeline, mField, "T", {time_scale, time_heatsource_scaling},
955 "HEAT_SOURCE", Sev::inform);
957 pipeline, mField, "U", {time_scale, time_bodyforce_scaling},
958 "BODY_FORCE", Sev::inform);
960 pipeline, mField, "T", vec_temp_ptr, "SETTEMP", Sev::inform);
961
963 };
964
965 auto add_domain_lhs_ops = [&](auto &pipeline) {
967 CHKERR addMatThermalBlockOps(pipeline, "MAT_THERMAL", block_thermal_params,
968 Sev::verbose);
969 CHKERR addMatThermoElasticBlockOps(pipeline, "MAT_THERMOELASTIC",
970 block_thermoelastic_params,
971 Sev::verbose);
973
974 auto hencky_common_data_ptr =
976 mField, pipeline, "U", "MAT_ELASTIC", Sev::inform, 1);
977 auto mat_D_ptr = hencky_common_data_ptr->matDPtr;
978 auto mat_grad_ptr = hencky_common_data_ptr->matGradPtr;
979
980 auto resistance = [heat_conductivity_ptr](const double, const double,
981 const double) {
982 return (1. / (*heat_conductivity_ptr));
983 };
984 auto capacity = [heat_capacity_ptr](const double, const double,
985 const double) {
986 return -(*heat_capacity_ptr);
987 };
988 pipeline.push_back(
989 new OpHdivHdiv("FLUX", "FLUX", resistance, mat_grad_ptr));
990 pipeline.push_back(
991 new OpHdivT("FLUX", "T", []() constexpr { return -1; }, true));
992
993 auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
994 pipeline.push_back(
995 new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
996
997 pipeline.push_back(
998 new OpHdivU("FLUX", "U", mat_flux_ptr, resistance, mat_grad_ptr));
999
1001 mField, pipeline, "U", "T", "MAT_ELASTIC", "MAT_THERMAL",
1002 "MAT_THERMOELASTIC", Sev::inform);
1003
1004 auto op_capacity = new OpCapacity("T", "T", capacity);
1005 op_capacity->feScalingFun = [](const FEMethod *fe_ptr) {
1006 return fe_ptr->ts_a;
1007 };
1008 pipeline.push_back(op_capacity);
1009
1010 auto vec_temp_ptr = boost::make_shared<VectorDouble>();
1011 pipeline.push_back(new OpCalculateScalarFieldValues("T", vec_temp_ptr));
1013 pipeline, mField, "T", vec_temp_ptr, "SETTEMP", Sev::verbose);
1014
1016 };
1017
1018 auto add_boundary_rhs_ops = [&](auto &pipeline) {
1020
1022
1024 pipeline, mField, "U", {time_scale, time_force_scaling}, "FORCE",
1025 Sev::inform);
1026
1027 // Temperature BC
1028
1029 using OpTemperatureBC =
1032 pipeline, mField, "FLUX", {time_scale, time_temperature_scaling},
1033 "TEMPERATURE", Sev::inform);
1034
1035 // Note: fluxes for temperature should be aggregated. Here separate is
1036 // NaturalMeshsetType<HEATFLUXSET>, we should add version with BLOCKSET,
1037 // convection, see example tutorials/vec-0/src/NaturalBoundaryBC.hpp.
1038 // and vec-0/elastic.cpp
1039
1040 using OpFluxBC =
1043 pipeline, mField, "TBC", {time_scale, time_flux_scaling}, "FLUX",
1044 Sev::inform);
1045
1047 using OpConvectionRhsBC =
1048 T::OpFlux<ThermoElasticOps::ConvectionBcType<BLOCKSET>, 1, 1>;
1049 using OpRadiationRhsBC =
1050 T::OpFlux<ThermoElasticOps::RadiationBcType<BLOCKSET>, 1, 1>;
1051 auto temp_bc_ptr = boost::make_shared<VectorDouble>();
1052 pipeline.push_back(new OpCalculateScalarFieldValues("TBC", temp_bc_ptr));
1053 T::AddFluxToPipeline<OpConvectionRhsBC>::add(
1054 pipeline, mField, "TBC", temp_bc_ptr, "CONVECTION", Sev::inform);
1055 T::AddFluxToPipeline<OpRadiationRhsBC>::add(
1056 pipeline, mField, "TBC", temp_bc_ptr, "RADIATION", Sev::inform);
1057
1058 auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
1059 pipeline.push_back(
1060 new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
1061
1062 // This is temporary implementation. It be moved to LinearFormsIntegrators,
1063 // however for hybridised case, what is goal of this changes, such function
1064 // is implemented for fluxes in broken space. Thus ultimately this operator
1065 // would be not needed.
1066
1067 struct OpTBCTimesNormalFlux
1068 : public FormsIntegrators<BoundaryEleOp>::Assembly<PETSC>::OpBase {
1069
1071
1072 OpTBCTimesNormalFlux(const std::string field_name,
1073 boost::shared_ptr<MatrixDouble> vec,
1074 boost::shared_ptr<Range> ents_ptr = nullptr)
1075 : OP(field_name, field_name, OP::OPROW, ents_ptr), sourceVec(vec) {}
1076
1080 // get integration weights
1081 auto t_w = OP::getFTensor0IntegrationWeight();
1082 // get base function values on rows
1083 auto t_row_base = row_data.getFTensor0N();
1084 // get normal vector
1085 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1086 // get vector values
1087 auto t_vec = getFTensor1FromMat<SPACE_DIM, 1>(*sourceVec);
1088 // loop over integration points
1089 for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
1090 // take into account Jacobian
1091 const double alpha = t_w * (t_vec(i) * t_normal(i));
1092 // loop over rows base functions
1093 int rr = 0;
1094 for (; rr != OP::nbRows; ++rr) {
1095 OP::locF[rr] += alpha * t_row_base;
1096 ++t_row_base;
1097 }
1098 for (; rr < OP::nbRowBaseFunctions; ++rr)
1099 ++t_row_base;
1100 ++t_w; // move to another integration weight
1101 ++t_vec;
1102 ++t_normal;
1103 }
1104 EntityType fe_type = OP::getNumeredEntFiniteElementPtr()->getEntType();
1105 if (fe_type == MBTRI) {
1106 OP::locF /= 2;
1107 }
1109 }
1110
1111 protected:
1112 boost::shared_ptr<MatrixDouble> sourceVec;
1113 };
1114 pipeline.push_back(new OpTBCTimesNormalFlux("TBC", mat_flux_ptr));
1115
1116 struct OpBaseNormalTimesTBC
1117 : public FormsIntegrators<BoundaryEleOp>::Assembly<PETSC>::OpBase {
1118
1120
1121 OpBaseNormalTimesTBC(const std::string field_name,
1122 boost::shared_ptr<VectorDouble> vec,
1123 boost::shared_ptr<Range> ents_ptr = nullptr)
1124 : OP(field_name, field_name, OP::OPROW, ents_ptr), sourceVec(vec) {}
1125
1129 // get integration weights
1130 auto t_w = OP::getFTensor0IntegrationWeight();
1131 // get base function values on rows
1132 auto t_row_base = row_data.getFTensor1N<3>();
1133 // get normal vector
1134 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1135 // get vector values
1136 auto t_vec = getFTensor0FromVec(*sourceVec);
1137 // loop over integration points
1138 for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
1139 // take into account Jacobian
1140 const double alpha = t_w * t_vec;
1141 // loop over rows base functions
1142 int rr = 0;
1143 for (; rr != OP::nbRows; ++rr) {
1144 OP::locF[rr] += alpha * (t_row_base(i) * t_normal(i));
1145 ++t_row_base;
1146 }
1147 for (; rr < OP::nbRowBaseFunctions; ++rr)
1148 ++t_row_base;
1149 ++t_w; // move to another integration weight
1150 ++t_vec;
1151 ++t_normal;
1152 }
1153 EntityType fe_type = OP::getNumeredEntFiniteElementPtr()->getEntType();
1154 if (fe_type == MBTRI) {
1155 OP::locF /= 2;
1156 }
1158 }
1159
1160 protected:
1161 boost::shared_ptr<VectorDouble> sourceVec;
1162 };
1163
1164 pipeline.push_back(new OpBaseNormalTimesTBC("FLUX", temp_bc_ptr));
1165
1167 };
1168
1169 auto add_boundary_lhs_ops = [&](auto &pipeline) {
1171
1173
1174 using T =
1175 NaturalBC<BoundaryEleOp>::template Assembly<PETSC>::BiLinearForm<GAUSS>;
1176 using OpConvectionLhsBC =
1177 T::OpFlux<ThermoElasticOps::ConvectionBcType<BLOCKSET>, 1, 1>;
1178 using OpRadiationLhsBC =
1179 T::OpFlux<ThermoElasticOps::RadiationBcType<BLOCKSET>, 1, 1>;
1180 auto temp_bc_ptr = boost::make_shared<VectorDouble>();
1181 pipeline.push_back(new OpCalculateScalarFieldValues("TBC", temp_bc_ptr));
1182 T::AddFluxToPipeline<OpConvectionLhsBC>::add(pipeline, mField, "TBC", "TBC",
1183 "CONVECTION", Sev::verbose);
1184 T::AddFluxToPipeline<OpRadiationLhsBC>::add(
1185 pipeline, mField, "TBC", "TBC", temp_bc_ptr, "RADIATION", Sev::verbose);
1186
1187 // This is temporary implementation. It be moved to LinearFormsIntegrators,
1188 // however for hybridised case, what is goal of this changes, such function
1189 // is implemented for fluxes in broken space. Thus ultimately this operator
1190 // would be not needed.
1191
1192 struct OpTBCTimesNormalFlux
1193 : public FormsIntegrators<BoundaryEleOp>::Assembly<PETSC>::OpBase {
1194
1196
1197 OpTBCTimesNormalFlux(const std::string row_field_name,
1198 const std::string col_field_name,
1199 boost::shared_ptr<Range> ents_ptr = nullptr)
1200 : OP(row_field_name, col_field_name, OP::OPROWCOL, ents_ptr) {
1201 this->sYmm = false;
1202 this->assembleTranspose = true;
1203 this->onlyTranspose = false;
1204 }
1205
1207 EntitiesFieldData::EntData &col_data) {
1209
1211
1212 // get integration weights
1213 auto t_w = OP::getFTensor0IntegrationWeight();
1214 // get base function values on rows
1215 auto t_row_base = row_data.getFTensor0N();
1216 // get normal vector
1217 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1218 // loop over integration points
1219 for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
1220 // loop over rows base functions
1221 auto a_mat_ptr = &*OP::locMat.data().begin();
1222 int rr = 0;
1223 for (; rr != OP::nbRows; rr++) {
1224 // take into account Jacobian
1225 const double alpha = t_w * t_row_base;
1226 // get column base functions values at gauss point gg
1227 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
1228 // loop over columns
1229 for (int cc = 0; cc != OP::nbCols; cc++) {
1230 // calculate element of local matrix
1231 // using L2 norm of normal vector for convenient area calculation
1232 // for quads, tris etc.
1233 *a_mat_ptr += alpha * (t_col_base(i) * t_normal(i));
1234 ++t_col_base;
1235 ++a_mat_ptr;
1236 }
1237 ++t_row_base;
1238 }
1239 for (; rr < OP::nbRowBaseFunctions; ++rr)
1240 ++t_row_base;
1241 ++t_normal;
1242 ++t_w; // move to another integration weight
1243 }
1244 EntityType fe_type = OP::getNumeredEntFiniteElementPtr()->getEntType();
1245 if (fe_type == MBTRI) {
1246 OP::locMat /= 2;
1247 }
1249 }
1250 };
1251
1252 pipeline.push_back(new OpTBCTimesNormalFlux("TBC", "FLUX"));
1253
1255 };
1256
1257 // Set BCs by eliminating degrees of freedom
1258 auto get_bc_hook_rhs = [&]() {
1260 mField, pipeline_mng->getDomainRhsFE(),
1261 {time_scale, time_displacement_scaling});
1262 return hook;
1263 };
1264 auto get_bc_hook_lhs = [&]() {
1266 mField, pipeline_mng->getDomainLhsFE(),
1267 {time_scale, time_displacement_scaling});
1268 return hook;
1269 };
1270
1271 pipeline_mng->getDomainRhsFE()->preProcessHook = get_bc_hook_rhs();
1272 pipeline_mng->getDomainLhsFE()->preProcessHook = get_bc_hook_lhs();
1273
1274 CHKERR add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
1275 CHKERR add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
1276 CHKERR add_boundary_rhs_ops(pipeline_mng->getOpBoundaryRhsPipeline());
1277 CHKERR add_boundary_lhs_ops(pipeline_mng->getOpBoundaryLhsPipeline());
1278
1280}
#define FTENSOR_INDEX(DIM, I)
@ H1
continuous field
Definition definitions.h:85
@ HDIV
field with continuous normal traction
Definition definitions.h:87
auto integration_rule
FTensor::Index< 'i', SPACE_DIM > i
auto commonDataFactory(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string block_name, Sev sev, double scale=1)
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
OpCalculateQdotQLhs_dU< SPACE_DIM, GAUSS, AssemblyDomainEleOp, IS_LARGE_STRAINS > OpHdivU
Integrate Lhs of flux term coupled to displacement field.
OpBaseDotT OpBaseDivFlux
Integrate Rhs base of temperature times divergence of flux (T)
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 >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 1, SPACE_DIM > OpHDivTemp
Integrate Rhs div flux base times temperature (T)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< SPACE_DIM > OpHdivT
Integrate Lhs div of base of flux times base of temperature (FLUX x T) and transpose of it,...
OpHdivFluxImpl< SPACE_DIM, IS_LARGE_STRAINS > OpHdivFlux
OpHdivHdivImpl< SPACE_DIM, IS_LARGE_STRAINS > OpHdivHdiv
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBaseDotT
Integrate Rhs base of temperature time heat capacity times heat rate (T)
constexpr auto field_name
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition radiation.cpp:29
Add operators pushing bases from local to physical configuration.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Class (Function) to enforce essential constrains.
Definition Essential.hpp:25
structure for User Loop Methods on finite elements
Assembly methods.
Definition Natural.hpp:65
Natural boundary conditions.
Definition Natural.hpp:57
virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate grad-grad operator.
Get vector field for H-div approximation.
Calculate divergence of vector field.
Get value at integration points for scalar field.
PipelineManager interface.
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode opThermoElasticFactoryDomainRhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HenckyOps::CommonData > elastic_common_ptr, boost::shared_ptr< ThermoElasticProblem::BlockedThermalParameters > thermal_common_ptr, boost::shared_ptr< ThermoElasticProblem::BlockedThermoElasticParameters > thermoelastic_common_ptr, Sev sev)
MoFEMErrorCode opThermoElasticFactoryDomainLhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string coupled_field_name, boost::shared_ptr< HenckyOps::CommonData > elastic_common_ptr, boost::shared_ptr< ThermoElasticProblem::BlockedThermalParameters > thermal_common_ptr, boost::shared_ptr< ThermoElasticProblem::BlockedThermoElasticParameters > thermoelastic_common_ptr, Sev sev)
MoFEMErrorCode addMatThermalBlockOps(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string block_name, boost::shared_ptr< BlockedThermalParameters > blockedParamsPtr, Sev sev)
MoFEMErrorCode addMatThermoElasticBlockOps(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string block_name, boost::shared_ptr< BlockedThermoElasticParameters > blockedParamsPtr, Sev sev)

◆ opThermoElasticFactoryDomainLhs() [1/2]

template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >
MoFEMErrorCode ThermoElasticProblem::opThermoElasticFactoryDomainLhs ( MoFEM::Interface & m_field,
boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > & pip,
std::string field_name,
std::string coupled_field_name,
boost::shared_ptr< HenckyOps::CommonData > elastic_common_ptr,
boost::shared_ptr< ThermoElasticProblem::BlockedThermalParameters > thermal_common_ptr,
boost::shared_ptr< ThermoElasticProblem::BlockedThermoElasticParameters > thermoelastic_common_ptr,
Sev sev )
inlineprivate
Examples
thermo_elastic.cpp.

Definition at line 255 of file thermo_elastic.cpp.

264 {
266
267 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
268 A>::template BiLinearForm<I>;
269 using OpKPiola = typename B::template OpGradTensorGrad<1, DIM, DIM, 1>;
270
272 auto vec_temp_ptr = boost::make_shared<VectorDouble>();
273 pip.push_back(new OpCalculateScalarFieldValues("T", vec_temp_ptr));
274 auto coeff_expansion_ptr = thermoelastic_common_ptr->getCoeffExpansionPtr();
275 auto ref_temp_ptr = thermoelastic_common_ptr->getRefTempPtr();
276 pip.push_back(
277 new typename H::template OpCalculateHenckyThermalStress<DIM, I, 0>(
278 "U", vec_temp_ptr, elastic_common_ptr, coeff_expansion_ptr,
279 ref_temp_ptr));
280 pip.push_back(new typename H::template OpHenckyTangent<DIM, I, 0>(
281 field_name, elastic_common_ptr));
282 pip.push_back(new OpKPiola(field_name, field_name,
283 elastic_common_ptr->getMatTangent()));
284 pip.push_back(new typename H::template OpCalculateHenckyThermalStressdT<
285 DIM, I, AssemblyDomainEleOp, 0>(
286 field_name, coupled_field_name, elastic_common_ptr,
287 coeff_expansion_ptr));
288
290 }
constexpr IntegrationType I
double H
Hardening.
Definition plastic.cpp:128
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition seepage.cpp:62

◆ opThermoElasticFactoryDomainLhs() [2/2]

template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >
MoFEMErrorCode ThermoElasticProblem::opThermoElasticFactoryDomainLhs ( MoFEM::Interface & m_field,
boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > & pip,
std::string field_name,
std::string coupled_field_name,
std::string elastic_block_name,
std::string thermal_block_name,
std::string thermoelastic_block_name,
Sev sev,
double scale = 1 )
inlineprivate

Definition at line 293 of file thermo_elastic.cpp.

298 {
300
301 auto elastic_common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
302 m_field, pip, field_name, elastic_block_name, sev, scale);
303 auto thermal_common_ptr = boost::make_shared<BlockedThermalParameters>();
304 CHKERR addMatThermalBlockOps(pip, thermal_block_name, thermal_common_ptr,
305 Sev::inform);
306 auto thermoelastic_common_ptr =
307 boost::make_shared<BlockedThermoElasticParameters>();
308 CHKERR addMatThermoElasticBlockOps(pip, thermoelastic_block_name,
309 thermoelastic_common_ptr, Sev::inform);
311 m_field, pip, field_name, coupled_field_name, elastic_common_ptr,
312 thermal_common_ptr, thermoelastic_common_ptr, sev);
313
315 }
double scale
Definition plastic.cpp:123

◆ opThermoElasticFactoryDomainRhs() [1/2]

template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >
MoFEMErrorCode ThermoElasticProblem::opThermoElasticFactoryDomainRhs ( MoFEM::Interface & m_field,
boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > & pip,
std::string field_name,
boost::shared_ptr< HenckyOps::CommonData > elastic_common_ptr,
boost::shared_ptr< ThermoElasticProblem::BlockedThermalParameters > thermal_common_ptr,
boost::shared_ptr< ThermoElasticProblem::BlockedThermoElasticParameters > thermoelastic_common_ptr,
Sev sev )
inlineprivate
Examples
thermo_elastic.cpp.

Definition at line 198 of file thermo_elastic.cpp.

207 {
209
210 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
211 A>::template LinearForm<I>;
213 auto vec_temp_ptr = boost::make_shared<VectorDouble>();
214 pip.push_back(new OpCalculateScalarFieldValues("T", vec_temp_ptr));
215 auto coeff_expansion_ptr = thermoelastic_common_ptr->getCoeffExpansionPtr();
216 auto ref_temp_ptr = thermoelastic_common_ptr->getRefTempPtr();
217 pip.push_back(
218 new typename H::template OpCalculateHenckyThermalStress<DIM, I, 0>(
219 "U", vec_temp_ptr, elastic_common_ptr, coeff_expansion_ptr,
220 ref_temp_ptr));
222 typename B::template OpGradTimesTensor<1, DIM, DIM>;
223 pip.push_back(new OpInternalForcePiola(
224 "U", elastic_common_ptr->getMatFirstPiolaStress()));
225
227 }
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition seepage.cpp:64

◆ opThermoElasticFactoryDomainRhs() [2/2]

template<int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp >
MoFEMErrorCode ThermoElasticProblem::opThermoElasticFactoryDomainRhs ( MoFEM::Interface & m_field,
boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > & pip,
std::string field_name,
std::string elastic_block_name,
std::string thermal_block_name,
std::string thermoelastic_block_name,
Sev sev,
double scale = 1 )
inlineprivate

Definition at line 230 of file thermo_elastic.cpp.

235 {
237
238 auto elastic_common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
239 m_field, pip, field_name, elastic_block_name, sev, scale);
240 auto thermal_common_ptr = boost::make_shared<BlockedThermalParameters>();
241 CHKERR addMatThermalBlockOps(pip, thermal_block_name, thermal_common_ptr,
242 Sev::inform);
243 auto thermoelastic_common_ptr =
244 boost::make_shared<BlockedThermoElasticParameters>();
245 CHKERR addMatThermoElasticBlockOps(pip, thermoelastic_block_name,
246 thermoelastic_common_ptr, Sev::inform);
248 m_field, pip, field_name, elastic_common_ptr, thermal_common_ptr,
249 thermoelastic_common_ptr, sev);
250
252 }

◆ runProblem()

MoFEMErrorCode ThermoElasticProblem::runProblem ( )

[Run problem]

Examples
thermo_elastic.cpp.

Definition at line 525 of file thermo_elastic.cpp.

525 {
529 CHKERR bC();
530 CHKERR OPs();
531 CHKERR tsSolve();
533}
MoFEMErrorCode getCommandLineParameters()
[Run problem]
MoFEMErrorCode setupProblem()
add fields
MoFEMErrorCode bC()
[Set up problem]
MoFEMErrorCode OPs()
[Boundary condition]
MoFEMErrorCode tsSolve()
[Push operators to pipeline]

◆ setupProblem()

MoFEMErrorCode ThermoElasticProblem::setupProblem ( )
private

add fields

[Get command line parameters]

[Set up problem]

Examples
thermo_elastic.cpp.

Definition at line 613 of file thermo_elastic.cpp.

613 {
616 // Add field
618 // Mechanical fields
620 CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE, SPACE_DIM);
621 // Temperature
622 constexpr auto flux_space = (SPACE_DIM == 2) ? HCURL : HDIV;
623 CHKERR simple->addDomainField("T", L2, base, 1);
624 CHKERR simple->addDomainField("FLUX", flux_space, base, 1);
625 CHKERR simple->addBoundaryField("FLUX", flux_space, base, 1);
626 CHKERR simple->addBoundaryField("TBC", L2, base, 1);
627
628 CHKERR simple->setFieldOrder("U", order_disp);
629 CHKERR simple->setFieldOrder("FLUX", order_flux);
630 CHKERR simple->setFieldOrder("T", order_temp);
631 CHKERR simple->setFieldOrder("TBC", order_temp);
632
633 CHKERR simple->setUp();
634
635 int coords_dim = 3;
636 CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
637 fieldEvalCoords.data(), &coords_dim,
638 &doEvalField);
639
640 tempFieldPtr = boost::make_shared<VectorDouble>();
641 fluxFieldPtr = boost::make_shared<MatrixDouble>();
642 dispFieldPtr = boost::make_shared<MatrixDouble>();
643 dispGradPtr = boost::make_shared<MatrixDouble>();
644 strainFieldPtr = boost::make_shared<MatrixDouble>();
645 stressFieldPtr = boost::make_shared<MatrixDouble>();
646
647 if (doEvalField) {
649 mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
650
651 CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree<SPACE_DIM>(
652 fieldEvalData, simple->getDomainFEName());
653
654 fieldEvalData->setEvalPoints(fieldEvalCoords.data(), 1);
655 auto no_rule = [](int, int, int) { return -1; };
656
657 auto field_eval_fe_ptr = fieldEvalData->feMethodPtr.lock();
658 field_eval_fe_ptr->getRuleHook = no_rule;
659
660 auto block_thermal_params = boost::make_shared<BlockedThermalParameters>();
661 auto block_thermoelastic_params =
662 boost::make_shared<BlockedThermoElasticParameters>();
663 auto coeff_expansion_ptr =
664 block_thermoelastic_params->getCoeffExpansionPtr();
665 auto ref_temp_ptr = block_thermoelastic_params->getRefTempPtr();
666
667 CHKERR addMatThermalBlockOps(field_eval_fe_ptr->getOpPtrVector(),
668 "MAT_THERMAL", block_thermal_params,
669 Sev::verbose);
671 field_eval_fe_ptr->getOpPtrVector(), "MAT_THERMOELASTIC",
672 block_thermoelastic_params, Sev::verbose);
673
675 field_eval_fe_ptr->getOpPtrVector(), {H1, HDIV});
676
677 auto hencky_common_data_ptr =
679 mField, field_eval_fe_ptr->getOpPtrVector(), "U", "MAT_ELASTIC",
680 Sev::inform);
681 auto mat_D_ptr = hencky_common_data_ptr->matDPtr;
682 auto dispGradPtr = hencky_common_data_ptr->matGradPtr;
683 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
684
685 field_eval_fe_ptr->getOpPtrVector().push_back(
687 field_eval_fe_ptr->getOpPtrVector().push_back(
689 field_eval_fe_ptr->getOpPtrVector().push_back(
691 field_eval_fe_ptr->getOpPtrVector().push_back(
693 dispGradPtr));
694
696
697 field_eval_fe_ptr->getOpPtrVector().push_back(
698 new
699 typename H::template OpCalculateHenckyThermalStress<SPACE_DIM, IT, 0>(
700 "U", tempFieldPtr, hencky_common_data_ptr, coeff_expansion_ptr,
701 ref_temp_ptr));
702 if (!IS_LARGE_STRAINS) {
703 field_eval_fe_ptr->getOpPtrVector().push_back(
705 hencky_common_data_ptr->getMatFirstPiolaStress(),
707 field_eval_fe_ptr->getOpPtrVector().push_back(
709 } else {
710 field_eval_fe_ptr->getOpPtrVector().push_back(
711 new typename H::template OpCalculateLogC<SPACE_DIM, IT>(
712 "U", hencky_common_data_ptr));
713 stressFieldPtr = hencky_common_data_ptr->getMatFirstPiolaStress();
714 strainFieldPtr = hencky_common_data_ptr->getMatLogC();
715 };
716 }
717
719}
FieldApproximationBase
approximation base
Definition definitions.h:58
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
@ L2
field with C-1 continuity
Definition definitions.h:88
@ HCURL
field with continuous tangents
Definition definitions.h:86
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Field evaluator interface.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Get values at integration pts for tensor field rank 1, i.e. vector field.
MoFEMErrorCode addDomainField(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition Simple.cpp:261
boost::shared_ptr< FieldEvaluatorInterface::SetPtsData > fieldEvalData
boost::shared_ptr< MatrixDouble > fluxFieldPtr
boost::shared_ptr< VectorDouble > tempFieldPtr
boost::shared_ptr< MatrixDouble > stressFieldPtr
std::array< double, 3 > fieldEvalCoords
boost::shared_ptr< MatrixDouble > dispGradPtr
boost::shared_ptr< MatrixDouble > dispFieldPtr
boost::shared_ptr< MatrixDouble > strainFieldPtr
constexpr bool IS_LARGE_STRAINS

◆ tsSolve()

MoFEMErrorCode ThermoElasticProblem::tsSolve ( )
private

[Push operators to pipeline]

[Solve]

Examples
thermo_elastic.cpp.

Definition at line 1284 of file thermo_elastic.cpp.

1284 {
1286
1289 ISManager *is_manager = mField.getInterface<ISManager>();
1290
1291 auto dm = simple->getDM();
1292 auto solver = pipeline_mng->createTSIM();
1293 auto snes_ctx_ptr = getDMSnesCtx(dm);
1294
1295 auto set_section_monitor = [&](auto solver) {
1297 SNES snes;
1298 CHKERR TSGetSNES(solver, &snes);
1299 CHKERR SNESMonitorSet(snes,
1300 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
1301 void *))MoFEMSNESMonitorFields,
1302 (void *)(snes_ctx_ptr.get()), nullptr);
1304 };
1305
1306 auto create_post_process_elements = [&]() {
1307 auto block_thermal_params = boost::make_shared<BlockedThermalParameters>();
1308
1309 auto block_thermoelastic_params =
1310 boost::make_shared<BlockedThermoElasticParameters>();
1311 auto coeff_expansion_ptr =
1312 block_thermoelastic_params->getCoeffExpansionPtr();
1313 auto ref_temp_ptr = block_thermoelastic_params->getRefTempPtr();
1314
1315 auto u_ptr = boost::make_shared<MatrixDouble>();
1316 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
1317 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
1318 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
1319 auto vec_temp_ptr = boost::make_shared<VectorDouble>();
1320 auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
1321
1322 auto push_domain_ops = [&](auto &pp_fe) {
1324 auto &pip = pp_fe->getOpPtrVector();
1325
1326 CHKERR addMatThermalBlockOps(pip, "MAT_THERMAL", block_thermal_params,
1327 Sev::verbose);
1329 pip, "MAT_THERMOELASTIC", block_thermoelastic_params, Sev::verbose);
1330
1332
1333 pip.push_back(new OpCalculateScalarFieldValues("T", vec_temp_ptr));
1334 pip.push_back(
1335 new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
1336
1337 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
1339 "U", mat_grad_ptr));
1340 auto elastic_common_ptr = commonDataFactory<SPACE_DIM, IT, DomainEleOp>(
1341 mField, pip, "U", "MAT_ELASTIC", Sev::inform);
1343 pip.push_back(
1344 new
1345 typename H::template OpCalculateHenckyThermalStress<SPACE_DIM, IT, 0>(
1346 "U", vec_temp_ptr, elastic_common_ptr, coeff_expansion_ptr,
1347 ref_temp_ptr));
1348 if (!IS_LARGE_STRAINS) {
1349 pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
1350 elastic_common_ptr->getMatFirstPiolaStress(), mat_stress_ptr));
1351 pip.push_back(
1352 new OpSymmetrizeTensor<SPACE_DIM>(mat_grad_ptr, mat_strain_ptr));
1353 } else {
1354 mat_stress_ptr = elastic_common_ptr->getMatFirstPiolaStress();
1355 mat_strain_ptr = elastic_common_ptr->getMatLogC();
1356 }
1357
1359 };
1360
1361 auto push_post_proc_ops = [&](auto &pp_fe) {
1363 auto &pip = pp_fe->getOpPtrVector();
1365
1366 if (!IS_LARGE_STRAINS) {
1367 pip.push_back(
1368
1369 new OpPPMap(
1370
1371 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
1372
1373 {{"T", vec_temp_ptr}},
1374
1375 {{"U", u_ptr}, {"FLUX", mat_flux_ptr}},
1376
1377 {},
1378
1379 {{"CAUCHY", mat_stress_ptr}, {"STRAIN", mat_strain_ptr}}
1380
1381 )
1382
1383 );
1384 } else {
1385 pip.push_back(
1386
1387 new OpPPMap(
1388
1389 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
1390
1391 {{"T", vec_temp_ptr}},
1392
1393 {{"U", u_ptr}, {"FLUX", mat_flux_ptr}},
1394
1395 {{"PIOLA", mat_stress_ptr}},
1396
1397 {{"HENCKY_STRAIN", mat_strain_ptr}}
1398
1399 )
1400
1401 );
1402 }
1403
1405 };
1406
1407 auto domain_post_proc = [&]() {
1408 if (do_output_domain == PETSC_FALSE)
1409 return boost::shared_ptr<PostProcEle>();
1410 auto pp_fe = boost::make_shared<PostProcEle>(mField);
1411 CHK_MOAB_THROW(push_domain_ops(pp_fe),
1412 "push domain ops to domain element");
1413 CHK_MOAB_THROW(push_post_proc_ops(pp_fe),
1414 "push post proc ops to domain element");
1415 return pp_fe;
1416 };
1417
1418 auto skin_post_proc = [&]() {
1419 if (do_output_skin == PETSC_FALSE)
1420 return boost::shared_ptr<SkinPostProcEle>();
1421 auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
1422 auto simple = mField.getInterface<Simple>();
1423 auto op_side = new OpLoopSide<SideEle>(mField, simple->getDomainFEName(),
1424 SPACE_DIM, Sev::verbose);
1425 CHK_MOAB_THROW(push_domain_ops(op_side),
1426 "push domain ops to side element");
1427 pp_fe->getOpPtrVector().push_back(op_side);
1428 CHK_MOAB_THROW(push_post_proc_ops(pp_fe),
1429 "push post proc ops to skin element");
1430 return pp_fe;
1431 };
1432
1433 return std::make_pair(domain_post_proc(), skin_post_proc());
1434 };
1435
1436 auto monitor_ptr = boost::make_shared<FEMethod>();
1437
1438 auto set_time_monitor = [&](auto dm, auto solver, auto domain_post_proc_fe,
1439 auto skin_post_proc_fe) {
1441 monitor_ptr->preProcessHook = [&]() {
1443
1444 if (save_every && (monitor_ptr->ts_step % save_every == 0)) {
1445 if (do_output_domain) {
1446 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
1447 domain_post_proc_fe,
1448 monitor_ptr->getCacheWeakPtr());
1449 CHKERR domain_post_proc_fe->writeFile(
1450 "out_" + boost::lexical_cast<std::string>(monitor_ptr->ts_step) +
1451 ".h5m");
1452 }
1453 if (do_output_skin) {
1454 CHKERR DMoFEMLoopFiniteElements(dm, simple->getBoundaryFEName(),
1455 skin_post_proc_fe,
1456 monitor_ptr->getCacheWeakPtr());
1457 CHKERR skin_post_proc_fe->writeFile(
1458 "out_skin_" +
1459 boost::lexical_cast<std::string>(monitor_ptr->ts_step) + ".h5m");
1460 }
1461 }
1462
1463 if (doEvalField) {
1464
1466 ->evalFEAtThePoint<SPACE_DIM>(
1467 fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
1468 simple->getDomainFEName(), fieldEvalData,
1470 MF_EXIST, QUIET);
1471
1472 if (atom_test) {
1473 auto eval_num_vec =
1474 createVectorMPI(mField.get_comm(), PETSC_DECIDE, 1);
1475 CHKERR VecZeroEntries(eval_num_vec);
1476 if (tempFieldPtr->size()) {
1477 CHKERR VecSetValue(eval_num_vec, 0, 1, ADD_VALUES);
1478 }
1479 CHKERR VecAssemblyBegin(eval_num_vec);
1480 CHKERR VecAssemblyEnd(eval_num_vec);
1481
1482 double eval_num;
1483 CHKERR VecSum(eval_num_vec, &eval_num);
1484 if (!(int)eval_num) {
1485 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1486 "atom test %d failed: did not find elements to evaluate "
1487 "the field, check the coordinates",
1488 atom_test);
1489 }
1490 }
1491
1492 if (tempFieldPtr->size()) {
1493 auto t_temp = getFTensor0FromVec(*tempFieldPtr);
1494 MOFEM_LOG("ThermoElasticSync", Sev::inform)
1495 << "Eval point T: " << t_temp;
1496 if (atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
1497 if (atom_test <= 3 && fabs(t_temp - 554.48) > 1e-2) {
1498 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1499 "atom test %d failed: wrong temperature value",
1500 atom_test);
1501 }
1502 if (atom_test == 4 && fabs(t_temp - 250) > 1e-2) {
1503 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1504 "atom test %d failed: wrong temperature value",
1505 atom_test);
1506 }
1507 if (atom_test == 5 && fabs(t_temp - 1) > 1e-2) {
1508 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1509 "atom test %d failed: wrong temperature value",
1510 atom_test);
1511 }
1512 }
1513 }
1514 if (fluxFieldPtr->size1()) {
1517 auto flux_mag = sqrt(t_flux(i) * t_flux(i));
1518 MOFEM_LOG("ThermoElasticSync", Sev::inform)
1519 << "Eval point FLUX magnitude: " << flux_mag;
1520 if (atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
1521 if (atom_test <= 3 && fabs(flux_mag - 27008.0) > 2e1) {
1522 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1523 "atom test %d failed: wrong flux value", atom_test);
1524 }
1525 if (atom_test == 4 && fabs(flux_mag - 150e3) > 1e-1) {
1526 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1527 "atom test %d failed: wrong flux value", atom_test);
1528 }
1529 if (atom_test == 5 && fabs(flux_mag) > 1e-6) {
1530 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1531 "atom test %d failed: wrong flux value", atom_test);
1532 }
1533 }
1534 }
1535 if (dispFieldPtr->size1()) {
1538 auto disp_mag = sqrt(t_disp(i) * t_disp(i));
1539 MOFEM_LOG("ThermoElasticSync", Sev::inform)
1540 << "Eval point U magnitude: " << disp_mag;
1541 if (atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
1542 if (atom_test == 1 && fabs(disp_mag - 0.00345) > 1e-5) {
1543 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1544 "atom test %d failed: wrong displacement value",
1545 atom_test);
1546 }
1547 if ((atom_test == 2 || atom_test == 3) &&
1548 fabs(disp_mag - 0.00265) > 1e-5) {
1549 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1550 "atom test %d failed: wrong displacement value",
1551 atom_test);
1552 }
1553 if ((atom_test == 5) &&
1554 fabs(t_disp(0) - (std::sqrt(std::exp(0.2)) - 1)) > 1e-5 &&
1555 fabs(t_disp(1) - (std::sqrt(std::exp(0.2)) - 1)) > 1e-5) {
1556 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1557 "atom test %d failed: wrong displacement value",
1558 atom_test);
1559 }
1560 }
1561 }
1562 if (strainFieldPtr->size1()) {
1564 auto t_strain =
1566 auto t_strain_trace = t_strain(i, i);
1567 if (atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
1568 if (atom_test == 1 && fabs(t_strain_trace - 0.00679) > 1e-5) {
1569 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1570 "atom test %d failed: wrong strain value", atom_test);
1571 }
1572 if ((atom_test == 2 || atom_test == 3) &&
1573 fabs(t_strain_trace - 0.00522) > 1e-5) {
1574 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1575 "atom test %d failed: wrong strain value", atom_test);
1576 }
1577 if ((atom_test == 5) && fabs(t_strain_trace - 0.2) > 1e-5) {
1578 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1579 "atom test %d failed: wrong strain value", atom_test);
1580 }
1581 }
1582 }
1583 if (stressFieldPtr->size1()) {
1584 double von_mises_stress;
1585 auto getVonMisesStress = [&](auto t_stress) {
1587 von_mises_stress = std::sqrt(
1588 0.5 *
1589 ((t_stress(0, 0) - t_stress(1, 1)) *
1590 (t_stress(0, 0) - t_stress(1, 1)) +
1591 (SPACE_DIM == 3 ? (t_stress(1, 1) - t_stress(2, 2)) *
1592 (t_stress(1, 1) - t_stress(2, 2))
1593 : 0) +
1594 (SPACE_DIM == 3 ? (t_stress(2, 2) - t_stress(0, 0)) *
1595 (t_stress(2, 2) - t_stress(0, 0))
1596 : 0) +
1597 6 * (t_stress(0, 1) * t_stress(0, 1) +
1598 (SPACE_DIM == 3 ? t_stress(1, 2) * t_stress(1, 2) : 0) +
1599 (SPACE_DIM == 3 ? t_stress(2, 0) * t_stress(2, 0) : 0))));
1601 };
1602 if (!IS_LARGE_STRAINS) {
1603 auto t_stress =
1605 CHKERR getVonMisesStress(t_stress);
1606 } else {
1607 auto t_stress =
1609 CHKERR getVonMisesStress(t_stress);
1610 }
1611 MOFEM_LOG("ThermoElasticSync", Sev::inform)
1612 << "Eval point von Mises Stress: " << von_mises_stress;
1613 if (atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
1614 if (atom_test == 1 && fabs(von_mises_stress - 523.0) > 5e-1) {
1615 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1616 "atom test %d failed: wrong von Mises stress value",
1617 atom_test);
1618 }
1619 if (atom_test == 2 && fabs(von_mises_stress - 16.3) > 5e-2) {
1620 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1621 "atom test %d failed: wrong von Mises stress value",
1622 atom_test);
1623 }
1624 if (atom_test == 3 && fabs(von_mises_stress - 14.9) > 5e-2) {
1625 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1626 "atom test %d failed: wrong von Mises stress value",
1627 atom_test);
1628 }
1629 if (atom_test == 5 && fabs(von_mises_stress) > 5e-2) {
1630 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1631 "atom test %d failed: wrong von Mises stress value",
1632 atom_test);
1633 }
1634 }
1635 }
1636
1638 }
1639
1641 };
1642 auto null = boost::shared_ptr<FEMethod>();
1643 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(), null,
1644 monitor_ptr, null);
1646 };
1647
1648 auto set_fieldsplit_preconditioner = [&](auto solver) {
1650
1651 SNES snes;
1652 CHKERR TSGetSNES(solver, &snes);
1653 KSP ksp;
1654 CHKERR SNESGetKSP(snes, &ksp);
1655 PC pc;
1656 CHKERR KSPGetPC(ksp, &pc);
1657 PetscBool is_pcfs = PETSC_FALSE;
1658 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1659
1660 // Setup fieldsplit (block) solver - optional: yes/no
1661 if (is_pcfs == PETSC_TRUE) {
1662 auto is_mng = mField.getInterface<ISManager>();
1663 auto name_prb = simple->getProblemName();
1664
1665 SmartPetscObj<IS> is_u;
1666 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "U", 0,
1667 SPACE_DIM, is_u);
1668 SmartPetscObj<IS> is_flux;
1669 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "FLUX", 0, 0,
1670 is_flux);
1671 SmartPetscObj<IS> is_T;
1672 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "T", 0, 0,
1673 is_T);
1674 SmartPetscObj<IS> is_TBC;
1675 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "TBC", 0, 0,
1676 is_TBC);
1677 IS is_tmp, is_tmp2;
1678 CHKERR ISExpand(is_T, is_flux, &is_tmp);
1679 CHKERR ISExpand(is_TBC, is_tmp, &is_tmp2);
1680 CHKERR ISDestroy(&is_tmp);
1681 auto is_TFlux = SmartPetscObj<IS>(is_tmp2);
1682
1683 CHKERR ISSort(is_u);
1684 CHKERR ISSort(is_TFlux);
1685 CHKERR PCFieldSplitSetIS(pc, PETSC_NULLPTR, is_TFlux);
1686 CHKERR PCFieldSplitSetIS(pc, PETSC_NULLPTR, is_u);
1687 }
1688
1690 };
1691
1692 auto B = createDMMatrix(dm);
1693 CHKERR TSSetIJacobian(solver, B, B, PETSC_NULLPTR, PETSC_NULLPTR);
1694 auto D = createDMVector(dm);
1695 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
1696 CHKERR TSSetSolution(solver, D);
1697 CHKERR TSSetFromOptions(solver);
1698
1699 CHKERR set_section_monitor(solver);
1700 CHKERR set_fieldsplit_preconditioner(solver);
1701
1702 auto [domain_post_proc_fe, skin_post_proc_fe] =
1703 create_post_process_elements();
1704 CHKERR set_time_monitor(dm, solver, domain_post_proc_fe, skin_post_proc_fe);
1705
1706 CHKERR TSSetUp(solver);
1707 CHKERR TSSolve(solver, NULL);
1708
1710}
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
@ QUIET
@ ROW
@ MF_EXIST
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition DMMoFEM.cpp:514
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:576
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1102
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:1059
double D
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:1046
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *ctx)
Sens monitor printing residual field by field.
Definition SnesCtx.cpp:592
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
static auto getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition DMMoFEM.hpp:1130
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Section manager is used to create indexes and sections.
Definition ISManager.hpp:23
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
intrusive_ptr for managing petsc objects

Member Data Documentation

◆ dispFieldPtr

boost::shared_ptr<MatrixDouble> ThermoElasticProblem::dispFieldPtr
private
Examples
thermo_elastic.cpp.

Definition at line 145 of file thermo_elastic.cpp.

◆ dispGradPtr

boost::shared_ptr<MatrixDouble> ThermoElasticProblem::dispGradPtr
private
Examples
thermo_elastic.cpp.

Definition at line 146 of file thermo_elastic.cpp.

◆ doEvalField

PetscBool ThermoElasticProblem::doEvalField
private
Examples
thermo_elastic.cpp.

Definition at line 139 of file thermo_elastic.cpp.

◆ fieldEvalCoords

std::array<double, 3> ThermoElasticProblem::fieldEvalCoords = {0.0, 0.0, 0.0}
private
Examples
thermo_elastic.cpp.

Definition at line 140 of file thermo_elastic.cpp.

140{0.0, 0.0, 0.0};

◆ fieldEvalData

boost::shared_ptr<FieldEvaluatorInterface::SetPtsData> ThermoElasticProblem::fieldEvalData
private
Examples
thermo_elastic.cpp.

Definition at line 141 of file thermo_elastic.cpp.

◆ fluxFieldPtr

boost::shared_ptr<MatrixDouble> ThermoElasticProblem::fluxFieldPtr
private
Examples
thermo_elastic.cpp.

Definition at line 144 of file thermo_elastic.cpp.

◆ mField

MoFEM::Interface& ThermoElasticProblem::mField
private
Examples
thermo_elastic.cpp.

Definition at line 137 of file thermo_elastic.cpp.

◆ strainFieldPtr

boost::shared_ptr<MatrixDouble> ThermoElasticProblem::strainFieldPtr
private
Examples
thermo_elastic.cpp.

Definition at line 147 of file thermo_elastic.cpp.

◆ stressFieldPtr

boost::shared_ptr<MatrixDouble> ThermoElasticProblem::stressFieldPtr
private
Examples
thermo_elastic.cpp.

Definition at line 148 of file thermo_elastic.cpp.

◆ tempFieldPtr

boost::shared_ptr<VectorDouble> ThermoElasticProblem::tempFieldPtr
private
Examples
thermo_elastic.cpp.

Definition at line 143 of file thermo_elastic.cpp.


The documentation for this struct was generated from the following file: