v0.14.0
ThermoElasticOps.hpp
/** \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);
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;
}
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;
}
}
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]
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) *
++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
> {
AddFluxToRhsPipelineImpl() = delete;
static MoFEMErrorCode add(
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));
}
};
CHKERR add_op(
m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
(boost::format("%s(.*)") % block_name).str()
))
);
}
};
template <AssemblyType A, IntegrationType I, typename OpBase>
struct AddFluxToLhsPipelineImpl<
OpFluxLhsImpl<ThermoElasticOps::SetTargetTemperature, 1, 1, A, I, OpBase>,
> {
AddFluxToLhsPipelineImpl() = delete;
static MoFEMErrorCode add(
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));
}
};
CHKERR add_op(
m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
(boost::format("%s(.*)") % block_name).str()
))
);
}
};
//! [Target temperature]
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
ThermoElasticOps::OpStressThermal::mDPtr
boost::shared_ptr< MatrixDouble > mDPtr
Definition: ThermoElasticOps.hpp:45
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
DEPRECATED
#define DEPRECATED
Definition: definitions.h:17
MoFEM::OpBaseImpl::locMat
MatrixDouble locMat
local entity block matrix
Definition: FormsIntegrators.hpp:249
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
ThermoElasticOps::OpStressThermal::coeffExpansionPtr
boost::shared_ptr< double > coeffExpansionPtr
Definition: ThermoElasticOps.hpp:46
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
ref_temp
double ref_temp
Definition: thermo_elastic.cpp:124
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
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
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::OpFluxLhsImpl
Definition: Natural.hpp:43
I
constexpr IntegrationType I
Definition: operators_tests.cpp:31
ThermoElasticOps::OpStressThermal::strainPtr
boost::shared_ptr< MatrixDouble > strainPtr
Definition: ThermoElasticOps.hpp:43
MoFEM::Sev
MoFEM::LogManager::SeverityLevel Sev
Definition: CoreTemplates.hpp:17
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
ThermoElasticOps::OpKCauchyThermoElasticity::mDPtr
boost::shared_ptr< MatrixDouble > mDPtr
Definition: ThermoElasticOps.hpp:19
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
DomainEleOp
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
Definition: child_and_parent.cpp:36
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: seepage.cpp:57
convert.type
type
Definition: convert.py:64
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
ThermoElasticOps::OpStressThermal::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate stress]
Definition: ThermoElasticOps.hpp:132
MoFEM::AddFluxToRhsPipelineImpl
Definition: Natural.hpp:46
ThermoElasticOps
Definition: ThermoElasticOps.hpp:7
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
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 >
MoFEM::OpBaseImpl::nbCols
int nbCols
number if dof on column
Definition: FormsIntegrators.hpp:237
DomainEleOp
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
MoFEM::OpBaseImpl::nbRows
int nbRows
number of dofs on rows
Definition: FormsIntegrators.hpp:236
ThermoElasticOps::OpStressThermal::stressPtr
boost::shared_ptr< MatrixDouble > stressPtr
Definition: ThermoElasticOps.hpp:47
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_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
ThermoElasticOps::OpKCauchyThermoElasticity::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ThermoElasticOps.hpp:61
ThermoElasticOps::OpStressThermal::tempPtr
boost::shared_ptr< VectorDouble > tempPtr
Definition: ThermoElasticOps.hpp:44
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
AssemblyDomainEleOp
FormsIntegrators< DomainEleOp >::Assembly< A >::OpBase AssemblyDomainEleOp
Definition: tensor_divergence_operator.cpp:59
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21