#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 3
#endif
IntegrationType::GAUSS;
GAUSS>::OpEssentialRhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
GAUSS>::OpEssentialLhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
1;
auto out_meshset = get_temp_meshset_ptr(moab);
CHKERR moab.add_entities(*out_meshset, r);
CHKERR moab.write_file(
name.c_str(),
"VTK",
"", out_meshset->get_ptr(), 1);
};
private:
boost::shared_ptr<FieldEvaluatorInterface::SetPtsData>
fieldEvalData;
struct BlockedThermalParameters
: public boost::enable_shared_from_this<BlockedThermalParameters> {
}
return boost::shared_ptr<double>(shared_from_this(), &
heatCapacity);
}
};
struct BlockedThermoElasticParameters
: public boost::enable_shared_from_this<BlockedThermoElasticParameters> {
return boost::shared_ptr<VectorDouble>(shared_from_this(),
}
return boost::shared_ptr<double>(shared_from_this(), &
refTemp);
}
};
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
std::string block_name,
boost::shared_ptr<BlockedThermalParameters> blockedParamsPtr,
Sev sev);
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
std::string block_name,
boost::shared_ptr<BlockedThermoElasticParameters> blockedParamsPtr,
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
boost::shared_ptr<HenckyOps::CommonData> elastic_common_ptr,
boost::shared_ptr<ThermoElasticProblem::BlockedThermalParameters>
thermal_common_ptr,
boost::shared_ptr<ThermoElasticProblem::BlockedThermoElasticParameters>
thermoelastic_common_ptr,
auto vec_temp_ptr = boost::make_shared<VectorDouble>();
auto coeff_expansion_ptr = thermoelastic_common_ptr->getCoeffExpansionPtr();
auto ref_temp_ptr = thermoelastic_common_ptr->getRefTempPtr();
pip.push_back(
new typename H::template OpCalculateHenckyThermalStress<DIM, I, 0>(
"U", vec_temp_ptr, elastic_common_ptr, coeff_expansion_ptr,
ref_temp_ptr));
"U", elastic_common_ptr->getMatFirstPiolaStress()));
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
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,
auto elastic_common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
auto thermal_common_ptr = boost::make_shared<BlockedThermalParameters>();
Sev::inform);
auto thermoelastic_common_ptr =
boost::make_shared<BlockedThermoElasticParameters>();
thermoelastic_common_ptr, Sev::inform);
CHKERR opThermoElasticFactoryDomainRhs<DIM, A, I, DomainEleOp>(
m_field, pip,
field_name, elastic_common_ptr, thermal_common_ptr,
thermoelastic_common_ptr, sev);
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
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,
using OpKPiola =
typename B::template OpGradTensorGrad<1, DIM, DIM, 1>;
auto vec_temp_ptr = boost::make_shared<VectorDouble>();
auto coeff_expansion_ptr = thermoelastic_common_ptr->getCoeffExpansionPtr();
auto ref_temp_ptr = thermoelastic_common_ptr->getRefTempPtr();
pip.push_back(
new typename H::template OpCalculateHenckyThermalStress<DIM, I, 0>(
"U", vec_temp_ptr, elastic_common_ptr, coeff_expansion_ptr,
ref_temp_ptr));
pip.push_back(new typename H::template OpHenckyTangent<DIM, I, 0>(
elastic_common_ptr->getMatTangent()));
pip.push_back(new typename H::template OpCalculateHenckyThermalStressdT<
field_name, coupled_field_name, elastic_common_ptr,
coeff_expansion_ptr));
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
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) {
auto elastic_common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
auto thermal_common_ptr = boost::make_shared<BlockedThermalParameters>();
Sev::inform);
auto thermoelastic_common_ptr =
boost::make_shared<BlockedThermoElasticParameters>();
thermoelastic_common_ptr, Sev::inform);
CHKERR opThermoElasticFactoryDomainLhs<DIM, A, I, DomainEleOp>(
m_field, pip,
field_name, coupled_field_name, elastic_common_ptr,
thermal_common_ptr, thermoelastic_common_ptr, sev);
}
};
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
std::string block_name,
boost::shared_ptr<BlockedThermalParameters> blockedParamsPtr,
Sev sev) {
OpMatThermalBlocks(boost::shared_ptr<double> conductivity_ptr,
boost::shared_ptr<double> capacity_ptr,
std::vector<const CubitMeshSets *> meshset_vec_ptr)
conductivityPtr(conductivity_ptr), capacityPtr(capacity_ptr) {
"Cannot get data from thermal block");
}
for (auto &b : blockData) {
if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
*conductivityPtr = b.conductivity;
*capacityPtr = b.capacity;
}
}
}
private:
double conductivity;
double capacity;
};
std::vector<BlockData> blockData;
std::vector<const CubitMeshSets *> meshset_vec_ptr,
for (
auto m : meshset_vec_ptr) {
std::vector<double> block_data;
CHKERR m->getAttributes(block_data);
if (block_data.size() < 2) {
"Expected that block has at least two attributes");
}
auto get_block_ents = [&]() {
m_field.
get_moab().get_entities_by_handle(
m->meshset, ents,
true);
return ents;
};
<<
m->getName() <<
": conductivity = " << block_data[0]
<< " capacity = " << block_data[1];
blockData.push_back({block_data[0], block_data[1], get_block_ents()});
}
}
boost::shared_ptr<double> conductivityPtr;
boost::shared_ptr<double> capacityPtr;
};
pipeline.push_back(new OpMatThermalBlocks(
blockedParamsPtr->getHeatConductivityPtr(),
blockedParamsPtr->getHeatCapacityPtr(),
mField, sev,
(boost::format("%s(.*)") % block_name).str()
))
));
}
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
std::string block_name,
boost::shared_ptr<BlockedThermoElasticParameters> blockedParamsPtr,
OpMatThermoElasticBlocks(boost::shared_ptr<VectorDouble> expansion_ptr,
boost::shared_ptr<double> ref_temp_ptr,
std::vector<const CubitMeshSets *> meshset_vec_ptr)
expansionPtr(expansion_ptr), refTempPtr(ref_temp_ptr) {
extractThermoElasticBlockData(m_field, meshset_vec_ptr, sev),
"Cannot get data from thermoelastic block");
}
for (auto &b : blockData) {
if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
*expansionPtr = b.expansion;
*refTempPtr = b.ref_temp;
}
}
}
private:
double ref_temp;
};
std::vector<BlockData> blockData;
std::vector<const CubitMeshSets *> meshset_vec_ptr,
Sev sev) {
for (
auto m : meshset_vec_ptr) {
std::vector<double> block_data;
CHKERR m->getAttributes(block_data);
if (block_data.size() < 2) {
"Expected that block has at least two attributes");
}
auto get_block_ents = [&]() {
m_field.
get_moab().get_entities_by_handle(
m->meshset, ents,
true);
return ents;
};
auto get_expansion = [&]() {
if (block_data.size() > 2) {
expansion[1] = block_data[2];
}
if (
SPACE_DIM == 3 && block_data.size() > 3) {
expansion[2] = block_data[3];
}
return expansion;
};
auto coeff_exp_vec = get_expansion();
<< " ref_temp = " << block_data[0]
<< " expansion = " << coeff_exp_vec;
blockData.push_back({block_data[0], coeff_exp_vec, get_block_ents()});
}
}
boost::shared_ptr<VectorDouble> expansionPtr;
boost::shared_ptr<double> refTempPtr;
};
pipeline.push_back(new OpMatThermoElasticBlocks(
blockedParamsPtr->getCoeffExpansionPtr(),
blockedParamsPtr->getRefTempPtr(),
mField, sev,
(boost::format("%s(.*)") % block_name).str()
))
));
}
}
auto get_command_line_parameters = [&]() {
PETSC_NULL);
PETSC_NULL);
PETSC_NULL);
PETSC_NULL);
PETSC_NULL);
PETSC_NULL);
PETSC_NULL);
} else {
}
PETSC_NULL);
};
CHKERR get_command_line_parameters();
}
CHKERR simple->addBoundaryField(
"FLUX", flux_space, base, 1);
auto no_rule = [](
int,
int,
int) {
return -1; };
field_eval_fe_ptr->getRuleHook = no_rule;
auto block_thermal_params = boost::make_shared<BlockedThermalParameters>();
auto block_thermoelastic_params =
boost::make_shared<BlockedThermoElasticParameters>();
auto coeff_expansion_ptr =
block_thermoelastic_params->getCoeffExpansionPtr();
auto ref_temp_ptr = block_thermoelastic_params->getRefTempPtr();
"MAT_THERMAL", block_thermal_params,
Sev::verbose);
field_eval_fe_ptr->getOpPtrVector(), "MAT_THERMOELASTIC",
block_thermoelastic_params, Sev::verbose);
field_eval_fe_ptr->getOpPtrVector(), {H1, HDIV});
auto hencky_common_data_ptr =
HenckyOps::commonDataFactory<SPACE_DIM, IT, DomainEleOp>(
mField, field_eval_fe_ptr->getOpPtrVector(),
"U",
"MAT_ELASTIC",
Sev::inform);
auto mat_D_ptr = hencky_common_data_ptr->matDPtr;
auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
field_eval_fe_ptr->getOpPtrVector().push_back(
field_eval_fe_ptr->getOpPtrVector().push_back(
field_eval_fe_ptr->getOpPtrVector().push_back(
field_eval_fe_ptr->getOpPtrVector().push_back(
field_eval_fe_ptr->getOpPtrVector().push_back(
new
typename H::template OpCalculateHenckyThermalStress<SPACE_DIM, IT, 0>(
"U",
tempFieldPtr, hencky_common_data_ptr, coeff_expansion_ptr,
ref_temp_ptr));
field_eval_fe_ptr->getOpPtrVector().push_back(
hencky_common_data_ptr->getMatFirstPiolaStress(),
field_eval_fe_ptr->getOpPtrVector().push_back(
} else {
field_eval_fe_ptr->getOpPtrVector().push_back(
new typename H::template OpCalculateLogC<SPACE_DIM, IT>(
"U", hencky_common_data_ptr));
};
}
}
auto get_skin = [&]() {
CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
return skin_ents;
};
auto filter_flux_blocks = [&](auto skin, bool temp_bc = false) {
auto remove_cubit_blocks = [&](
auto c) {
) {
skin = subtract(skin, ents);
}
};
auto remove_named_blocks = [&](
auto n) {
(boost::format(
"%s(.*)") %
n).str()
))
) {
skin = subtract(skin, ents);
}
};
if (!temp_bc) {
"remove_cubit_blocks");
"remove_named_blocks");
}
"remove_cubit_blocks");
return skin;
};
auto filter_true_skin = [&](auto skin) {
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
return boundary_ents;
};
auto remove_flux_ents = filter_true_skin(filter_flux_blocks(get_skin()));
auto remove_temp_bc_ents =
filter_true_skin(filter_flux_blocks(get_skin(), true));
remove_flux_ents);
remove_temp_bc_ents);
MOFEM_LOG(
"SYNC", Sev::noisy) << remove_flux_ents << endl;
MOFEM_LOG(
"SYNC", Sev::noisy) << remove_temp_bc_ents << endl;
#ifndef NDEBUG
remove_flux_ents);
remove_temp_bc_ents);
#endif
simple->getProblemName(),
"FLUX", remove_flux_ents);
simple->getProblemName(),
"TBC", remove_temp_bc_ents);
auto set_init_temp = [](boost::shared_ptr<FieldEntity> field_entity_ptr) {
return 0;
};
"T");
simple->getProblemName(),
"U");
simple->getProblemName(),
"FLUX",
false);
}
auto boundary_marker =
bc_mng->getMergedBlocksMarker(vector<string>{"HEATFLUX"});
};
auto block_thermal_params = boost::make_shared<BlockedThermalParameters>();
auto heat_conductivity_ptr = block_thermal_params->getHeatConductivityPtr();
auto heat_capacity_ptr = block_thermal_params->getHeatCapacityPtr();
auto block_thermoelastic_params =
boost::make_shared<BlockedThermoElasticParameters>();
auto coeff_expansion_ptr = block_thermoelastic_params->getCoeffExpansionPtr();
auto ref_temp_ptr = block_thermoelastic_params->getRefTempPtr();
auto time_scale =
boost::make_shared<TimeScale>("", false, [](const double) { return 1; });
auto def_time_scale = [time_scale](const double time) {
return (!time_scale->argFileScale) ? time : 1;
};
auto def_file_name = [time_scale](
const std::string &&
name) {
return (!time_scale->argFileScale) ?
name :
"";
};
auto time_bodyforce_scaling = boost::make_shared<TimeScale>(
def_file_name("bodyforce_scale.txt"), false, def_time_scale);
auto time_heatsource_scaling = boost::make_shared<TimeScale>(
def_file_name("heatsource_scale.txt"), false, def_time_scale);
auto time_temperature_scaling = boost::make_shared<TimeScale>(
def_file_name("temperature_bc_scale.txt"), false, def_time_scale);
auto time_displacement_scaling = boost::make_shared<TimeScale>(
def_file_name("displacement_bc_scale.txt"), false, def_time_scale);
auto time_flux_scaling = boost::make_shared<TimeScale>(
def_file_name("flux_bc_scale.txt"), false, def_time_scale);
auto time_force_scaling = boost::make_shared<TimeScale>(
def_file_name("force_bc_scale.txt"), false, def_time_scale);
auto add_domain_rhs_ops = [&](auto &pipeline) {
Sev::inform);
block_thermoelastic_params, Sev::inform);
auto hencky_common_data_ptr =
HenckyOps::commonDataFactory<SPACE_DIM, IT, DomainEleOp>(
mField, pipeline,
"U",
"MAT_ELASTIC", Sev::inform);
auto mat_D_ptr = hencky_common_data_ptr->matDPtr;
auto mat_grad_ptr = hencky_common_data_ptr->matGradPtr;
auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
auto vec_temp_ptr = boost::make_shared<VectorDouble>();
auto vec_temp_dot_ptr = boost::make_shared<VectorDouble>();
auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
auto vec_temp_div_ptr = boost::make_shared<VectorDouble>();
pipeline.push_back(
"FLUX", vec_temp_div_ptr));
pipeline.push_back(
CHKERR opThermoElasticFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
mField, pipeline,
"U",
"MAT_ELASTIC",
"MAT_THERMAL",
"MAT_THERMOELASTIC", Sev::inform);
auto resistance = [heat_conductivity_ptr](
const double,
const double,
return (1. / (*heat_conductivity_ptr));
};
auto capacity = [heat_capacity_ptr](
const double,
const double,
return -(*heat_capacity_ptr);
};
return -1.;
};
pipeline.push_back(
new OpHdivFlux(
"FLUX", mat_flux_ptr, resistance, mat_grad_ptr));
pipeline.push_back(
new OpHDivTemp(
"FLUX", vec_temp_ptr, unity));
pipeline.push_back(
new OpBaseDivFlux(
"T", vec_temp_div_ptr, unity));
pipeline.push_back(
new OpBaseDotT(
"T", vec_temp_dot_ptr, capacity));
pipeline,
mField,
"T", {time_scale, time_heatsource_scaling},
"HEAT_SOURCE", Sev::inform);
pipeline,
mField,
"U", {time_scale, time_bodyforce_scaling},
"BODY_FORCE", Sev::inform);
pipeline,
mField,
"T", vec_temp_ptr,
"SETTEMP", Sev::inform);
};
auto add_domain_lhs_ops = [&](auto &pipeline) {
Sev::verbose);
block_thermoelastic_params,
Sev::verbose);
auto hencky_common_data_ptr =
HenckyOps::commonDataFactory<SPACE_DIM, IT, DomainEleOp>(
mField, pipeline,
"U",
"MAT_ELASTIC", Sev::inform, 1);
auto mat_D_ptr = hencky_common_data_ptr->matDPtr;
auto mat_grad_ptr = hencky_common_data_ptr->matGradPtr;
auto resistance = [heat_conductivity_ptr](
const double,
const double,
return (1. / (*heat_conductivity_ptr));
};
auto capacity = [heat_capacity_ptr](
const double,
const double,
return -(*heat_capacity_ptr);
};
pipeline.push_back(
new OpHdivHdiv(
"FLUX",
"FLUX", resistance, mat_grad_ptr));
pipeline.push_back(
new OpHdivT(
"FLUX",
"T", []()
constexpr {
return -1; },
true));
auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(
pipeline.push_back(
new OpHdivU(
"FLUX",
"U", mat_flux_ptr, resistance, mat_grad_ptr));
CHKERR opThermoElasticFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
mField, pipeline,
"U",
"T",
"MAT_ELASTIC",
"MAT_THERMAL",
"MAT_THERMOELASTIC", Sev::inform);
auto op_capacity =
new OpCapacity(
"T",
"T", capacity);
op_capacity->feScalingFun = [](
const FEMethod *fe_ptr) {
return fe_ptr->ts_a;
};
pipeline.push_back(op_capacity);
pipeline,
mField,
"T",
"SETTEMP",
[](
const FEMethod *fe_ptr) {
return fe_ptr->ts_dt; }, Sev::verbose);
};
auto add_boundary_rhs_ops = [&](auto &pipeline) {
pipeline,
mField,
"U", {time_scale, time_force_scaling},
"FORCE",
Sev::inform);
pipeline,
mField,
"FLUX", {time_scale, time_temperature_scaling},
"TEMPERATURE", Sev::inform);
using OpFluxBC =
pipeline,
mField,
"TBC", {time_scale, time_flux_scaling},
"FLUX",
Sev::inform);
using OpConvectionRhsBC =
T::OpFlux<ThermoElasticOps::ConvectionBcType<BLOCKSET>, 1, 1>;
using OpRadiationRhsBC =
T::OpFlux<ThermoElasticOps::RadiationBcType<BLOCKSET>, 1, 1>;
auto temp_bc_ptr = boost::make_shared<VectorDouble>();
T::AddFluxToPipeline<OpConvectionRhsBC>::add(
pipeline,
mField,
"TBC", temp_bc_ptr,
"CONVECTION", Sev::inform);
T::AddFluxToPipeline<OpRadiationRhsBC>::add(
pipeline,
mField,
"TBC", temp_bc_ptr,
"RADIATION", Sev::inform);
auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(
struct OpTBCTimesNormalFlux
OpTBCTimesNormalFlux(
const std::string
field_name,
boost::shared_ptr<MatrixDouble> vec,
boost::shared_ptr<Range> ents_ptr = nullptr)
auto t_w = OP::getFTensor0IntegrationWeight();
auto t_normal = OP::getFTensor1NormalsAtGaussPts();
auto t_vec = getFTensor1FromMat<SPACE_DIM, 1>(*sourceVec);
for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
const double alpha = t_w * (t_vec(
i) * t_normal(
i));
int rr = 0;
for (; rr != OP::nbRows; ++rr) {
OP::locF[rr] += alpha * t_row_base;
++t_row_base;
}
for (; rr < OP::nbRowBaseFunctions; ++rr)
++t_row_base;
++t_w;
++t_vec;
++t_normal;
}
EntityType fe_type = OP::getNumeredEntFiniteElementPtr()->getEntType();
if (fe_type == MBTRI) {
OP::locF /= 2;
}
}
protected:
boost::shared_ptr<MatrixDouble> sourceVec;
};
pipeline.push_back(new OpTBCTimesNormalFlux("TBC", mat_flux_ptr));
struct OpBaseNormalTimesTBC
OpBaseNormalTimesTBC(
const std::string
field_name,
boost::shared_ptr<VectorDouble> vec,
boost::shared_ptr<Range> ents_ptr = nullptr)
auto t_w = OP::getFTensor0IntegrationWeight();
auto t_normal = OP::getFTensor1NormalsAtGaussPts();
auto t_vec = getFTensor0FromVec(*sourceVec);
for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
const double alpha = t_w * t_vec;
int rr = 0;
for (; rr != OP::nbRows; ++rr) {
OP::locF[rr] += alpha * (t_row_base(
i) * t_normal(
i));
++t_row_base;
}
for (; rr < OP::nbRowBaseFunctions; ++rr)
++t_row_base;
++t_w;
++t_vec;
++t_normal;
}
EntityType fe_type = OP::getNumeredEntFiniteElementPtr()->getEntType();
if (fe_type == MBTRI) {
OP::locF /= 2;
}
}
protected:
boost::shared_ptr<VectorDouble> sourceVec;
};
pipeline.push_back(new OpBaseNormalTimesTBC("FLUX", temp_bc_ptr));
};
auto add_boundary_lhs_ops = [&](auto &pipeline) {
using T =
using OpConvectionLhsBC =
T::OpFlux<ThermoElasticOps::ConvectionBcType<BLOCKSET>, 1, 1>;
using OpRadiationLhsBC =
T::OpFlux<ThermoElasticOps::RadiationBcType<BLOCKSET>, 1, 1>;
auto temp_bc_ptr = boost::make_shared<VectorDouble>();
T::AddFluxToPipeline<OpConvectionLhsBC>::add(pipeline,
mField,
"TBC",
"TBC",
"CONVECTION", Sev::verbose);
T::AddFluxToPipeline<OpRadiationLhsBC>::add(
pipeline,
mField,
"TBC",
"TBC", temp_bc_ptr,
"RADIATION", Sev::verbose);
struct OpTBCTimesNormalFlux
OpTBCTimesNormalFlux(const std::string row_field_name,
const std::string col_field_name,
boost::shared_ptr<Range> ents_ptr = nullptr)
:
OP(row_field_name, col_field_name,
OP::OPROWCOL, ents_ptr) {
this->sYmm = false;
this->assembleTranspose = true;
this->onlyTranspose = false;
}
auto t_w = OP::getFTensor0IntegrationWeight();
auto t_normal = OP::getFTensor1NormalsAtGaussPts();
for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
auto a_mat_ptr = &*OP::locMat.data().begin();
int rr = 0;
for (; rr != OP::nbRows; rr++) {
const double alpha = t_w * t_row_base;
for (int cc = 0; cc != OP::nbCols; cc++) {
*a_mat_ptr += alpha * (t_col_base(
i) * t_normal(
i));
++t_col_base;
++a_mat_ptr;
}
++t_row_base;
}
for (; rr < OP::nbRowBaseFunctions; ++rr)
++t_row_base;
++t_normal;
++t_w;
}
EntityType fe_type = OP::getNumeredEntFiniteElementPtr()->getEntType();
if (fe_type == MBTRI) {
OP::locMat /= 2;
}
}
};
pipeline.push_back(new OpTBCTimesNormalFlux("TBC", "FLUX"));
};
auto get_bc_hook_rhs = [&]() {
mField, pipeline_mng->getDomainRhsFE(),
{time_scale, time_displacement_scaling});
};
auto get_bc_hook_lhs = [&]() {
mField, pipeline_mng->getDomainLhsFE(),
{time_scale, time_displacement_scaling});
};
pipeline_mng->getDomainRhsFE()->preProcessHook = get_bc_hook_rhs();
pipeline_mng->getDomainLhsFE()->preProcessHook = get_bc_hook_lhs();
CHKERR add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
CHKERR add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
CHKERR add_boundary_rhs_ops(pipeline_mng->getOpBoundaryRhsPipeline());
CHKERR add_boundary_lhs_ops(pipeline_mng->getOpBoundaryLhsPipeline());
}
auto solver = pipeline_mng->createTSIM();
auto set_section_monitor = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
(void *)(snes_ctx_ptr.get()), nullptr);
};
auto create_post_process_elements = [&]() {
auto block_thermal_params = boost::make_shared<BlockedThermalParameters>();
auto block_thermoelastic_params =
boost::make_shared<BlockedThermoElasticParameters>();
auto coeff_expansion_ptr =
block_thermoelastic_params->getCoeffExpansionPtr();
auto ref_temp_ptr = block_thermoelastic_params->getRefTempPtr();
auto u_ptr = boost::make_shared<MatrixDouble>();
auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
auto vec_temp_ptr = boost::make_shared<VectorDouble>();
auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
auto push_domain_ops = [&](auto &pp_fe) {
auto &pip = pp_fe->getOpPtrVector();
Sev::verbose);
pip, "MAT_THERMOELASTIC", block_thermoelastic_params, Sev::verbose);
pip.push_back(
"U", mat_grad_ptr));
auto elastic_common_ptr = commonDataFactory<SPACE_DIM, IT, DomainEleOp>(
mField, pip,
"U",
"MAT_ELASTIC", Sev::inform);
pip.push_back(
new
typename H::template OpCalculateHenckyThermalStress<SPACE_DIM, IT, 0>(
"U", vec_temp_ptr, elastic_common_ptr, coeff_expansion_ptr,
ref_temp_ptr));
elastic_common_ptr->getMatFirstPiolaStress(), mat_stress_ptr));
pip.push_back(
} else {
mat_stress_ptr = elastic_common_ptr->getMatFirstPiolaStress();
mat_strain_ptr = elastic_common_ptr->getMatLogC();
}
};
auto push_post_proc_ops = [&](auto &pp_fe) {
auto &pip = pp_fe->getOpPtrVector();
pip.push_back(
pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
{{"T", vec_temp_ptr}},
{{"U", u_ptr}, {"FLUX", mat_flux_ptr}},
{},
{{"CAUCHY", mat_stress_ptr}, {"STRAIN", mat_strain_ptr}}
)
);
} else {
pip.push_back(
pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
{{"T", vec_temp_ptr}},
{{"U", u_ptr}, {"FLUX", mat_flux_ptr}},
{{"PIOLA", mat_stress_ptr}},
{{"HENCKY_STRAIN", mat_strain_ptr}}
)
);
}
};
auto domain_post_proc = [&]() {
return boost::shared_ptr<PostProcEle>();
auto pp_fe = boost::make_shared<PostProcEle>(mField);
"push domain ops to domain element");
"push post proc ops to domain element");
return pp_fe;
};
auto skin_post_proc = [&]() {
return boost::shared_ptr<SkinPostProcEle>();
auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
"push domain ops to side element");
pp_fe->getOpPtrVector().push_back(op_side);
"push post proc ops to skin element");
return pp_fe;
};
return std::make_pair(domain_post_proc(), skin_post_proc());
};
auto monitor_ptr = boost::make_shared<FEMethod>();
auto set_time_monitor = [&](auto dm, auto solver, auto domain_post_proc_fe,
auto skin_post_proc_fe) {
monitor_ptr->preProcessHook = [&]() {
domain_post_proc_fe,
monitor_ptr->getCacheWeakPtr());
CHKERR domain_post_proc_fe->writeFile(
"out_" + boost::lexical_cast<std::string>(monitor_ptr->ts_step) +
".h5m");
}
skin_post_proc_fe,
monitor_ptr->getCacheWeakPtr());
CHKERR skin_post_proc_fe->writeFile(
"out_skin_" +
boost::lexical_cast<std::string>(monitor_ptr->ts_step) + ".h5m");
}
}
->evalFEAtThePoint<SPACE_DIM>(
fieldEvalCoords.data(), 1e-12,
simple->getProblemName(),
simple->getDomainFEName(), fieldEvalData,
auto eval_num_vec =
CHKERR VecZeroEntries(eval_num_vec);
if (tempFieldPtr->size()) {
CHKERR VecSetValue(eval_num_vec, 0, 1, ADD_VALUES);
}
CHKERR VecAssemblyBegin(eval_num_vec);
CHKERR VecAssemblyEnd(eval_num_vec);
double eval_num;
CHKERR VecSum(eval_num_vec, &eval_num);
if (!(int)eval_num) {
"atom test %d failed: did not find elements to evaluate "
"the field, check the coordinates",
}
}
if (tempFieldPtr->size()) {
auto t_temp = getFTensor0FromVec(*tempFieldPtr);
<< "Eval point T: " << t_temp;
if (
atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
if (
atom_test <= 3 && fabs(t_temp - 554.48) > 1e-2) {
"atom test %d failed: wrong temperature value",
}
if (
atom_test == 4 && fabs(t_temp - 250) > 1e-2) {
"atom test %d failed: wrong temperature value",
}
if (
atom_test == 5 && fabs(t_temp - 1) > 1e-2) {
"atom test %d failed: wrong temperature value",
}
}
}
if (fluxFieldPtr->size1()) {
auto t_flux = getFTensor1FromMat<SPACE_DIM>(*fluxFieldPtr);
auto flux_mag = sqrt(t_flux(
i) * t_flux(
i));
<< "Eval point FLUX magnitude: " << flux_mag;
if (
atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
if (
atom_test <= 3 && fabs(flux_mag - 27008.0) > 2e1) {
"atom test %d failed: wrong flux value",
atom_test);
}
if (
atom_test == 4 && fabs(flux_mag - 150e3) > 1e-1) {
"atom test %d failed: wrong flux value",
atom_test);
}
if (
atom_test == 5 && fabs(flux_mag) > 1e-6) {
"atom test %d failed: wrong flux value",
atom_test);
}
}
}
if (dispFieldPtr->size1()) {
auto t_disp = getFTensor1FromMat<SPACE_DIM>(*dispFieldPtr);
auto disp_mag = sqrt(t_disp(
i) * t_disp(
i));
<< "Eval point U magnitude: " << disp_mag;
if (
atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
if (
atom_test == 1 && fabs(disp_mag - 0.00345) > 1e-5) {
"atom test %d failed: wrong displacement value",
}
fabs(disp_mag - 0.00265) > 1e-5) {
"atom test %d failed: wrong displacement value",
}
fabs(t_disp(0) - (std::sqrt(std::exp(0.2)) - 1)) > 1e-5 &&
fabs(t_disp(1) - (std::sqrt(std::exp(0.2)) - 1)) > 1e-5) {
"atom test %d failed: wrong displacement value",
}
}
}
if (strainFieldPtr->size1()) {
auto t_strain =
getFTensor2SymmetricFromMat<SPACE_DIM>(*strainFieldPtr);
auto t_strain_trace = t_strain(
i,
i);
if (
atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
if (
atom_test == 1 && fabs(t_strain_trace - 0.00679) > 1e-5) {
"atom test %d failed: wrong strain value",
atom_test);
}
fabs(t_strain_trace - 0.00522) > 1e-5) {
"atom test %d failed: wrong strain value",
atom_test);
}
if ((
atom_test == 5) && fabs(t_strain_trace - 0.2) > 1e-5) {
"atom test %d failed: wrong strain value",
atom_test);
}
}
}
if (stressFieldPtr->size1()) {
double von_mises_stress;
auto getVonMisesStress = [&](auto t_stress) {
von_mises_stress = std::sqrt(
0.5 *
((t_stress(0, 0) - t_stress(1, 1)) *
(t_stress(0, 0) - t_stress(1, 1)) +
(
SPACE_DIM == 3 ? (t_stress(1, 1) - t_stress(2, 2)) *
(t_stress(1, 1) - t_stress(2, 2))
: 0) +
(
SPACE_DIM == 3 ? (t_stress(2, 2) - t_stress(0, 0)) *
(t_stress(2, 2) - t_stress(0, 0))
: 0) +
6 * (t_stress(0, 1) * t_stress(0, 1) +
(
SPACE_DIM == 3 ? t_stress(1, 2) * t_stress(1, 2) : 0) +
(
SPACE_DIM == 3 ? t_stress(2, 0) * t_stress(2, 0) : 0))));
};
auto t_stress =
getFTensor2SymmetricFromMat<SPACE_DIM>(*stressFieldPtr);
CHKERR getVonMisesStress(t_stress);
} else {
auto t_stress =
getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*stressFieldPtr);
CHKERR getVonMisesStress(t_stress);
}
<< "Eval point von Mises Stress: " << von_mises_stress;
if (
atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
if (
atom_test == 1 && fabs(von_mises_stress - 523.0) > 5e-1) {
"atom test %d failed: wrong von Mises stress value",
}
if (
atom_test == 2 && fabs(von_mises_stress - 16.3) > 5e-2) {
"atom test %d failed: wrong von Mises stress value",
}
if (
atom_test == 3 && fabs(von_mises_stress - 14.9) > 5e-2) {
"atom test %d failed: wrong von Mises stress value",
}
if (
atom_test == 5 && fabs(von_mises_stress) > 5e-2) {
"atom test %d failed: wrong von Mises stress value",
}
}
}
}
};
auto null = boost::shared_ptr<FEMethod>();
monitor_ptr, null);
};
auto set_fieldsplit_preconditioner = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
PC pc;
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs == PETSC_TRUE) {
auto name_prb =
simple->getProblemName();
CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"U", 0,
CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"FLUX", 0, 0,
is_flux);
CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"T", 0, 0,
is_T);
CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"TBC", 0, 0,
is_TBC);
IS is_tmp, is_tmp2;
CHKERR ISExpand(is_T, is_flux, &is_tmp);
CHKERR ISExpand(is_TBC, is_tmp, &is_tmp2);
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_TFlux);
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
}
};
CHKERR TSSetFromOptions(solver);
CHKERR set_section_monitor(solver);
CHKERR set_fieldsplit_preconditioner(solver);
auto [domain_post_proc_fe, skin_post_proc_fe] =
create_post_process_elements();
CHKERR set_time_monitor(dm, solver, domain_post_proc_fe, skin_post_proc_fe);
}
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
const char param_file[] = "param_file.petsc";
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "ThermoElastic"));
LogManager::setLog("ThermoElastic");
core_log->add_sink(
LogManager::createSink(LogManager::getStrmSync(), "ThermoElasticSync"));
LogManager::setLog("ThermoElasticSync");
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
}
}
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)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#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, RowColData rc=RowColData::COL)
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, RowColData rc=RowColData::COL)
Get smart vector from DM.
#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
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.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
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)
OpCalculateQdotQLhs_dU< SPACE_DIM, GAUSS, AssemblyDomainEleOp, IS_LARGE_STRAINS > OpHdivU
Integrate Lhs of flux term coupled to displacement field.
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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 1, SPACE_DIM > OpHDivTemp
Integrate Rhs div flux base times temperature (T)
OpBaseDotT OpBaseDivFlux
Integrate Rhs base of temperature times divergence of flux (T)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBaseDotT
Integrate Rhs base of temperature time heat capacity times heat rate (T)
OpHdivHdivImpl< SPACE_DIM, IS_LARGE_STRAINS > OpHdivHdiv
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
constexpr IntegrationType I
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
PostProcBrokenMeshInMoab< BoundaryEle > SkinPostProcEle
#define EXECUTABLE_DIMENSION
ElementsAndOps< SPACE_DIM >::SideEle SideEle
constexpr auto field_name
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
EssentialBC< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpEssentialLhs< HeatFluxCubitBcData, 3, SPACE_DIM > OpEssentialFluxLhs
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
EssentialBC< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpEssentialRhs< HeatFluxCubitBcData, 3, SPACE_DIM > OpEssentialFluxRhs
[Natural boundary conditions]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, -1 > OpKPiola
[Only used for dynamics]
BoundaryNaturalBC::OpFlux< NaturalTemperatureMeshsets, 3, SPACE_DIM > OpTemperatureBC
FTensor::Index< 'm', 3 > m
Add operators pushing bases from local to physical configuration.
Boundary condition manager for finite element problem setup.
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.
auto getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Class (Function) to enforce essential constrains.
Field evaluator interface.
Definition of the heat flux bc data structure.
Section manager is used to create indexes and sections.
Interface for managing meshsets containing materials and boundary conditions.
Natural boundary conditions.
Get vector field for H-div approximation.
Calculate divergence of vector field.
Specialization for double precision scalar field values calculation.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Specialization for MatrixDouble vector field values calculation.
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
Operator for symmetrizing tensor fields.
Template struct for dimension-specific finite element types.
PipelineManager interface.
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
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, 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)
constexpr AssemblyType AT
constexpr IntegrationType IT
#define FINITE_DEFORMATION_FLAG
NaturalBC< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS > BoundaryNaturalBC
[Body and heat source]
constexpr bool IS_LARGE_STRAINS
BoundaryNaturalBC::OpFlux< NaturalForceMeshsets, 1, SPACE_DIM > OpForce
NaturalBC< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS > DomainNaturalBCRhs
[Thermal problem]
double default_heat_capacity
DomainNaturalBCRhs::OpFlux< NaturalMeshsetType< BLOCKSET >, 1, 1 > OpHeatSource
double default_young_modulus
[Natural boundary conditions]
double default_coeff_expansion
PetscBool do_output_domain
double default_heat_conductivity
NaturalBC< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS > DomainNaturalBCLhs
DomainNaturalBCRhs::OpFlux< NaturalMeshsetType< BLOCKSET >, 1, SPACE_DIM > OpBodyForce
double default_poisson_ratio
using ScalerFunTwoArgs = boost::function<
double(
double,
double)>;
struct BlockedThermalParameters
: public boost::enable_shared_from_this<BlockedThermalParameters> {
double heatConductivity;
double heatConductivityScaling;
double heatCapacity;
double heatCapacityScaling;
inline auto getHeatConductivityPtr() {
return boost::shared_ptr<double>(shared_from_this(), &heatConductivity);
}
inline auto getHeatConductivityScalingPtr() {
return boost::shared_ptr<double>(shared_from_this(),
&heatConductivityScaling);
}
inline auto getHeatCapacityPtr() {
return boost::shared_ptr<double>(shared_from_this(), &heatCapacity);
}
inline auto getHeatCapacityScalingPtr() {
return boost::shared_ptr<double>(shared_from_this(), &heatCapacityScaling);
}
};
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
std::string block_name,
boost::shared_ptr<BlockedThermalParameters> blockedParamsPtr,
double default_thermal_conductivity_scale,
double default_thermal_capacity_scale, Sev sev,
ScalerFunTwoArgs thermal_conductivity_scaling_func,
ScalerFunTwoArgs heat_capacity_scaling_func) {
double local_heat_conductivity_scale = default_thermal_conductivity_scale;
double local_heat_capacity_scale = default_thermal_capacity_scale;
OpMatThermalBlocks(boost::shared_ptr<double> conductivity_ptr,
boost::shared_ptr<double> capacity_ptr,
boost::shared_ptr<double> heat_conductivity_scaling_ptr,
boost::shared_ptr<double> heat_capacity_scaling_ptr,
double local_heat_conductivity,
double local_heat_capacity,
double local_heat_conductivity_scale,
double local_heat_capacity_scale,
const ScalerFunTwoArgs &thermal_conductivity_scaling,
const ScalerFunTwoArgs &heat_capacity_scaling,
std::vector<const CubitMeshSets *> meshset_vec_ptr)
conductivityPtr(conductivity_ptr), capacityPtr(capacity_ptr),
conductivityScalingPtr(heat_conductivity_scaling_ptr),
capacityScalingPtr(heat_capacity_scaling_ptr),
defaultHeatConductivity(local_heat_conductivity),
defaultHeatCapacity(local_heat_capacity),
defaultHeatConductivityScale(local_heat_conductivity_scale),
defaultHeatCapacityScale(local_heat_capacity_scale),
thermalConductivityScaling(thermal_conductivity_scaling),
heatCapacityScaling(heat_capacity_scaling) {
"Cannot get data from thermal block");
}
EntitiesFieldData::EntData &data) {
auto scale_param = [this](auto parameter, double scaling) {
return parameter * scaling;
};
for (auto &b : blockData) {
if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
*conductivityScalingPtr =
thermalConductivityScaling(b.conductivity, b.capacity);
*conductivityPtr =
scale_param(b.conductivity, *conductivityScalingPtr);
*capacityScalingPtr = heatCapacityScaling(b.conductivity, b.capacity);
*capacityPtr = scale_param(b.capacity, *capacityScalingPtr);
}
}
*conductivityScalingPtr = thermalConductivityScaling(
defaultHeatConductivity / defaultHeatConductivityScale,
defaultHeatCapacity /
defaultHeatCapacityScale);
*conductivityPtr =
scale_param(defaultHeatConductivity, *conductivityScalingPtr);
*capacityScalingPtr = heatCapacityScaling(
defaultHeatConductivity / defaultHeatConductivityScale,
defaultHeatCapacity /
defaultHeatCapacityScale);
*capacityPtr = scale_param(defaultHeatCapacity, *capacityScalingPtr);
}
private:
double defaultHeatConductivity;
double defaultHeatCapacity;
double defaultHeatConductivityScale;
double defaultHeatCapacityScale;
ScalerFunTwoArgs thermalConductivityScaling;
ScalerFunTwoArgs heatCapacityScaling;
double conductivity;
double capacity;
};
std::vector<BlockData> blockData;
std::vector<const CubitMeshSets *> meshset_vec_ptr,
Sev sev) {
for (
auto m : meshset_vec_ptr) {
std::vector<double> block_data;
CHKERR m->getAttributes(block_data);
if (block_data.size() < 2) {
"Expected that block has at least two attributes");
}
auto get_block_ents = [&]() {
m_field.
get_moab().get_entities_by_handle(
m->meshset, ents,
true);
return ents;
};
<<
m->getName() <<
": conductivity = " << block_data[0]
<< " capacity = " << block_data[1];
blockData.push_back({block_data[0], block_data[1], get_block_ents()});
}
}
boost::shared_ptr<double> conductivityPtr;
boost::shared_ptr<double> conductivityScalingPtr;
boost::shared_ptr<double> capacityPtr;
boost::shared_ptr<double> capacityScalingPtr;
};
pipeline.push_back(new OpMatThermalBlocks(
blockedParamsPtr->getHeatConductivityPtr(),
blockedParamsPtr->getHeatCapacityPtr(),
blockedParamsPtr->getHeatConductivityScalingPtr(),
blockedParamsPtr->getHeatCapacityScalingPtr(), local_heat_conductivity,
local_heat_capacity, local_heat_conductivity_scale,
local_heat_capacity_scale, thermal_conductivity_scaling_func,
heat_capacity_scaling_func, m_field, sev,
m_field.
getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
(boost::format("%s(.*)") % block_name).str()
))
));
}
struct BlockedThermoElasticParameters
: public boost::enable_shared_from_this<BlockedThermoElasticParameters> {
double refTemp;
inline auto getCoeffExpansionPtr() {
return boost::shared_ptr<VectorDouble>(shared_from_this(), &coeffExpansion);
}
inline auto getRefTempPtr() {
return boost::shared_ptr<double>(shared_from_this(), &refTemp);
}
};
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
std::string block_name,
boost::shared_ptr<BlockedThermoElasticParameters> blockedParamsPtr,
OpMatThermoElasticBlocks(boost::shared_ptr<VectorDouble> expansion_ptr,
boost::shared_ptr<double> ref_temp_ptr,
double local_coeff_expansion,
Sev sev,
std::vector<const CubitMeshSets *> meshset_vec_ptr)
expansionPtr(expansion_ptr), refTempPtr(ref_temp_ptr),
defaultCoeffExpansion(local_coeff_expansion),
defaultRefTemp(local_ref_temp) {
extractThermoElasticBlockData(m_field, meshset_vec_ptr, sev),
"Cannot get data from thermoelastic block");
}
EntitiesFieldData::EntData &data) {
for (auto &b : blockData) {
if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
*expansionPtr = b.expansion;
*refTempPtr = b.ref_temp;
}
}
*refTempPtr = defaultRefTemp;
}
private:
double defaultCoeffExpansion;
double defaultRefTemp;
double ref_temp;
};
std::vector<BlockData> blockData;
std::vector<const CubitMeshSets *> meshset_vec_ptr, Sev sev) {
for (
auto m : meshset_vec_ptr) {
std::vector<double> block_data;
CHKERR m->getAttributes(block_data);
if (block_data.size() < 2) {
"Expected that block has at least two attributes");
}
auto get_block_ents = [&]() {
m_field.
get_moab().get_entities_by_handle(
m->meshset, ents,
true);
return ents;
};
auto get_expansion = [&]() {
if (block_data.size() > 2) {
expansion[1] = block_data[2];
}
if (
SPACE_DIM == 3 && block_data.size() > 3) {
expansion[2] = block_data[3];
}
return expansion;
};
auto coeff_exp_vec = get_expansion();
<< " ref_temp = " << block_data[0]
<< " expansion = " << coeff_exp_vec;
blockData.push_back({block_data[0], coeff_exp_vec, get_block_ents()});
}
}
boost::shared_ptr<VectorDouble> expansionPtr;
boost::shared_ptr<double> refTempPtr;
};
pipeline.push_back(new OpMatThermoElasticBlocks(
blockedParamsPtr->getCoeffExpansionPtr(),
blockedParamsPtr->getRefTempPtr(), local_coeff_expansion, local_ref_temp,
m_field, sev,
m_field.
getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
(boost::format("%s(.*)") % block_name).str()
))
));
}
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<MatrixDouble> mDptr,
boost::shared_ptr<VectorDouble> coeff_expansion_ptr);
EntitiesFieldData::EntData &col_data);
private:
boost::shared_ptr<MatrixDouble>
mDPtr;
};
boost::shared_ptr<MatrixDouble> strain_ptr,
boost::shared_ptr<VectorDouble> temp_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr,
boost::shared_ptr<VectorDouble> coeff_expansion_ptr,
boost::shared_ptr<MatrixDouble> stress_ptr);
boost::shared_ptr<VectorDouble> temp_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr,
boost::shared_ptr<VectorDouble> coeff_expansion_ptr,
boost::shared_ptr<double> ref_temp_ptr,
boost::shared_ptr<MatrixDouble> stress_ptr);
private:
boost::shared_ptr<VectorDouble>
tempPtr;
boost::shared_ptr<MatrixDouble>
mDPtr;
};
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<MatrixDouble> mDptr,
boost::shared_ptr<VectorDouble> coeff_expansion_ptr)
mDPtr(mDptr), coeffExpansionPtr(coeff_expansion_ptr) {
sYmm = false;
}
OpKCauchyThermoElasticity::iNtegrate(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
const auto nb_integration_pts = row_data.getN().size1();
const auto nb_row_base_functions = row_data.getN().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_row_diff_base = row_data.getFTensor1DiffN<
SPACE_DIM>();
auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mDPtr);
t_coeff_exp(d, d) = (*coeffExpansionPtr)[
d];
}
t_eigen_strain(
i,
j) = (t_D(
i,
j,
k,
l) * t_coeff_exp(
k,
l));
for (auto gg = 0; gg != nb_integration_pts; ++gg) {
double alpha = getMeasure() * t_w;
auto rr = 0;
auto t_mat =
DataLayoutTraits<DataLayout::CoeffsByGauss>>(
auto t_col_base = col_data.getFTensor0N(gg, 0);
(t_row_diff_base(
j) * t_eigen_strain(
i,
j)) * (t_col_base * alpha);
++t_mat;
++t_col_base;
}
++t_row_diff_base;
}
for (; rr != nb_row_base_functions; ++rr)
++t_row_diff_base;
++t_w;
}
}
OpStressThermal::OpStressThermal(
const std::string
field_name, boost::shared_ptr<MatrixDouble> strain_ptr,
boost::shared_ptr<VectorDouble> temp_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr,
boost::shared_ptr<VectorDouble> coeff_expansion_ptr,
boost::shared_ptr<MatrixDouble> stress_ptr)
tempPtr(temp_ptr), mDPtr(m_D_ptr), coeffExpansionPtr(coeff_expansion_ptr),
stressPtr(stress_ptr) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
OpStressThermal::OpStressThermal(
boost::shared_ptr<MatrixDouble> strain_ptr,
boost::shared_ptr<VectorDouble> temp_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr,
boost::shared_ptr<VectorDouble> coeff_expansion_ptr,
boost::shared_ptr<double> ref_temp_ptr,
boost::shared_ptr<MatrixDouble> stress_ptr)
tempPtr(temp_ptr), mDPtr(m_D_ptr), coeffExpansionPtr(coeff_expansion_ptr),
refTempPtr(ref_temp_ptr), stressPtr(stress_ptr) {}
const auto nb_gauss_pts = getGaussPts().size2();
auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mDPtr);
auto t_strain = getFTensor2SymmetricFromMat<SPACE_DIM>(*strainPtr);
auto t_stress = getFTensor2SymmetricFromMat<SPACE_DIM>(*stressPtr);
auto t_temp = getFTensor0FromVec(*tempPtr);
t_coeff_exp(d, d) = (*coeffExpansionPtr)[
d];
}
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
(t_strain(
k,
l) - t_coeff_exp(
k,
l) * (t_temp - (*refTempPtr)));
++t_strain;
++t_stress;
++t_temp;
}
}
struct SetTargetTemperature;
}
template <AssemblyType A, IntegrationType I,
typename OpBase>
template <AssemblyType A, IntegrationType I, typename OpBase>
OpBase>;
template <AssemblyType A, IntegrationType I, typename OpBase>
MoFEM::OpFluxRhsImpl<ThermoElasticOps::SetTargetTemperature, 1, 1, A, I,
OpBase>,
A, I, OpBase
> {
AddFluxToRhsPipelineImpl() = delete;
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
boost::shared_ptr<VectorDouble> temp_ptr, std::string block_name,
Sev sev
) {
using OP_SOURCE = typename FormsIntegrators<OpBase>::template Assembly<
using OP_TEMP = typename FormsIntegrators<OpBase>::template Assembly<
auto add_op = [&](auto &&meshset_vec_ptr) {
for (
auto m : meshset_vec_ptr) {
std::vector<double> block_data;
m->getAttributes(block_data);
if (block_data.size() < 2) {
"Expected two parameters");
}
double target_temperature = block_data[0];
double beta =
block_data[1];
auto ents_ptr = boost::make_shared<Range>();
*(ents_ptr), true);
<<
"Add " << *
m <<
" target temperature " << target_temperature
<< " penalty " << beta;
pipeline.push_back(new OP_SOURCE(
[target_temperature, beta](double, double, double) {
return target_temperature * beta;
},
ents_ptr));
pipeline.push_back(new OP_TEMP(
[beta](double, double, double) { return -beta; }, ents_ptr));
}
};
m_field.
getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
(boost::format("%s(.*)") % block_name).str()
))
);
}
};
template <AssemblyType A, IntegrationType I,
typename OpBase>
struct AddFluxToLhsPipelineImpl<
> {
AddFluxToLhsPipelineImpl() = delete;
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
boost::shared_ptr<VectorDouble> temp_ptr, std::string block_name, Sev sev
) {
using OP_MASS = typename FormsIntegrators<OpBase>::template Assembly<
auto add_op = [&](auto &&meshset_vec_ptr) {
for (
auto m : meshset_vec_ptr) {
std::vector<double> block_data;
m->getAttributes(block_data);
if (block_data.size() != 2) {
"Expected two parameters");
}
double beta =
block_data[1];
auto ents_ptr = boost::make_shared<Range>();
*(ents_ptr), true);
<<
"Add " << *
m <<
" penalty " << beta;
pipeline.push_back(new OP_MASS(
[beta](double, double, double) { return -beta; }, ents_ptr));
}
};
m_field.
getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
(boost::format("%s(.*)") % block_name).str()
))
);
}
};
template <AssemblyType A, IntegrationType I>
boost::shared_ptr<VectorDouble> dot_T_ptr,
boost::shared_ptr<VectorDouble> T_ptr,
boost::shared_ptr<MatrixDouble> grad_T_ptr,
boost::shared_ptr<double> initial_T_ptr,
boost::shared_ptr<double> peak_T_ptr)
dotTPtr(dot_T_ptr), TPtr(T_ptr), gradTPtr(grad_T_ptr),
initialTPtr(initial_T_ptr), peakTPtr(peak_T_ptr) {}
const double vol = getMeasure();
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_base = data.getFTensor0N();
auto t_diff_base = data.getFTensor1DiffN<
SPACE_DIM>();
#ifndef NDEBUG
if (data.getDiffN().size1() != data.getN().size1())
if (data.getDiffN().size2() != data.getN().size2() *
SPACE_DIM) {
<< "Side " << rowSide << " " << CN::EntityTypeName(rowType);
MOFEM_LOG(
"SELF", Sev::error) << data.getN();
MOFEM_LOG(
"SELF", Sev::error) << data.getDiffN();
}
#endif
const double alpha = t_w * vol;
const double set_T = init_T(*initialTPtr, *peakTPtr, t_coords(0),
t_coords(1), t_coords(2));
int bb = 0;
locF[bb] += (t_base * alpha) * (t_T - set_T);
++t_base;
++t_diff_base;
}
++t_base;
++t_diff_base;
}
++t_T;
++t_coords;
++t_w;
}
}
private:
boost::shared_ptr<VectorDouble> dotTPtr;
boost::shared_ptr<VectorDouble> TPtr;
boost::shared_ptr<MatrixDouble> gradTPtr;
boost::shared_ptr<MatrixDouble> gradQPtr;
boost::shared_ptr<double> initialTPtr;
boost::shared_ptr<double> peakTPtr;
};
template <AssemblyType A, IntegrationType I>
boost::shared_ptr<VectorDouble> T_ptr)
TPtr(T_ptr) {
sYmm = false;
}
EntitiesFieldData::EntData &col_data) {
const double vol = getMeasure();
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_row_base = row_data.getFTensor0N();
auto t_row_diff_base = row_data.getFTensor1DiffN<
SPACE_DIM>();
#ifndef NDEBUG
if (row_data.getDiffN().size1() != row_data.getN().size1())
if (row_data.getDiffN().size2() != row_data.getN().size2() *
SPACE_DIM) {
<< "Side " << rowSide << " " << CN::EntityTypeName(rowType);
MOFEM_LOG(
"SELF", Sev::error) << row_data.getN();
MOFEM_LOG(
"SELF", Sev::error) << row_data.getDiffN();
}
if (col_data.getDiffN().size1() != col_data.getN().size1())
if (col_data.getDiffN().size2() != col_data.getN().size2() *
SPACE_DIM) {
<< "Side " << rowSide << " " << CN::EntityTypeName(rowType);
MOFEM_LOG(
"SELF", Sev::error) << col_data.getN();
MOFEM_LOG(
"SELF", Sev::error) << col_data.getDiffN();
}
#endif
const double alpha = t_w * vol;
int rr = 0;
auto t_col_base = col_data.getFTensor0N(gg, 0);
auto t_col_diff_base = col_data.getFTensor1DiffN<
SPACE_DIM>(gg, 0);
for (
int cc = 0; cc !=
nbCols; ++cc) {
locMat(rr, cc) += (t_row_base * t_col_base * alpha);
++t_col_base;
++t_col_diff_base;
}
++t_row_base;
++t_row_diff_base;
}
++t_row_base;
++t_row_diff_base;
}
++t_T;
++t_w;
++t_coords;
}
}
private:
boost::shared_ptr<VectorDouble> TPtr;
};
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
auto getFTensor1FromMat(M &data, int rr=0, int cc=0)
Get tensor rank 1 (vector) form data matrix.
MoFEM::LogManager::SeverityLevel Sev
static auto getFTensor0FromVec(V &data)
Get tensor rank 0 (scalar) form data vector.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
int nbRows
number of dofs on rows
int nbIntegrationPts
number of integration points
MatrixDouble locMat
local entity block matrix
virtual MoFEMErrorCode iNtegrate(EntData &, EntData &)
Integrate grad-grad operator.
int nbCols
number if dof on column
int nbRowBaseFunctions
number or row base functions
boost::shared_ptr< MatrixDouble > mDPtr
OpKCauchyThermoElasticity(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< MatrixDouble > mDptr, boost::shared_ptr< VectorDouble > coeff_expansion_ptr)
boost::shared_ptr< VectorDouble > coeffExpansionPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
DEPRECATED OpStressThermal(const std::string field_name, boost::shared_ptr< MatrixDouble > strain_ptr, boost::shared_ptr< VectorDouble > temp_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr, boost::shared_ptr< VectorDouble > coeff_expansion_ptr, boost::shared_ptr< MatrixDouble > stress_ptr)
boost::shared_ptr< VectorDouble > coeffExpansionPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate stress]
boost::shared_ptr< MatrixDouble > strainPtr
boost::shared_ptr< MatrixDouble > stressPtr
boost::shared_ptr< double > refTempPtr
boost::shared_ptr< VectorDouble > tempPtr
boost::shared_ptr< MatrixDouble > mDPtr
