848 {
850
853
856
857 auto boundary_marker =
858 bc_mng->getMergedBlocksMarker(vector<string>{"HEATFLUX"});
859
862 };
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
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
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) {
904 Sev::inform);
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
923 pipeline.push_back(
926 "FLUX", vec_temp_div_ptr));
927 pipeline.push_back(
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,
936 return (1. / (*heat_conductivity_ptr));
937 };
938
939 auto capacity = [heat_capacity_ptr](
const double,
const double,
941 return -(*heat_capacity_ptr);
942 };
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
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) {
968 Sev::verbose);
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,
982 return (1. / (*heat_conductivity_ptr));
983 };
984 auto capacity = [heat_capacity_ptr](
const double,
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(
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>();
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
1028
1032 pipeline,
mField,
"FLUX", {time_scale, time_temperature_scaling},
1033 "TEMPERATURE", Sev::inform);
1034
1035
1036
1037
1038
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>();
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(
1061
1062
1063
1064
1065
1066
1067 struct OpTBCTimesNormalFlux
1069
1071
1072 OpTBCTimesNormalFlux(
const std::string
field_name,
1073 boost::shared_ptr<MatrixDouble> vec,
1074 boost::shared_ptr<Range> ents_ptr = nullptr)
1076
1080
1081 auto t_w = OP::getFTensor0IntegrationWeight();
1082
1084
1085 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1086
1088
1089 for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
1090
1091 const double alpha = t_w * (t_vec(
i) * t_normal(
i));
1092
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;
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
1118
1120
1121 OpBaseNormalTimesTBC(
const std::string
field_name,
1122 boost::shared_ptr<VectorDouble> vec,
1123 boost::shared_ptr<Range> ents_ptr = nullptr)
1125
1129
1130 auto t_w = OP::getFTensor0IntegrationWeight();
1131
1133
1134 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1135
1137
1138 for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
1139
1140 const double alpha = t_w * t_vec;
1141
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;
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 =
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);
1186
1187
1188
1189
1190
1191
1192 struct OpTBCTimesNormalFlux
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;
1204 }
1205
1209
1211
1212
1213 auto t_w = OP::getFTensor0IntegrationWeight();
1214
1216
1217 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1218
1219 for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
1220
1221 auto a_mat_ptr = &*OP::locMat.data().begin();
1222 int rr = 0;
1223 for (; rr != OP::nbRows; rr++) {
1224
1225 const double alpha = t_w * t_row_base;
1226
1228
1229 for (int cc = 0; cc != OP::nbCols; cc++) {
1230
1231
1232
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;
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
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)
@ HDIV
field with continuous normal traction
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
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.
structure for User Loop Methods on finite elements
Natural boundary conditions.
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)