v0.14.0
Note
Prerequisites of this tutorial include VEC-0: Linear elasticity

Note
Intended learning outcome:
• multi-physics
• multi-field problem
/**
* \file thermo_elastic.cpp
* \example thermo_elastic.cpp
*
* Thermo plasticity
*
*/
#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 3
#endif
#include <MoFEM.hpp>
using namespace MoFEM;
constexpr int SPACE_DIM =
EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
using BoundaryEle =
//! [Linear elastic problem]
0>; //< Elastic stiffness matrix
SPACE_DIM>; //< Elastic internal forces
//! [Linear elastic problem]
//! [Thermal problem]
/**
* @brief Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)
*
*/
/**
* @brief Integrate Lhs div of base of flux time base of temperature (FLUX x T)
* and transpose of it, i.e. (T x FLAX)
*
*/
GAUSS>::OpMixDivTimesScalar<SPACE_DIM>;
/**
* @brief Integrate Lhs base of temperature times (heat capacity) times base of
* temperature (T x T)
*
*/
/**
* @brief Integrating Rhs flux base (1/k) flux (FLUX)
*/
GAUSS>::OpBaseTimesVector<3, SPACE_DIM, 1>;
/**
* @brief Integrate Rhs div flux base times temperature (T)
*
*/
GAUSS>::OpMixDivTimesU<3, 1, SPACE_DIM>;
/**
* @brief Integrate Rhs base of temperature time heat capacity times heat rate
* (T)
*
*/
GAUSS>::OpBaseTimesScalar<1>;
/**
* @brief Integrate Rhs base of temperature times divergence of flux (T)
*
*/
//! [Thermal problem]
//! [Body and heat source]
using OpBodyForce =
using OpHeatSource =
//! [Body and heat source]
//! [Natural boundary conditions]
//! [Natural boundary conditions]
//! [Essential boundary conditions (Least square approach)]
GAUSS>::OpEssentialRhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
GAUSS>::OpEssentialLhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
//! [Essential boundary conditions (Least square approach)]
double default_poisson_ratio = 0.25;
double ref_temp = 0.0;
double init_temp = 0.0;
1; // Force / (time temperature ) or Power /
// (length temperature) // Time unit is hour. force unit kN
double default_heat_capacity = 1; // length^2/(time^2 temperature) // length is
// millimeter time is hour
int order = 2; //< default approximation order
#include <ThermoElasticOps.hpp> //< additional coupling opearyors
using namespace ThermoElasticOps; //< name space of coupling operators
auto save_range = [](moab::Interface &moab, const std::string name,
const Range r) {
auto out_meshset = get_temp_meshset_ptr(moab);
CHKERR moab.write_file(name.c_str(), "VTK", "", out_meshset->get_ptr(), 1);
};
ThermoElasticProblem(MoFEM::Interface &m_field) : mField(m_field) {}
MoFEMErrorCode runProblem();
private:
PetscBool doEvalField;
std::array<double, SPACE_DIM> fieldEvalCoords;
boost::shared_ptr<FieldEvaluatorInterface::SetPtsData> fieldEvalData;
boost::shared_ptr<VectorDouble> scalarFieldPtr;
boost::shared_ptr<MatrixDouble> vectorFieldPtr;
boost::shared_ptr<MatrixDouble> tensorFieldPtr;
MoFEMErrorCode createCommonData(); //< read global data from command line
MoFEMErrorCode bC(); //< read boundary conditions
MoFEMErrorCode OPs(); //< add operators to pipeline
MoFEMErrorCode tsSolve(); //< time solver
struct BlockedParameters
: public boost::enable_shared_from_this<BlockedParameters> {
double coeffExpansion;
double heatConductivity;
double heatCapacity;
inline auto getDPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &mD);
}
inline auto getCoeffExpansionPtr() {
return boost::shared_ptr<double>(shared_from_this(), &coeffExpansion);
}
inline auto getHeatConductivityPtr() {
return boost::shared_ptr<double>(shared_from_this(), &heatConductivity);
}
inline auto getHeatCapacityPtr() {
return boost::shared_ptr<double>(shared_from_this(), &heatCapacity);
}
};
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
std::string block_elastic_name, std::string block_thermal_name,
boost::shared_ptr<BlockedParameters> blockedParamsPtr, Sev sev);
};
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
std::string block_elastic_name, std::string block_thermal_name,
boost::shared_ptr<BlockedParameters> blockedParamsPtr, Sev sev) {
struct OpMatElasticBlocks : public DomainEleOp {
OpMatElasticBlocks(boost::shared_ptr<MatrixDouble> m, double bulk_modulus_K,
Sev sev,
std::vector<const CubitMeshSets *> meshset_vec_ptr)
: DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), matDPtr(m),
bulkModulusKDefault(bulk_modulus_K),
shearModulusGDefault(shear_modulus_G) {
CHK_THROW_MESSAGE(extractElasticBlockData(m_field, meshset_vec_ptr, sev),
"Can not get data from block");
}
MoFEMErrorCode doWork(int side, EntityType type,
for (auto &b : blockData) {
if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
CHKERR getMatDPtr(matDPtr, b.bulkModulusK, b.shearModulusG);
}
}
CHKERR getMatDPtr(matDPtr, bulkModulusKDefault, shearModulusGDefault);
}
private:
boost::shared_ptr<MatrixDouble> matDPtr;
struct BlockData {
double bulkModulusK;
double shearModulusG;
Range blockEnts;
};
double bulkModulusKDefault;
double shearModulusGDefault;
std::vector<BlockData> blockData;
extractElasticBlockData(MoFEM::Interface &m_field,
std::vector<const CubitMeshSets *> meshset_vec_ptr,
Sev sev) {
for (auto m : meshset_vec_ptr) {
MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Elastic Block") << *m;
std::vector<double> block_data;
CHKERR m->getAttributes(block_data);
if (block_data.size() < 2) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Expected that block has two attributes");
}
auto get_block_ents = [&]() {
Range ents;
m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
return ents;
};
double young_modulus = block_data[0];
double poisson_ratio = block_data[1];
double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Elastic Block")
<< m->getName() << ": E = " << young_modulus
<< " nu = " << poisson_ratio;
blockData.push_back(
{bulk_modulus_K, shear_modulus_G, get_block_ents()});
}
}
MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
double bulk_modulus_K, double shear_modulus_G) {
//! [Calculate elasticity tensor]
auto set_material_stiffness = [&]() {
double A = (SPACE_DIM == 2)
: 1;
auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
t_D(i, j, k, l) =
2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
t_kd(k, l);
};
//! [Calculate elasticity tensor]
constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
mat_D_ptr->resize(size_symm * size_symm, 1);
set_material_stiffness();
}
};
double default_bulk_modulus_K =
double default_shear_modulus_G =
pipeline.push_back(new OpMatElasticBlocks(
blockedParamsPtr->getDPtr(), default_bulk_modulus_K,
default_bulk_modulus_K, mField, sev,
// Get blockset using regular expression
mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
(boost::format("%s(.*)") % block_elastic_name).str()
))
));
struct OpMatThermalBlocks : public DomainEleOp {
OpMatThermalBlocks(boost::shared_ptr<double> expansion_ptr,
boost::shared_ptr<double> conductivity_ptr,
boost::shared_ptr<double> capacity_ptr,
MoFEM::Interface &m_field, Sev sev,
std::vector<const CubitMeshSets *> meshset_vec_ptr)
expansionPtr(expansion_ptr), conductivityPtr(conductivity_ptr),
capacityPtr(capacity_ptr) {
CHK_THROW_MESSAGE(extractThermallockData(m_field, meshset_vec_ptr, sev),
"Can not get data from block");
}
MoFEMErrorCode doWork(int side, EntityType type,
for (auto &b : blockData) {
if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
*expansionPtr = b.expansion;
*conductivityPtr = b.conductivity;
*capacityPtr = b.capacity;
}
}
*expansionPtr = default_coeff_expansion;
*conductivityPtr = default_heat_conductivity;
*capacityPtr = default_heat_capacity;
}
private:
struct BlockData {
double expansion;
double conductivity;
double capacity;
Range blockEnts;
};
std::vector<BlockData> blockData;
extractThermallockData(MoFEM::Interface &m_field,
std::vector<const CubitMeshSets *> meshset_vec_ptr,
Sev sev) {
for (auto m : meshset_vec_ptr) {
MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Thermal Block") << *m;
std::vector<double> block_data;
CHKERR m->getAttributes(block_data);
if (block_data.size() < 3) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Expected that block has two attributes");
}
auto get_block_ents = [&]() {
Range ents;
m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
return ents;
};
MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Thermal Block")
<< m->getName() << ": expansion = " << block_data[0]
<< " conductivity = " << block_data[1] << " capacity "
<< block_data[2];
blockData.push_back(
{block_data[0], block_data[1], block_data[2], get_block_ents()});
}
}
boost::shared_ptr<double> expansionPtr;
boost::shared_ptr<double> conductivityPtr;
boost::shared_ptr<double> capacityPtr;
};
pipeline.push_back(new OpMatThermalBlocks(
blockedParamsPtr->getCoeffExpansionPtr(),
blockedParamsPtr->getHeatConductivityPtr(),
blockedParamsPtr->getHeatCapacityPtr(), mField, sev,
// Get blockset using regular expression
mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
(boost::format("%s(.*)") % block_thermal_name).str()
))
));
}
//! [Run problem]
CHKERR setupProblem();
CHKERR createCommonData();
CHKERR bC();
CHKERR OPs();
CHKERR tsSolve();
}
//! [Run problem]
//! [Set up problem]
Simple *simple = mField.getInterface<Simple>();
// Mechanical fields
// Temperature
constexpr auto flux_space = (SPACE_DIM == 2) ? HCURL : HDIV;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
CHKERR simple->setFieldOrder("U", order + 1);
CHKERR simple->setFieldOrder("FLUX", order + 1);
CHKERR simple->setFieldOrder("T", order);
CHKERR simple->setUp();
int coords_dim = SPACE_DIM;
CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
fieldEvalCoords.data(), &coords_dim,
scalarFieldPtr = boost::make_shared<VectorDouble>();
vectorFieldPtr = boost::make_shared<MatrixDouble>();
tensorFieldPtr = boost::make_shared<MatrixDouble>();
if (doEvalField) {
fieldEvalData =
mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
if constexpr (SPACE_DIM == 3) {
CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree3D(
fieldEvalData, simple->getDomainFEName());
} else {
CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree2D(
fieldEvalData, simple->getDomainFEName());
}
fieldEvalData->setEvalPoints(fieldEvalCoords.data(), 1);
auto no_rule = [](int, int, int) { return -1; };
auto field_eval_fe_ptr = fieldEvalData->feMethodPtr.lock();
field_eval_fe_ptr->getRuleHook = no_rule;
field_eval_fe_ptr->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("T", scalarFieldPtr));
field_eval_fe_ptr->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>("U", vectorFieldPtr));
field_eval_fe_ptr->getOpPtrVector().push_back(
"U", tensorFieldPtr));
}
}
//! [Set up problem]
//! [Create common data]
auto get_command_line_parameters = [&]() {
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
&default_young_modulus, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
&default_poisson_ratio, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-coeff_expansion",
&default_coeff_expansion, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-ref_temp", &ref_temp,
PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-init_temp", &init_temp,
PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-capacity",
&default_heat_capacity, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-conductivity",
MOFEM_LOG("ThermoElastic", Sev::inform)
<< "Young modulus " << default_young_modulus;
MOFEM_LOG("ThermoElastic", Sev::inform)
<< "Poisson ratio " << default_poisson_ratio;
MOFEM_LOG("ThermoElastic", Sev::inform)
<< "Coeff of expansion " << default_coeff_expansion;
MOFEM_LOG("ThermoElastic", Sev::inform)
<< "Capacity " << default_heat_capacity;
MOFEM_LOG("ThermoElastic", Sev::inform)
<< "Heat conductivity " << default_heat_conductivity;
MOFEM_LOG("ThermoElastic", Sev::inform)
<< "Reference temperature " << ref_temp;
MOFEM_LOG("ThermoElastic", Sev::inform)
<< "Initial temperature " << init_temp;
};
CHKERR get_command_line_parameters();
}
//! [Create common data]
//! [Boundary condition]
MOFEM_LOG("SYNC", Sev::noisy) << "bC";
MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
auto simple = mField.getInterface<Simple>();
auto bc_mng = mField.getInterface<BcManager>();
CHKERR bc_mng->removeBlockDOFsOnEntities<DisplacementCubitBcData>(
simple->getProblemName(), "U");
CHKERR bc_mng->pushMarkDOFsOnEntities<HeatFluxCubitBcData>(
simple->getProblemName(), "FLUX", false);
auto get_skin = [&]() {
Range body_ents;
CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
Skinner skin(&mField.get_moab());
Range skin_ents;
CHKERR skin.find_skin(0, body_ents, false, skin_ents);
return skin_ents;
};
auto filter_flux_blocks = [&](auto skin) {
auto remove_cubit_blocks = [&](auto c) {
for (auto m :
mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(c)
) {
Range ents;
CHKERR mField.get_moab().get_entities_by_dimension(
m->getMeshset(), SPACE_DIM - 1, ents, true);
skin = subtract(skin, ents);
}
};
auto remove_named_blocks = [&](auto n) {
for (auto m : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
std::regex(
(boost::format("%s(.*)") % n).str()
))
) {
Range ents;
CHKERR mField.get_moab().get_entities_by_dimension(
m->getMeshset(), SPACE_DIM - 1, ents, true);
skin = subtract(skin, ents);
}
};
CHK_THROW_MESSAGE(remove_cubit_blocks(NODESET | TEMPERATURESET),
"remove_cubit_blocks");
CHK_THROW_MESSAGE(remove_cubit_blocks(SIDESET | HEATFLUXSET),
"remove_cubit_blocks");
CHK_THROW_MESSAGE(remove_named_blocks("TEMPERATURE"),
"remove_named_blocks");
CHK_THROW_MESSAGE(remove_named_blocks("HEATFLUX"), "remove_named_blocks");
return skin;
};
auto filter_true_skin = [&](auto skin) {
Range boundary_ents;
ParallelComm *pcomm =
ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
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()));
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
remove_flux_ents);
MOFEM_LOG("SYNC", Sev::noisy) << remove_flux_ents << endl;
MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
#ifndef NDEBUG
mField.get_moab(),
(boost::format("flux_remove_%d.vtk") % mField.get_comm_rank()).str(),
remove_flux_ents);
#endif
CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
simple->getProblemName(), "FLUX", remove_flux_ents);
auto set_init_temp = [](boost::shared_ptr<FieldEntity> field_entity_ptr) {
field_entity_ptr->getEntFieldData()[0] = init_temp;
return 0;
};
CHKERR mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(set_init_temp,
"T");
}
//! [Boundary condition]
//! [Push operators to pipeline]
MOFEM_LOG("SYNC", Sev::noisy) << "OPs";
MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
auto pipeline_mng = mField.getInterface<PipelineManager>();
auto simple = mField.getInterface<Simple>();
auto bc_mng = mField.getInterface<BcManager>();
auto boundary_marker =
bc_mng->getMergedBlocksMarker(vector<string>{"HEATFLUX"});
auto integration_rule = [](int, int, int approx_order) {
return 2 * approx_order;
};
CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
auto block_params = boost::make_shared<BlockedParameters>();
auto mDPtr = block_params->getDPtr();
auto coeff_expansion_ptr = block_params->getCoeffExpansionPtr();
auto heat_conductivity_ptr = block_params->getHeatConductivityPtr();
auto heat_capacity_ptr = block_params->getHeatCapacityPtr();
// Default time scaling of BCs and sources, from command line arguments
auto time_scale = boost::make_shared<TimeScale>();
// Files which define scaling for separate variables, if provided
auto time_bodyforce_scaling =
boost::make_shared<TimeScale>("bodyforce_scale.txt");
auto time_heatsource_scaling =
boost::make_shared<TimeScale>("heatsource_scale.txt");
auto time_temperature_scaling =
boost::make_shared<TimeScale>("temperature_bc_scale.txt");
auto time_displacement_scaling =
boost::make_shared<TimeScale>("displacement_bc_scale.txt");
auto time_flux_scaling = boost::make_shared<TimeScale>("flux_bc_scale.txt");
auto time_force_scaling =
boost::make_shared<TimeScale>("force_bc_scale.txt");
auto add_domain_rhs_ops = [&](auto &pipeline) {
Sev::inform);
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(new OpCalculateScalarFieldValues("T", vec_temp_ptr));
pipeline.push_back(
new OpCalculateScalarFieldValuesDot("T", vec_temp_dot_ptr));
"FLUX", vec_temp_div_ptr));
pipeline.push_back(
new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
pipeline.push_back(
pipeline.push_back(new OpStressThermal(mat_strain_ptr, vec_temp_ptr, mDPtr,
coeff_expansion_ptr,
mat_stress_ptr));
pipeline.push_back(new OpSetBc("FLUX", true, boundary_marker));
pipeline.push_back(new OpInternalForceCauchy(
"U", mat_stress_ptr,
[](double, double, double) constexpr { return 1; }));
auto resistance = [heat_conductivity_ptr](const double, const double,
const double) {
return (1. / (*heat_conductivity_ptr));
};
// negative value is consequence of symmetric system
auto capacity = [heat_capacity_ptr](const double, const double,
const double) {
return -(*heat_capacity_ptr);
};
auto unity = [](const double, const double, const double) constexpr {
return -1.;
};
pipeline.push_back(new OpHdivFlux("FLUX", mat_flux_ptr, resistance));
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.push_back(new OpUnSetBc("FLUX"));
// Set source terms
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);
pipeline.push_back(new OpSetBc("FLUX", true, boundary_marker));
pipeline.push_back(new OpKCauchy("U", "U", mDPtr));
"U", "T", mDPtr, coeff_expansion_ptr));
auto resistance = [heat_conductivity_ptr](const double, const double,
const double) {
return (1. / (*heat_conductivity_ptr));
};
auto capacity = [heat_capacity_ptr](const double, const double,
const double) {
return -(*heat_capacity_ptr);
};
pipeline.push_back(new OpHdivHdiv("FLUX", "FLUX", resistance));
pipeline.push_back(new OpHdivT(
"FLUX", "T", []() constexpr { return -1; }, true));
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.push_back(new OpUnSetBc("FLUX"));
auto vec_temp_ptr = boost::make_shared<VectorDouble>();
pipeline.push_back(new OpCalculateScalarFieldValues("T", vec_temp_ptr));
pipeline, mField, "T", vec_temp_ptr, "SETTEMP", Sev::verbose);
};
auto add_boundary_rhs_ops = [&](auto &pipeline) {
pipeline.push_back(new OpSetBc("FLUX", true, boundary_marker));
// Set BCs using the least squares imposition
pipeline_mng->getOpBoundaryRhsPipeline(), mField, "U",
{time_scale, time_force_scaling}, "FORCE", Sev::inform);
pipeline_mng->getOpBoundaryRhsPipeline(), mField, "FLUX",
{time_scale, time_temperature_scaling}, "TEMPERATURE", Sev::inform);
pipeline.push_back(new OpUnSetBc("FLUX"));
auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(
new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
mField, pipeline, simple->getProblemName(), "FLUX", mat_flux_ptr,
{time_scale, time_flux_scaling});
};
auto add_boundary_lhs_ops = [&](auto &pipeline) {
mField, pipeline, simple->getProblemName(), "FLUX");
};
// Set BCs by eliminating degrees of freedom
auto get_bc_hook_rhs = [&]() {
mField, pipeline_mng->getDomainRhsFE(),
{time_scale, time_displacement_scaling});
return hook;
};
auto get_bc_hook_lhs = [&]() {
mField, pipeline_mng->getDomainLhsFE(),
{time_scale, time_displacement_scaling});
return hook;
};
pipeline_mng->getDomainRhsFE()->preProcessHook = get_bc_hook_rhs();
pipeline_mng->getDomainLhsFE()->preProcessHook = get_bc_hook_lhs();
}
//! [Push operators to pipeline]
//! [Solve]
Simple *simple = mField.getInterface<Simple>();
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
ISManager *is_manager = mField.getInterface<ISManager>();
auto dm = simple->getDM();
auto solver = pipeline_mng->createTSIM();
auto snes_ctx_ptr = getDMSnesCtx(dm);
auto set_section_monitor = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
CHKERR SNESMonitorSet(snes,
(MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
(void *)(snes_ctx_ptr.get()), nullptr);
};
auto create_post_process_element = [&]() {
auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
auto block_params = boost::make_shared<BlockedParameters>();
auto mDPtr = block_params->getDPtr();
auto coeff_expansion_ptr = block_params->getCoeffExpansionPtr();
"MAT_THERMAL", block_params, Sev::verbose);
post_proc_fe->getOpPtrVector(), {H1, HDIV});
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>();
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("T", vec_temp_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
auto u_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpStressThermal(mat_strain_ptr, vec_temp_ptr, mDPtr,
coeff_expansion_ptr, mat_stress_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"T", vec_temp_ptr}},
{{"U", u_ptr}, {"FLUX", mat_flux_ptr}},
{},
{{"STRAIN", mat_strain_ptr}, {"STRESS", mat_stress_ptr}}
)
);
return post_proc_fe;
};
auto monitor_ptr = boost::make_shared<FEMethod>();
auto post_proc_fe = create_post_process_element();
auto set_time_monitor = [&](auto dm, auto solver) {
monitor_ptr->preProcessHook = [&]() {
CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
post_proc_fe,
monitor_ptr->getCacheWeakPtr());
CHKERR post_proc_fe->writeFile(
"out_" + boost::lexical_cast<std::string>(monitor_ptr->ts_step) +
".h5m");
if (doEvalField) {
if constexpr (SPACE_DIM == 3) {
CHKERR mField.getInterface<FieldEvaluatorInterface>()
->evalFEAtThePoint3D(
fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
simple->getDomainFEName(), fieldEvalData,
mField.get_comm_rank(), mField.get_comm_rank(), nullptr,
} else {
CHKERR mField.getInterface<FieldEvaluatorInterface>()
->evalFEAtThePoint2D(
fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
simple->getDomainFEName(), fieldEvalData,
mField.get_comm_rank(), mField.get_comm_rank(), nullptr,
}
if (scalarFieldPtr->size()) {
auto t_temp = getFTensor0FromVec(*scalarFieldPtr);
MOFEM_LOG("ThermoElasticSync", Sev::inform)
<< "Eval point T: " << t_temp;
}
if (vectorFieldPtr->size1()) {
auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vectorFieldPtr);
MOFEM_LOG("ThermoElasticSync", Sev::inform)
<< "Eval point U magnitude: " << sqrt(t_disp(i) * t_disp(i));
}
if (tensorFieldPtr->size1()) {
getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*tensorFieldPtr);
MOFEM_LOG("ThermoElasticSync", Sev::inform)
}
MOFEM_LOG_SYNCHRONISE(mField.get_comm());
}
};
auto null = boost::shared_ptr<FEMethod>();
CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(), null,
monitor_ptr, null);
};
auto set_fieldsplit_preconditioner = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
PC pc;
CHKERR KSPGetPC(ksp, &pc);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
// Setup fieldsplit (block) solver - optional: yes/no
if (is_pcfs == PETSC_TRUE) {
auto bc_mng = mField.getInterface<BcManager>();
auto is_mng = mField.getInterface<ISManager>();
auto name_prb = simple->getProblemName();
CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "U", 0,
SPACE_DIM, is_u);
CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "FLUX", 0, 0,
is_flux);
CHKERR is_mng->isCreateProblemFieldAndRank(name_prb, ROW, "T", 0, 0,
is_T);
IS is_tmp;
CHKERR ISExpand(is_T, is_flux, &is_tmp);
auto is_TFlux = SmartPetscObj<IS>(is_tmp);
auto is_all_bc = bc_mng->getBlockIS(name_prb, "HEATFLUX", "FLUX", 0, 0);
int is_all_bc_size;
CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
MOFEM_LOG("ThermoElastic", Sev::inform)
<< "Field split block size " << is_all_bc_size;
if (is_all_bc_size) {
IS is_tmp2;
CHKERR ISDifference(is_TFlux, is_all_bc, &is_tmp2);
is_TFlux = SmartPetscObj<IS>(is_tmp2);
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
is_all_bc); // boundary block
}
CHKERR ISSort(is_u);
CHKERR ISSort(is_TFlux);
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_TFlux);
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
}
};
auto D = createDMVector(dm);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR TSSetSolution(solver, D);
CHKERR TSSetFromOptions(solver);
CHKERR set_section_monitor(solver);
CHKERR set_fieldsplit_preconditioner(solver);
CHKERR set_time_monitor(dm, solver);
CHKERR TSSetUp(solver);
CHKERR TSSolve(solver, NULL);
}
//! [Solve]
static char help[] = "...\n\n";
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
// Add logging channel for example
auto core_log = logging::core::get();
LogManager::setLog("ThermoElastic");
MOFEM_LOG_TAG("ThermoElastic", "ThermoElastic");
LogManager::setLog("ThermoElasticSync");
MOFEM_LOG_TAG("ThermoElasticSync", "ThermoElasticSync");
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database interface
//! [Create MoFEM]
CHKERR simple->getOptions();
//! [ThermoElasticProblem]
ThermoElasticProblem ex(m_field);
CHKERR ex.runProblem();
//! [ThermoElasticProblem]
}
}
/** \file ThermoElasticOps.hpp
* \example ThermoElasticOps.hpp
*/
namespace ThermoElasticOps {
struct OpKCauchyThermoElasticity : public AssemblyDomainEleOp {
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<MatrixDouble> mDptr,
boost::shared_ptr<double> coeff_expansion_ptr);
private:
boost::shared_ptr<MatrixDouble> mDPtr;
boost::shared_ptr<double> coeffExpansionPtr;
};
struct OpStressThermal : public DomainEleOp {
/**
* @deprecated do not use this constructor
*/
boost::shared_ptr<MatrixDouble> strain_ptr,
boost::shared_ptr<VectorDouble> temp_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr,
boost::shared_ptr<double> coeff_expansion_ptr,
boost::shared_ptr<MatrixDouble> stress_ptr);
DEPRECATED OpStressThermal(boost::shared_ptr<MatrixDouble> strain_ptr,
boost::shared_ptr<VectorDouble> temp_ptr,
boost::shared_ptr<MatrixDouble> m_D_ptr,
boost::shared_ptr<double> coeff_expansion_ptr,
boost::shared_ptr<MatrixDouble> stress_ptr);
MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
private:
boost::shared_ptr<MatrixDouble> strainPtr;
boost::shared_ptr<VectorDouble> tempPtr;
boost::shared_ptr<MatrixDouble> mDPtr;
boost::shared_ptr<double> coeffExpansionPtr;
boost::shared_ptr<MatrixDouble> stressPtr;
};
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<MatrixDouble> mDptr,
boost::shared_ptr<double> coeff_expansion_ptr)
: AssemblyDomainEleOp(row_field_name, col_field_name,
DomainEleOp::OPROWCOL),
mDPtr(mDptr), coeffExpansionPtr(coeff_expansion_ptr) {
sYmm = false;
}
OpKCauchyThermoElasticity::iNtegrate(EntitiesFieldData::EntData &row_data,
const auto nb_integration_pts = row_data.getN().size1();
const auto nb_row_base_functions = row_data.getN().size2();
auto t_w = getFTensor0IntegrationWeight();
constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mDPtr);
t_eigen_strain(i, j) = (t_D(i, j, k, l) * t_kd(k, l)) * (*coeffExpansionPtr);
for (auto gg = 0; gg != nb_integration_pts; ++gg) {
double alpha = getMeasure() * t_w;
auto rr = 0;
for (; rr != AssemblyDomainEleOp::nbRows / SPACE_DIM; ++rr) {
auto t_mat = getFTensor1FromMat<SPACE_DIM, 1>(locMat, rr * SPACE_DIM);
auto t_col_base = col_data.getFTensor0N(gg, 0);
for (auto cc = 0; cc != AssemblyDomainEleOp::nbCols; cc++) {
t_mat(i) -=
(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<double> coeff_expansion_ptr,
boost::shared_ptr<MatrixDouble> stress_ptr)
: DomainEleOp(field_name, DomainEleOp::OPROW), strainPtr(strain_ptr),
tempPtr(temp_ptr), mDPtr(m_D_ptr), coeffExpansionPtr(coeff_expansion_ptr),
stressPtr(stress_ptr) {
// Operator is only executed for vertices
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<double> coeff_expansion_ptr,
boost::shared_ptr<MatrixDouble> stress_ptr)
: DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), strainPtr(strain_ptr),
tempPtr(temp_ptr), mDPtr(m_D_ptr), coeffExpansionPtr(coeff_expansion_ptr),
stressPtr(stress_ptr) {}
//! [Calculate stress]
MoFEMErrorCode OpStressThermal::doWork(int side, EntityType type,
EntData &data) {
const auto nb_gauss_pts = getGaussPts().size2();
stressPtr->resize((SPACE_DIM * (SPACE_DIM + 1)) / 2, nb_gauss_pts);
constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
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);
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
t_stress(i, j) =
t_D(i, j, k, l) * (t_strain(k, l) - t_kd(k, l) * (t_temp - ref_temp) *
(*coeffExpansionPtr));
++t_strain;
++t_stress;
++t_temp;
}
}
//! [Calculate stress]
struct SetTargetTemperature;
} // namespace ThermoElasticOps
//! [Target temperature]
template <AssemblyType A, IntegrationType I, typename OpBase>
struct MoFEM::OpFluxRhsImpl<ThermoElasticOps::SetTargetTemperature, 1, 1, A, I,
template <AssemblyType A, IntegrationType I, typename OpBase>
struct MoFEM::OpFluxLhsImpl<ThermoElasticOps::SetTargetTemperature, 1, 1, A, I,
OpBase>;
template <AssemblyType A, IntegrationType I, typename OpBase>
MoFEM::OpFluxRhsImpl<ThermoElasticOps::SetTargetTemperature, 1, 1, A, I,
OpBase>,
A, I, OpBase
> {
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
MoFEM::Interface &m_field, const std::string field_name,
boost::shared_ptr<VectorDouble> temp_ptr, std::string block_name, Sev sev
) {
using OP_SOURCE = typename FormsIntegrators<OpBase>::template Assembly<
A>::template LinearForm<I>::template OpSource<1, 1>;
using OP_TEMP = typename FormsIntegrators<OpBase>::template Assembly<
A>::template LinearForm<I>::template OpBaseTimesScalar<1>;
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) {
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"Expected two parameters");
}
double target_temperature = block_data[0];
double beta =
block_data[1]; // Set temperature parameter [ W/K * (1/m^3)]
auto ents_ptr = boost::make_shared<Range>();
CHKERR m_field.get_moab().get_entities_by_handle(m->meshset,
*(ents_ptr), true);
MOFEM_TAG_AND_LOG("WORLD", sev, "SetTargetTemperature")
<< "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(
field_name, temp_ptr,
[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>
OpFluxLhsImpl<ThermoElasticOps::SetTargetTemperature, 1, 1, A, I, OpBase>,
> {
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
MoFEM::Interface &m_field, const std::string field_name,
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) {
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"Expected two parameters");
}
double beta =
block_data[1]; // Set temperature parameter [ W/K * (1/m^3)]
auto ents_ptr = boost::make_shared<Range>();
CHKERR m_field.get_moab().get_entities_by_handle(m->meshset,
*(ents_ptr), true);
MOFEM_TAG_AND_LOG("WORLD", sev, "SetTargetTemperature")
<< "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()
))
);
}
};
//! [Target temperature]
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
ThermoElasticOps::OpStressThermal::mDPtr
boost::shared_ptr< MatrixDouble > mDPtr
Definition: ThermoElasticOps.hpp:45
default_coeff_expansion
double default_coeff_expansion
Definition: thermo_elastic.cpp:127
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
SIDESET
@ SIDESET
Definition: definitions.h:147
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
ThermoElasticProblem::tsSolve
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
Definition: thermo_elastic.cpp:885
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
DEPRECATED
#define DEPRECATED
Definition: definitions.h:17
MOFEM_LOG_SEVERITY_SYNC
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
Definition: LogManager.hpp:352
MoFEM::OpBaseImpl::locMat
MatrixDouble locMat
local entity block matrix
Definition: FormsIntegrators.hpp:239
OpHdivHdiv
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, SPACE_DIM > OpHdivHdiv
[Linear elastic problem]
Definition: thermo_elastic.cpp:47
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:172
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:596
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
help
static char help[]
Definition: activate_deactivate_dofs.cpp:13
OpHDivTemp
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 1, SPACE_DIM > OpHDivTemp
Integrate Rhs div flux base times temperature (T)
Definition: thermo_elastic.cpp:76
ThermoElasticOps::OpStressThermal::coeffExpansionPtr
boost::shared_ptr< double > coeffExpansionPtr
Definition: ThermoElasticOps.hpp:46
MoFEM::EssentialBC
Essential boundary conditions.
Definition: Essential.hpp:111
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
ref_temp
double ref_temp
Definition: thermo_elastic.cpp:124
ThermoElasticProblem::bC
MoFEMErrorCode bC()
[Create common data]
Definition: thermo_elastic.cpp:556
MoFEM::OpCalculateScalarFieldValuesDot
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
Definition: UserDataOperators.hpp:273
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:104
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
ThermoElasticProblem::setupProblem
MoFEMErrorCode setupProblem()
Definition: thermo_elastic.cpp:448
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
bulk_modulus_K
double bulk_modulus_K
Definition: dynamic_first_order_con_law.cpp:96
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:527
MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
default_poisson_ratio
double default_poisson_ratio
Definition: thermo_elastic.cpp:123
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
OpBaseDotT
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBaseDotT
Integrate Rhs base of temperature time heat capacity times heat rate (T)
Definition: thermo_elastic.cpp:84
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
ThermoElasticOps::OpStressThermal::OpStressThermal
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< double > coeff_expansion_ptr, boost::shared_ptr< MatrixDouble > stress_ptr)
Definition: ThermoElasticOps.hpp:109
sdf.r
int r
Definition: sdf.py:8
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::OpFluxLhsImpl
Definition: Natural.hpp:43
ROW
@ ROW
Definition: definitions.h:123
I
constexpr IntegrationType I
Definition: operators_tests.cpp:31
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
approx_order
static constexpr int approx_order
Definition: prism_polynomial_approximation.cpp:14
ThermoElasticOps::OpStressThermal::strainPtr
boost::shared_ptr< MatrixDouble > strainPtr
Definition: ThermoElasticOps.hpp:43
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
NODESET
@ NODESET
Definition: definitions.h:146
MoFEM::Sev
MoFEM::LogManager::SeverityLevel Sev
Definition: CoreTemplates.hpp:17
ThermoElasticOps.hpp
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:170
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1018
ThermoElasticOps::OpKCauchyThermoElasticity::mDPtr
boost::shared_ptr< MatrixDouble > mDPtr
Definition: ThermoElasticOps.hpp:19
MoFEM::EssentialPreProc< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:91
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
default_heat_conductivity
double default_heat_conductivity
Definition: thermo_elastic.cpp:128
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: seepage.cpp:57
ThermoElasticProblem
Definition: thermo_elastic.cpp:153
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
MoFEM::OpCalculateHVecVectorField
Get vector field for H-div approximation.
Definition: UserDataOperators.hpp:2116
MoFEM::get_temp_meshset_ptr
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1857
MoFEM::FieldEvaluatorInterface
Field evaluator interface.
Definition: FieldEvaluator.hpp:21
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
double
MoFEM::OpEssentialLhsImpl
Enforce essential constrains on lhs.
Definition: Essential.hpp:81
convert.type
type
Definition: convert.py:64
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:302
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
MoFEM::LogManager::getStrmSync
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
Definition: LogManager.cpp:348
OpBaseDivFlux
OpBaseDotT OpBaseDivFlux
Integrate Rhs base of temperature times divergence of flux (T)
Definition: thermo_elastic.cpp:90
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:173
ThermoElasticOps::OpStressThermal::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate stress]
Definition: ThermoElasticOps.hpp:132
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:128
OpKCauchy
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Linear elastic problem]
Definition: thermo_elastic.cpp:34
Definition: Natural.hpp:49
doEvalField
PetscBool doEvalField
Definition: incompressible_elasticity.cpp:41
ThermoElasticProblem::createCommonData
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: thermo_elastic.cpp:509
init_temp
double init_temp
Definition: thermo_elastic.cpp:125
MoFEM::LogManager::SeverityLevel
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
Definition: Natural.hpp:46
ThermoElasticOps
Definition: ThermoElasticOps.hpp:7
size_symm
constexpr auto size_symm
Definition: plastic.cpp:42
ThermoElasticOps::OpKCauchyThermoElasticity::coeffExpansionPtr
boost::shared_ptr< double > coeffExpansionPtr
Definition: ThermoElasticOps.hpp:20
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
BiLinearForm
main
int main(int argc, char *argv[])
Definition: activate_deactivate_dofs.cpp:15
OpInternalForceCauchy
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: thermo_elastic.cpp:38
ThermoElasticProblem::OPs
MoFEMErrorCode OPs()
[Boundary condition]
Definition: thermo_elastic.cpp:665
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', SPACE_DIM >
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:503
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1536
convert.n
n
Definition: convert.py:82
integration_rule
auto integration_rule
Definition: free_surface.cpp:185
MoFEMErrorCode addMatBlockOps(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string block_elastic_name, std::string block_thermal_name, boost::shared_ptr< BlockedParameters > blockedParamsPtr, Sev sev)
Definition: thermo_elastic.cpp:206
MoFEM::OpUnSetBc
Definition: FormsIntegrators.hpp:49
MoFEM::OpBaseImpl::nbCols
int nbCols
number if dof on column
Definition: FormsIntegrators.hpp:227
Range
DomainEleOp
MoFEM::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition: DeprecatedPetsc.hpp:192
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Definition: LogManager.hpp:362
save_range
auto save_range
Definition: thermo_elastic.cpp:144
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::NaturalBC::Assembly
Assembly methods.
Definition: Natural.hpp:65
MoFEM::OpEssentialRhsImpl
Enforce essential constrains on rhs.
Definition: Essential.hpp:65
MoFEM::OpBaseImpl::nbRows
int nbRows
number of dofs on rows
Definition: FormsIntegrators.hpp:226
shear_modulus_G
double shear_modulus_G
Definition: dynamic_first_order_con_law.cpp:97
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
default_young_modulus
double default_young_modulus
[Essential boundary conditions (Least square approach)]
Definition: thermo_elastic.cpp:122
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1094
ThermoElasticOps::OpStressThermal::stressPtr
boost::shared_ptr< MatrixDouble > stressPtr
Definition: ThermoElasticOps.hpp:47
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
ThermoElasticOps::OpKCauchyThermoElasticity::OpKCauchyThermoElasticity
OpKCauchyThermoElasticity(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< MatrixDouble > mDptr, boost::shared_ptr< double > coeff_expansion_ptr)
Definition: ThermoElasticOps.hpp:50
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::DMMoFEMTSSetMonitor
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMoFEM.cpp:1060
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::OpSetBc
Set indices on entities on finite element.
Definition: FormsIntegrators.hpp:38
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
TEMPERATURESET
@ TEMPERATURESET
Definition: definitions.h:155
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::BlockData
Definition: MeshsetsManager.cpp:752
MoFEM::OpCalculateHdivVectorDivergence
Calculate divergence of vector field.
Definition: UserDataOperators.hpp:2215
ThermoElasticOps::OpKCauchyThermoElasticity::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ThermoElasticOps.hpp:61
MoFEM::MoFEMSNESMonitorFields
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition: SnesCtx.cpp:232
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
MoFEM::EssentialBC::Assembly
Assembly methods.
Definition: Essential.hpp:119
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
OpCapacity
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)
Definition: thermo_elastic.cpp:63
ThermoElasticOps::OpStressThermal::tempPtr
boost::shared_ptr< VectorDouble > tempPtr
Definition: ThermoElasticOps.hpp:44
MoFEM::HeatFluxCubitBcData
Definition of the heat flux bc data structure.
Definition: BCData.hpp:427
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
default_heat_capacity
double default_heat_capacity
Definition: thermo_elastic.cpp:131
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
MoFEM::getDMSnesCtx
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:1046
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
HEATFLUXSET
@ HEATFLUXSET
Definition: definitions.h:156
QUIET
@ QUIET
Definition: definitions.h:208
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEM::MeshsetsManager::getCubitMeshsetPtr
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
Definition: MeshsetsManager.cpp:575
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
MoFEM::SmartPetscObj< IS >
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
OpHdivT
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< SPACE_DIM > OpHdivT
Integrate Lhs div of base of flux time base of temperature (FLUX x T) and transpose of it,...
Definition: thermo_elastic.cpp:55
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:590
MF_EXIST
@ MF_EXIST
Definition: definitions.h:100
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
ThermoElasticProblem::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: thermo_elastic.cpp:436
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
OpHdivFlux
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 3, SPACE_DIM, 1 > OpHdivFlux
Integrating Rhs flux base (1/k) flux (FLUX)
Definition: thermo_elastic.cpp:69