9#ifndef EXECUTABLE_DIMENSION
10 #define EXECUTABLE_DIMENSION 3
13#ifndef FINITE_DEFORMATION_FLAG
14 #define FINITE_DEFORMATION_FLAG true
46 IntegrationType::GAUSS;
83 GAUSS>::OpEssentialRhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
86 GAUSS>::OpEssentialLhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
121auto save_range = [](moab::Interface &moab,
const std::string name,
125 CHKERR moab.add_entities(*out_meshset, r);
126 CHKERR moab.write_file(name.c_str(),
"VTK",
"", out_meshset->get_ptr(), 1);
158 :
public boost::enable_shared_from_this<BlockedThermalParameters> {
167 return boost::shared_ptr<double>(shared_from_this(), &
heatCapacity);
172 :
public boost::enable_shared_from_this<BlockedThermoElasticParameters> {
177 return boost::shared_ptr<VectorDouble>(shared_from_this(),
182 return boost::shared_ptr<double>(shared_from_this(), &
refTemp);
187 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
188 std::string block_name,
189 boost::shared_ptr<BlockedThermalParameters> blockedParamsPtr,
Sev sev);
192 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
193 std::string block_name,
194 boost::shared_ptr<BlockedThermoElasticParameters> blockedParamsPtr,
197 template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
200 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
202 boost::shared_ptr<HenckyOps::CommonData> elastic_common_ptr,
203 boost::shared_ptr<ThermoElasticProblem::BlockedThermalParameters>
205 boost::shared_ptr<ThermoElasticProblem::BlockedThermoElasticParameters>
206 thermoelastic_common_ptr,
213 auto vec_temp_ptr = boost::make_shared<VectorDouble>();
215 auto coeff_expansion_ptr = thermoelastic_common_ptr->getCoeffExpansionPtr();
216 auto ref_temp_ptr = thermoelastic_common_ptr->getRefTempPtr();
218 new typename H::template OpCalculateHenckyThermalStress<DIM, I, 0>(
219 "U", vec_temp_ptr, elastic_common_ptr, coeff_expansion_ptr,
224 "U", elastic_common_ptr->getMatFirstPiolaStress()));
229 template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
232 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
233 std::string
field_name, std::string elastic_block_name,
234 std::string thermal_block_name, std::string thermoelastic_block_name,
240 auto thermal_common_ptr = boost::make_shared<BlockedThermalParameters>();
243 auto thermoelastic_common_ptr =
244 boost::make_shared<BlockedThermoElasticParameters>();
246 thermoelastic_common_ptr, Sev::inform);
248 m_field, pip,
field_name, elastic_common_ptr, thermal_common_ptr,
249 thermoelastic_common_ptr, sev);
254 template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
257 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
258 std::string
field_name, std::string coupled_field_name,
259 boost::shared_ptr<HenckyOps::CommonData> elastic_common_ptr,
260 boost::shared_ptr<ThermoElasticProblem::BlockedThermalParameters>
262 boost::shared_ptr<ThermoElasticProblem::BlockedThermoElasticParameters>
263 thermoelastic_common_ptr,
269 using OpKPiola =
typename B::template OpGradTensorGrad<1, DIM, DIM, 1>;
272 auto vec_temp_ptr = boost::make_shared<VectorDouble>();
274 auto coeff_expansion_ptr = thermoelastic_common_ptr->getCoeffExpansionPtr();
275 auto ref_temp_ptr = thermoelastic_common_ptr->getRefTempPtr();
277 new typename H::template OpCalculateHenckyThermalStress<DIM, I, 0>(
278 "U", vec_temp_ptr, elastic_common_ptr, coeff_expansion_ptr,
280 pip.push_back(
new typename H::template OpHenckyTangent<DIM, I, 0>(
283 elastic_common_ptr->getMatTangent()));
284 pip.push_back(
new typename H::template OpCalculateHenckyThermalStressdT<
286 field_name, coupled_field_name, elastic_common_ptr,
287 coeff_expansion_ptr));
292 template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
295 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
296 std::string
field_name, std::string coupled_field_name,
297 std::string elastic_block_name, std::string thermal_block_name,
298 std::string thermoelastic_block_name,
Sev sev,
double scale = 1) {
303 auto thermal_common_ptr = boost::make_shared<BlockedThermalParameters>();
306 auto thermoelastic_common_ptr =
307 boost::make_shared<BlockedThermoElasticParameters>();
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);
319 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
320 std::string block_name,
321 boost::shared_ptr<BlockedThermalParameters> blockedParamsPtr,
Sev sev) {
325 OpMatThermalBlocks(boost::shared_ptr<double> conductivity_ptr,
326 boost::shared_ptr<double> capacity_ptr,
328 std::vector<const CubitMeshSets *> meshset_vec_ptr)
330 conductivityPtr(conductivity_ptr), capacityPtr(capacity_ptr) {
332 "Cannot get data from thermal block");
339 for (
auto &b : blockData) {
341 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
342 *conductivityPtr = b.conductivity;
343 *capacityPtr = b.capacity;
361 std::vector<BlockData> blockData;
365 std::vector<const CubitMeshSets *> meshset_vec_ptr,
369 for (
auto m : meshset_vec_ptr) {
371 std::vector<double> block_data;
372 CHKERR m->getAttributes(block_data);
373 if (block_data.size() < 2) {
375 "Expected that block has at least two attributes");
377 auto get_block_ents = [&]() {
380 m_field.
get_moab().get_entities_by_handle(
m->meshset, ents,
true);
385 <<
m->getName() <<
": conductivity = " << block_data[0]
386 <<
" capacity = " << block_data[1];
388 blockData.push_back({block_data[0], block_data[1], get_block_ents()});
394 boost::shared_ptr<double> conductivityPtr;
395 boost::shared_ptr<double> capacityPtr;
398 pipeline.push_back(
new OpMatThermalBlocks(
399 blockedParamsPtr->getHeatConductivityPtr(),
400 blockedParamsPtr->getHeatCapacityPtr(),
mField, sev,
405 (boost::format(
"%s(.*)") % block_name).str()
415 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
416 std::string block_name,
417 boost::shared_ptr<BlockedThermoElasticParameters> blockedParamsPtr,
421 struct OpMatThermoElasticBlocks :
public DomainEleOp {
422 OpMatThermoElasticBlocks(boost::shared_ptr<VectorDouble> expansion_ptr,
423 boost::shared_ptr<double> ref_temp_ptr,
425 std::vector<const CubitMeshSets *> meshset_vec_ptr)
427 expansionPtr(expansion_ptr), refTempPtr(ref_temp_ptr) {
429 extractThermoElasticBlockData(m_field, meshset_vec_ptr, sev),
430 "Cannot get data from thermoelastic block");
437 for (
auto &b : blockData) {
439 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
440 *expansionPtr = b.expansion;
441 *refTempPtr = b.ref_temp;
459 std::vector<BlockData> blockData;
463 std::vector<const CubitMeshSets *> meshset_vec_ptr,
Sev sev) {
466 for (
auto m : meshset_vec_ptr) {
468 std::vector<double> block_data;
469 CHKERR m->getAttributes(block_data);
470 if (block_data.size() < 2) {
472 "Expected that block has at least two attributes");
474 auto get_block_ents = [&]() {
477 m_field.
get_moab().get_entities_by_handle(
m->meshset, ents,
true);
481 auto get_expansion = [&]() {
483 if (block_data.size() > 2) {
484 expansion[1] = block_data[2];
486 if (
SPACE_DIM == 3 && block_data.size() > 3) {
487 expansion[2] = block_data[3];
492 auto coeff_exp_vec = get_expansion();
495 <<
" ref_temp = " << block_data[0]
496 <<
" expansion = " << coeff_exp_vec;
498 blockData.push_back({block_data[0], coeff_exp_vec, get_block_ents()});
504 boost::shared_ptr<VectorDouble> expansionPtr;
505 boost::shared_ptr<double> refTempPtr;
508 pipeline.push_back(
new OpMatThermoElasticBlocks(
509 blockedParamsPtr->getCoeffExpansionPtr(),
510 blockedParamsPtr->getRefTempPtr(),
mField, sev,
515 (boost::format(
"%s(.*)") % block_name).str()
540 auto get_command_line_parameters = [&]() {
606 CHKERR get_command_line_parameters();
624 CHKERR simple->addDomainField(
"FLUX", flux_space, base, 1);
625 CHKERR simple->addBoundaryField(
"FLUX", flux_space, base, 1);
655 auto no_rule = [](int, int, int) {
return -1; };
658 field_eval_fe_ptr->getRuleHook = no_rule;
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();
668 "MAT_THERMAL", block_thermal_params,
671 field_eval_fe_ptr->getOpPtrVector(),
"MAT_THERMOELASTIC",
672 block_thermoelastic_params, Sev::verbose);
675 field_eval_fe_ptr->getOpPtrVector(), {H1, HDIV});
677 auto hencky_common_data_ptr =
679 mField, field_eval_fe_ptr->getOpPtrVector(),
"U",
"MAT_ELASTIC",
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>();
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(
697 field_eval_fe_ptr->getOpPtrVector().push_back(
699 typename H::template OpCalculateHenckyThermalStress<SPACE_DIM, IT, 0>(
700 "U",
tempFieldPtr, hencky_common_data_ptr, coeff_expansion_ptr,
703 field_eval_fe_ptr->getOpPtrVector().push_back(
705 hencky_common_data_ptr->getMatFirstPiolaStress(),
707 field_eval_fe_ptr->getOpPtrVector().push_back(
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();
732 auto get_skin = [&]() {
737 CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
741 auto filter_flux_blocks = [&](
auto skin,
bool temp_bc =
false) {
742 auto remove_cubit_blocks = [&](
auto c) {
752 skin = subtract(skin, ents);
757 auto remove_named_blocks = [&](
auto n) {
762 (boost::format(
"%s(.*)") %
n).str()
770 skin = subtract(skin, ents);
776 "remove_cubit_blocks");
778 "remove_named_blocks");
781 "remove_cubit_blocks");
790 ParallelComm *pcomm =
792 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
793 PSTATUS_NOT, -1, &boundary_ents);
794 return boundary_ents;
798 auto remove_temp_bc_ents =
804 remove_temp_bc_ents);
806 MOFEM_LOG(
"SYNC", Sev::noisy) << remove_flux_ents << endl;
809 MOFEM_LOG(
"SYNC", Sev::noisy) << remove_temp_bc_ents << endl;
822 remove_temp_bc_ents);
827 simple->getProblemName(),
"FLUX", remove_flux_ents);
829 simple->getProblemName(),
"TBC", remove_temp_bc_ents);
831 auto set_init_temp = [](boost::shared_ptr<FieldEntity> field_entity_ptr) {
839 simple->getProblemName(),
"U");
841 simple->getProblemName(),
"FLUX",
false);
857 auto boundary_marker =
858 bc_mng->getMergedBlocksMarker(vector<string>{
"HEATFLUX"});
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();
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();
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;
883 auto def_file_name = [time_scale](
const std::string &&name) {
884 return (!time_scale->argFileScale) ? name :
"";
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);
901 auto add_domain_rhs_ops = [&](
auto &pipeline) {
906 block_thermoelastic_params, Sev::inform);
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>();
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>();
926 "FLUX", vec_temp_div_ptr));
931 mField, pipeline,
"U",
"MAT_ELASTIC",
"MAT_THERMAL",
932 "MAT_THERMOELASTIC", Sev::inform);
934 auto resistance = [heat_conductivity_ptr](
const double,
const double,
936 return (1. / (*heat_conductivity_ptr));
939 auto capacity = [heat_capacity_ptr](
const double,
const double,
941 return -(*heat_capacity_ptr);
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));
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);
965 auto add_domain_lhs_ops = [&](
auto &pipeline) {
970 block_thermoelastic_params,
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;
980 auto resistance = [heat_conductivity_ptr](
const double,
const double,
982 return (1. / (*heat_conductivity_ptr));
984 auto capacity = [heat_capacity_ptr](
const double,
const double,
986 return -(*heat_capacity_ptr);
989 new OpHdivHdiv(
"FLUX",
"FLUX", resistance, mat_grad_ptr));
991 new OpHdivT(
"FLUX",
"T", []()
constexpr {
return -1; },
true));
993 auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
998 new OpHdivU(
"FLUX",
"U", mat_flux_ptr, resistance, mat_grad_ptr));
1001 mField, pipeline,
"U",
"T",
"MAT_ELASTIC",
"MAT_THERMAL",
1002 "MAT_THERMOELASTIC", Sev::inform);
1004 auto op_capacity =
new OpCapacity(
"T",
"T", capacity);
1005 op_capacity->feScalingFun = [](
const FEMethod *fe_ptr) {
1006 return fe_ptr->ts_a;
1008 pipeline.push_back(op_capacity);
1010 auto vec_temp_ptr = boost::make_shared<VectorDouble>();
1013 pipeline,
mField,
"T", vec_temp_ptr,
"SETTEMP", Sev::verbose);
1018 auto add_boundary_rhs_ops = [&](
auto &pipeline) {
1024 pipeline,
mField,
"U", {time_scale, time_force_scaling},
"FORCE",
1032 pipeline,
mField,
"FLUX", {time_scale, time_temperature_scaling},
1033 "TEMPERATURE", Sev::inform);
1043 pipeline,
mField,
"TBC", {time_scale, time_flux_scaling},
"FLUX",
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>();
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);
1058 auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
1067 struct OpTBCTimesNormalFlux
1072 OpTBCTimesNormalFlux(
const std::string
field_name,
1073 boost::shared_ptr<MatrixDouble> vec,
1074 boost::shared_ptr<Range> ents_ptr =
nullptr)
1081 auto t_w = OP::getFTensor0IntegrationWeight();
1085 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1089 for (
int gg = 0; gg != OP::nbIntegrationPts; gg++) {
1091 const double alpha = t_w * (t_vec(
i) * t_normal(
i));
1094 for (; rr != OP::nbRows; ++rr) {
1095 OP::locF[rr] += alpha * t_row_base;
1098 for (; rr < OP::nbRowBaseFunctions; ++rr)
1104 EntityType fe_type = OP::getNumeredEntFiniteElementPtr()->getEntType();
1105 if (fe_type == MBTRI) {
1112 boost::shared_ptr<MatrixDouble> sourceVec;
1114 pipeline.push_back(
new OpTBCTimesNormalFlux(
"TBC", mat_flux_ptr));
1116 struct OpBaseNormalTimesTBC
1121 OpBaseNormalTimesTBC(
const std::string
field_name,
1122 boost::shared_ptr<VectorDouble> vec,
1123 boost::shared_ptr<Range> ents_ptr =
nullptr)
1130 auto t_w = OP::getFTensor0IntegrationWeight();
1134 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1138 for (
int gg = 0; gg != OP::nbIntegrationPts; gg++) {
1140 const double alpha = t_w * t_vec;
1143 for (; rr != OP::nbRows; ++rr) {
1144 OP::locF[rr] += alpha * (t_row_base(
i) * t_normal(
i));
1147 for (; rr < OP::nbRowBaseFunctions; ++rr)
1153 EntityType fe_type = OP::getNumeredEntFiniteElementPtr()->getEntType();
1154 if (fe_type == MBTRI) {
1161 boost::shared_ptr<VectorDouble> sourceVec;
1164 pipeline.push_back(
new OpBaseNormalTimesTBC(
"FLUX", temp_bc_ptr));
1169 auto add_boundary_lhs_ops = [&](
auto &pipeline) {
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>();
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);
1192 struct OpTBCTimesNormalFlux
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) {
1202 this->assembleTranspose =
true;
1203 this->onlyTranspose =
false;
1213 auto t_w = OP::getFTensor0IntegrationWeight();
1217 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1219 for (
int gg = 0; gg != OP::nbIntegrationPts; gg++) {
1221 auto a_mat_ptr = &*OP::locMat.data().begin();
1223 for (; rr != OP::nbRows; rr++) {
1225 const double alpha = t_w * t_row_base;
1229 for (
int cc = 0; cc != OP::nbCols; cc++) {
1233 *a_mat_ptr += alpha * (t_col_base(
i) * t_normal(
i));
1239 for (; rr < OP::nbRowBaseFunctions; ++rr)
1244 EntityType fe_type = OP::getNumeredEntFiniteElementPtr()->getEntType();
1245 if (fe_type == MBTRI) {
1252 pipeline.push_back(
new OpTBCTimesNormalFlux(
"TBC",
"FLUX"));
1258 auto get_bc_hook_rhs = [&]() {
1260 mField, pipeline_mng->getDomainRhsFE(),
1261 {time_scale, time_displacement_scaling});
1264 auto get_bc_hook_lhs = [&]() {
1266 mField, pipeline_mng->getDomainLhsFE(),
1267 {time_scale, time_displacement_scaling});
1271 pipeline_mng->getDomainRhsFE()->preProcessHook = get_bc_hook_rhs();
1272 pipeline_mng->getDomainLhsFE()->preProcessHook = get_bc_hook_lhs();
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());
1291 auto dm =
simple->getDM();
1292 auto solver = pipeline_mng->createTSIM();
1295 auto set_section_monitor = [&](
auto solver) {
1298 CHKERR TSGetSNES(solver, &snes);
1299 CHKERR SNESMonitorSet(snes,
1302 (
void *)(snes_ctx_ptr.get()),
nullptr);
1306 auto create_post_process_elements = [&]() {
1307 auto block_thermal_params = boost::make_shared<BlockedThermalParameters>();
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();
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>();
1322 auto push_domain_ops = [&](
auto &pp_fe) {
1324 auto &pip = pp_fe->getOpPtrVector();
1329 pip,
"MAT_THERMOELASTIC", block_thermoelastic_params, Sev::verbose);
1339 "U", mat_grad_ptr));
1341 mField, pip,
"U",
"MAT_ELASTIC", Sev::inform);
1345 typename H::template OpCalculateHenckyThermalStress<SPACE_DIM, IT, 0>(
1346 "U", vec_temp_ptr, elastic_common_ptr, coeff_expansion_ptr,
1350 elastic_common_ptr->getMatFirstPiolaStress(), mat_stress_ptr));
1354 mat_stress_ptr = elastic_common_ptr->getMatFirstPiolaStress();
1355 mat_strain_ptr = elastic_common_ptr->getMatLogC();
1361 auto push_post_proc_ops = [&](
auto &pp_fe) {
1363 auto &pip = pp_fe->getOpPtrVector();
1371 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
1373 {{
"T", vec_temp_ptr}},
1375 {{
"U", u_ptr}, {
"FLUX", mat_flux_ptr}},
1379 {{
"CAUCHY", mat_stress_ptr}, {
"STRAIN", mat_strain_ptr}}
1389 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
1391 {{
"T", vec_temp_ptr}},
1393 {{
"U", u_ptr}, {
"FLUX", mat_flux_ptr}},
1395 {{
"PIOLA", mat_stress_ptr}},
1397 {{
"HENCKY_STRAIN", mat_strain_ptr}}
1407 auto domain_post_proc = [&]() {
1409 return boost::shared_ptr<PostProcEle>();
1410 auto pp_fe = boost::make_shared<PostProcEle>(mField);
1412 "push domain ops to domain element");
1414 "push post proc ops to domain element");
1418 auto skin_post_proc = [&]() {
1420 return boost::shared_ptr<SkinPostProcEle>();
1421 auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
1426 "push domain ops to side element");
1427 pp_fe->getOpPtrVector().push_back(op_side);
1429 "push post proc ops to skin element");
1433 return std::make_pair(domain_post_proc(), skin_post_proc());
1436 auto monitor_ptr = boost::make_shared<FEMethod>();
1438 auto set_time_monitor = [&](
auto dm,
auto solver,
auto domain_post_proc_fe,
1439 auto skin_post_proc_fe) {
1441 monitor_ptr->preProcessHook = [&]() {
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) +
1456 monitor_ptr->getCacheWeakPtr());
1457 CHKERR skin_post_proc_fe->writeFile(
1459 boost::lexical_cast<std::string>(monitor_ptr->ts_step) +
".h5m");
1466 ->evalFEAtThePoint<SPACE_DIM>(
1467 fieldEvalCoords.data(), 1e-12,
simple->getProblemName(),
1468 simple->getDomainFEName(), fieldEvalData,
1469 mField.get_comm_rank(), mField.get_comm_rank(),
nullptr,
1475 CHKERR VecZeroEntries(eval_num_vec);
1476 if (tempFieldPtr->size()) {
1477 CHKERR VecSetValue(eval_num_vec, 0, 1, ADD_VALUES);
1479 CHKERR VecAssemblyBegin(eval_num_vec);
1480 CHKERR VecAssemblyEnd(eval_num_vec);
1483 CHKERR VecSum(eval_num_vec, &eval_num);
1484 if (!(
int)eval_num) {
1486 "atom test %d failed: did not find elements to evaluate "
1487 "the field, check the coordinates",
1492 if (tempFieldPtr->size()) {
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) {
1499 "atom test %d failed: wrong temperature value",
1502 if (
atom_test == 4 && fabs(t_temp - 250) > 1e-2) {
1504 "atom test %d failed: wrong temperature value",
1507 if (
atom_test == 5 && fabs(t_temp - 1) > 1e-2) {
1509 "atom test %d failed: wrong temperature value",
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) {
1523 "atom test %d failed: wrong flux value",
atom_test);
1525 if (
atom_test == 4 && fabs(flux_mag - 150e3) > 1e-1) {
1527 "atom test %d failed: wrong flux value",
atom_test);
1529 if (
atom_test == 5 && fabs(flux_mag) > 1e-6) {
1531 "atom test %d failed: wrong flux value",
atom_test);
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) {
1544 "atom test %d failed: wrong displacement value",
1548 fabs(disp_mag - 0.00265) > 1e-5) {
1550 "atom test %d failed: wrong displacement value",
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) {
1557 "atom test %d failed: wrong displacement value",
1562 if (strainFieldPtr->size1()) {
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) {
1570 "atom test %d failed: wrong strain value",
atom_test);
1573 fabs(t_strain_trace - 0.00522) > 1e-5) {
1575 "atom test %d failed: wrong strain value",
atom_test);
1577 if ((
atom_test == 5) && fabs(t_strain_trace - 0.2) > 1e-5) {
1579 "atom test %d failed: wrong strain value",
atom_test);
1583 if (stressFieldPtr->size1()) {
1584 double von_mises_stress;
1585 auto getVonMisesStress = [&](
auto t_stress) {
1587 von_mises_stress = std::sqrt(
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))
1594 (
SPACE_DIM == 3 ? (t_stress(2, 2) - t_stress(0, 0)) *
1595 (t_stress(2, 2) - t_stress(0, 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))));
1605 CHKERR getVonMisesStress(t_stress);
1609 CHKERR getVonMisesStress(t_stress);
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) {
1616 "atom test %d failed: wrong von Mises stress value",
1619 if (
atom_test == 2 && fabs(von_mises_stress - 16.3) > 5e-2) {
1621 "atom test %d failed: wrong von Mises stress value",
1624 if (
atom_test == 3 && fabs(von_mises_stress - 14.9) > 5e-2) {
1626 "atom test %d failed: wrong von Mises stress value",
1629 if (
atom_test == 5 && fabs(von_mises_stress) > 5e-2) {
1631 "atom test %d failed: wrong von Mises stress value",
1642 auto null = boost::shared_ptr<FEMethod>();
1648 auto set_fieldsplit_preconditioner = [&](
auto solver) {
1652 CHKERR TSGetSNES(solver, &snes);
1654 CHKERR SNESGetKSP(snes, &ksp);
1656 CHKERR KSPGetPC(ksp, &pc);
1657 PetscBool is_pcfs = PETSC_FALSE;
1658 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1661 if (is_pcfs == PETSC_TRUE) {
1662 auto is_mng = mField.getInterface<
ISManager>();
1663 auto name_prb =
simple->getProblemName();
1666 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"U", 0,
1669 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"FLUX", 0, 0,
1672 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"T", 0, 0,
1675 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"TBC", 0, 0,
1678 CHKERR ISExpand(is_T, is_flux, &is_tmp);
1679 CHKERR ISExpand(is_TBC, is_tmp, &is_tmp2);
1680 CHKERR ISDestroy(&is_tmp);
1685 CHKERR PCFieldSplitSetIS(pc, PETSC_NULLPTR, is_TFlux);
1686 CHKERR PCFieldSplitSetIS(pc, PETSC_NULLPTR, is_u);
1693 CHKERR TSSetIJacobian(solver,
B,
B, PETSC_NULLPTR, PETSC_NULLPTR);
1696 CHKERR TSSetSolution(solver,
D);
1697 CHKERR TSSetFromOptions(solver);
1699 CHKERR set_section_monitor(solver);
1700 CHKERR set_fieldsplit_preconditioner(solver);
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);
1707 CHKERR TSSolve(solver, NULL);
1718 const char param_file[] =
"param_file.petsc";
1722 auto core_log = logging::core::get();
1736 DMType dm_name =
"DMMOFEM";
1741 moab::Core mb_instance;
1742 moab::Interface &moab = mb_instance;
static auto filter_true_skin(MoFEM::Interface &m_field, Range &&skin)
Operators for Piola transformation of the thermal conductivity.
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define FTENSOR_INDEX(DIM, I)
Implementation of thermal convection bc.
Thermal operators agnostic to small/large deformations.
Implementation of thermal radiation bc.
void simple(double P1[], double P2[], double P3[], double c[], const int N)
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
[Define dimension]
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
#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()
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#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
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
auto createDMVector(DM dm)
Get smart vector from DM.
auto createDMMatrix(DM dm)
Get smart matrix from DM.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
SeverityLevel
Severity levels.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#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
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double n
refractive index of diffusive medium
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)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasVector< double > VectorDouble
implementation of Data Operators for Forces and Sources
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.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *ctx)
Sens monitor printing residual field by field.
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)
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
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.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
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)
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
constexpr IntegrationType I
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
ElementsAndOps< SPACE_DIM >::SideEle SideEle
constexpr auto field_name
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
FTensor::Index< 'm', 3 > m
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Definition of the displacement bc data structure.
Data on single entity (This is passed as argument to DataOperator::doWork)
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.
structure for User Loop Methods on finite elements
Field evaluator interface.
Definition of the heat flux bc data structure.
Section manager is used to create indexes and sections.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
Interface for managing meshsets containing materials and boundary conditions.
Natural boundary conditions.
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 field rank 0, i.e. vector field.
Get values at integration pts for tensor field rank 1, i.e. vector field.
Enforce essential constrains on lhs.
Enforce essential constrains on rhs.
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
PipelineManager interface.
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
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.
MoFEMErrorCode getOptions()
get options
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
auto getHeatCapacityPtr()
auto getHeatConductivityPtr()
VectorDouble coeffExpansion
auto getCoeffExpansionPtr()
boost::shared_ptr< FieldEvaluatorInterface::SetPtsData > fieldEvalData
MoFEMErrorCode getCommandLineParameters()
[Run problem]
MoFEMErrorCode runProblem()
[Run problem]
boost::shared_ptr< MatrixDouble > fluxFieldPtr
boost::shared_ptr< VectorDouble > tempFieldPtr
boost::shared_ptr< MatrixDouble > stressFieldPtr
MoFEM::Interface & mField
ThermoElasticProblem(MoFEM::Interface &m_field)
std::array< double, 3 > fieldEvalCoords
boost::shared_ptr< MatrixDouble > dispGradPtr
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 setupProblem()
add fields
boost::shared_ptr< MatrixDouble > dispFieldPtr
MoFEMErrorCode bC()
[Set up problem]
MoFEMErrorCode OPs()
[Boundary condition]
boost::shared_ptr< MatrixDouble > strainFieldPtr
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
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)
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)
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)
#define FINITE_DEFORMATION_FLAG
constexpr AssemblyType AT
constexpr IntegrationType IT
static char help[]
[Solve]
constexpr bool IS_LARGE_STRAINS
#define EXECUTABLE_DIMENSION
double default_heat_capacity
double default_young_modulus
[Essential boundary conditions (Least square approach)]
double default_coeff_expansion
PetscBool do_output_domain
double default_heat_conductivity
double default_poisson_ratio