v0.14.0
Loading...
Searching...
No Matches
ThermoElasticOps.hpp
Go to the documentation of this file.
1
2
3/** \file ThermoElasticOps.hpp
4 * \example ThermoElasticOps.hpp
5 */
6
7namespace ThermoElasticOps {
8
11 const std::string row_field_name, const std::string col_field_name,
12 boost::shared_ptr<MatrixDouble> mDptr,
13 boost::shared_ptr<double> coeff_expansion_ptr);
14
15 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
16 EntitiesFieldData::EntData &col_data);
17
18private:
19 boost::shared_ptr<MatrixDouble> mDPtr;
20 boost::shared_ptr<double> coeffExpansionPtr;
21};
22
24 OpStressThermal(const std::string field_name,
25 boost::shared_ptr<MatrixDouble> strain_ptr,
26 boost::shared_ptr<VectorDouble> temp_ptr,
27 boost::shared_ptr<MatrixDouble> m_D_ptr,
28 boost::shared_ptr<double> coeff_expansion_ptr,
29 boost::shared_ptr<MatrixDouble> stress_ptr);
30 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
31
32private:
33 boost::shared_ptr<MatrixDouble> strainPtr;
34 boost::shared_ptr<VectorDouble> tempPtr;
35 boost::shared_ptr<MatrixDouble> mDPtr;
36 boost::shared_ptr<double> coeffExpansionPtr;
37 boost::shared_ptr<MatrixDouble> stressPtr;
38};
39
41 const std::string row_field_name, const std::string col_field_name,
42 boost::shared_ptr<MatrixDouble> mDptr,
43 boost::shared_ptr<double> coeff_expansion_ptr)
44 : AssemblyDomainEleOp(row_field_name, col_field_name,
45 DomainEleOp::OPROWCOL),
46 mDPtr(mDptr), coeffExpansionPtr(coeff_expansion_ptr) {
47 sYmm = false;
48}
49
50MoFEMErrorCode
51OpKCauchyThermoElasticity::iNtegrate(EntitiesFieldData::EntData &row_data,
52 EntitiesFieldData::EntData &col_data) {
54
56
57 const auto nb_integration_pts = row_data.getN().size1();
58 const auto nb_row_base_functions = row_data.getN().size2();
59 auto t_w = getFTensor0IntegrationWeight();
60
61 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
62 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
63 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mDPtr);
64
70 t_eigen_strain(i, j) = (t_D(i, j, k, l) * t_kd(k, l)) * (*coeffExpansionPtr);
71
72 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
73
74 double alpha = getMeasure() * t_w;
75 auto rr = 0;
76 for (; rr != AssemblyDomainEleOp::nbRows / SPACE_DIM; ++rr) {
77 auto t_mat = getFTensor1FromMat<SPACE_DIM, 1>(locMat, rr * SPACE_DIM);
78 auto t_col_base = col_data.getFTensor0N(gg, 0);
79 for (auto cc = 0; cc != AssemblyDomainEleOp::nbCols; cc++) {
80
81 t_mat(i) -=
82 (t_row_diff_base(j) * t_eigen_strain(i, j)) * (t_col_base * alpha);
83
84 ++t_mat;
85 ++t_col_base;
86 }
87
88 ++t_row_diff_base;
89 }
90 for (; rr != nb_row_base_functions; ++rr)
91 ++t_row_diff_base;
92
93 ++t_w;
94 }
95
97}
98
100 boost::shared_ptr<MatrixDouble> strain_ptr,
101 boost::shared_ptr<VectorDouble> temp_ptr,
102 boost::shared_ptr<MatrixDouble> m_D_ptr,
103 boost::shared_ptr<double> coeff_expansion_ptr,
104 boost::shared_ptr<MatrixDouble> stress_ptr)
105 : DomainEleOp(field_name, DomainEleOp::OPROW), strainPtr(strain_ptr),
106 tempPtr(temp_ptr), mDPtr(m_D_ptr), coeffExpansionPtr(coeff_expansion_ptr),
107 stressPtr(stress_ptr) {
108 // Operator is only executed for vertices
109 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
110}
111
112//! [Calculate stress]
113MoFEMErrorCode OpStressThermal::doWork(int side, EntityType type,
114 EntData &data) {
116 const auto nb_gauss_pts = getGaussPts().size2();
117 stressPtr->resize((SPACE_DIM * (SPACE_DIM + 1)) / 2, nb_gauss_pts);
118
119 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
120 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mDPtr);
121 auto t_strain = getFTensor2SymmetricFromMat<SPACE_DIM>(*strainPtr);
122 auto t_stress = getFTensor2SymmetricFromMat<SPACE_DIM>(*stressPtr);
123 auto t_temp = getFTensor0FromVec(*tempPtr);
128 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
129 t_stress(i, j) =
130 t_D(i, j, k, l) * (t_strain(k, l) - t_kd(k, l) * (t_temp - ref_temp) *
132 ++t_strain;
133 ++t_stress;
134 ++t_temp;
135 }
137}
138//! [Calculate stress]
139
140struct SetTargetTemperature;
141
142} // namespace ThermoElasticOps
143
144//! [Target temperature]
145
146template <AssemblyType A, IntegrationType I, typename OpBase>
147struct MoFEM::OpFluxRhsImpl<ThermoElasticOps::SetTargetTemperature, 1, 1, A, I,
148 OpBase>;
149
150template <AssemblyType A, IntegrationType I, typename OpBase>
151struct MoFEM::OpFluxLhsImpl<ThermoElasticOps::SetTargetTemperature, 1, 1, A, I,
152 OpBase>;
153
154template <AssemblyType A, IntegrationType I, typename OpBase>
156
157 MoFEM::OpFluxRhsImpl<ThermoElasticOps::SetTargetTemperature, 1, 1, A, I,
158 OpBase>,
159 A, I, OpBase
160
161 > {
162
164
166
167 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
168 MoFEM::Interface &m_field, const std::string field_name,
169 boost::shared_ptr<VectorDouble> temp_ptr, std::string block_name, Sev sev
170
171 ) {
173
174 using OP_SOURCE = typename FormsIntegrators<OpBase>::template Assembly<
175 A>::template LinearForm<I>::template OpSource<1, 1>;
176 using OP_TEMP = typename FormsIntegrators<OpBase>::template Assembly<
177 A>::template LinearForm<I>::template OpBaseTimesScalar<1>;
178
179 auto add_op = [&](auto &&meshset_vec_ptr) {
181 for (auto m : meshset_vec_ptr) {
182 std::vector<double> block_data;
183 m->getAttributes(block_data);
184 if (block_data.size() != 2) {
185 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
186 "Expected two parameters");
187 }
188 double target_temperature = block_data[0];
189 double beta =
190 block_data[1]; // Set temperature parameter [ W/K * (1/m^3)]
191 auto ents_ptr = boost::make_shared<Range>();
192 CHKERR m_field.get_moab().get_entities_by_handle(m->meshset,
193 *(ents_ptr), true);
194
195 MOFEM_TAG_AND_LOG("WORLD", sev, "SetTargetTemperature")
196 << "Add " << *m << " target temperature " << target_temperature
197 << " penalty " << beta;
198
199 pipeline.push_back(new OP_SOURCE(
201 [target_temperature, beta](double, double, double) {
202 return target_temperature * beta;
203 },
204 ents_ptr));
205 pipeline.push_back(new OP_TEMP(
206 field_name, temp_ptr,
207 [beta](double, double, double) { return -beta; }, ents_ptr));
208 }
210 };
211
212 CHKERR add_op(
213
214 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
215
216 (boost::format("%s(.*)") % block_name).str()
217
218 ))
219
220 );
221
222 MOFEM_LOG_CHANNEL("WORLD");
223
225 }
226};
227
228template <AssemblyType A, IntegrationType I, typename OpBase>
229struct AddFluxToLhsPipelineImpl<
230
231 OpFluxLhsImpl<ThermoElasticOps::SetTargetTemperature, 1, 1, A, I, OpBase>,
232 A, I, OpBase
233
234 > {
235
237
238 static MoFEMErrorCode add(
239
240 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
241 MoFEM::Interface &m_field, const std::string field_name,
242 boost::shared_ptr<VectorDouble> temp_ptr, std::string block_name, Sev sev
243
244 ) {
246
247 using OP_MASS = typename FormsIntegrators<OpBase>::template Assembly<
249
250 auto add_op = [&](auto &&meshset_vec_ptr) {
252 for (auto m : meshset_vec_ptr) {
253 std::vector<double> block_data;
254 m->getAttributes(block_data);
255 if (block_data.size() != 2) {
256 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
257 "Expected two parameters");
258 }
259 double beta =
260 block_data[1]; // Set temperature parameter [ W/K * (1/m^3)]
261 auto ents_ptr = boost::make_shared<Range>();
262 CHKERR m_field.get_moab().get_entities_by_handle(m->meshset,
263 *(ents_ptr), true);
264
265 MOFEM_TAG_AND_LOG("WORLD", sev, "SetTargetTemperature")
266 << "Add " << *m << " penalty " << beta;
267
268 pipeline.push_back(new OP_MASS(
270 [beta](double, double, double) { return -beta; }, ents_ptr));
271 }
273 };
274
275 CHKERR add_op(
276
277 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
278
279 (boost::format("%s(.*)") % block_name).str()
280
281 ))
282
283 );
284
285 MOFEM_LOG_CHANNEL("WORLD");
286
288 }
289};
290
291//! [Target temperature]
292
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
constexpr int SPACE_DIM
Kronecker Delta class.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FTensor::Index< 'm', SPACE_DIM > m
constexpr auto t_kd
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
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
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
constexpr IntegrationType I
constexpr AssemblyType A
constexpr auto field_name
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)
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)
virtual moab::Interface & get_moab()=0
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Interface for managing meshsets containing materials and boundary conditions.
int nbRows
number of dofs on rows
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
boost::shared_ptr< double > coeffExpansionPtr
boost::shared_ptr< MatrixDouble > mDPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
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)
boost::shared_ptr< double > coeffExpansionPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate stress]
boost::shared_ptr< MatrixDouble > strainPtr
boost::shared_ptr< MatrixDouble > stressPtr
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)
boost::shared_ptr< VectorDouble > tempPtr
boost::shared_ptr< MatrixDouble > mDPtr
double ref_temp