v0.15.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<VectorDouble> 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<VectorDouble> coeffExpansionPtr;
21};
22
24 /**
25 * @deprecated do not use this constructor
26 */
28 OpStressThermal(const std::string field_name,
29 boost::shared_ptr<MatrixDouble> strain_ptr,
30 boost::shared_ptr<VectorDouble> temp_ptr,
31 boost::shared_ptr<MatrixDouble> m_D_ptr,
32 boost::shared_ptr<VectorDouble> coeff_expansion_ptr,
33 boost::shared_ptr<MatrixDouble> stress_ptr);
34
35 OpStressThermal(boost::shared_ptr<MatrixDouble> strain_ptr,
36 boost::shared_ptr<VectorDouble> temp_ptr,
37 boost::shared_ptr<MatrixDouble> m_D_ptr,
38 boost::shared_ptr<VectorDouble> coeff_expansion_ptr,
39 boost::shared_ptr<double> ref_temp_ptr,
40 boost::shared_ptr<MatrixDouble> stress_ptr);
41
42 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
43
44private:
45 boost::shared_ptr<MatrixDouble> strainPtr;
46 boost::shared_ptr<VectorDouble> tempPtr;
47 boost::shared_ptr<MatrixDouble> mDPtr;
48 boost::shared_ptr<VectorDouble> coeffExpansionPtr;
49 boost::shared_ptr<double> refTempPtr;
50 boost::shared_ptr<MatrixDouble> stressPtr;
51};
52
54 const std::string row_field_name, const std::string col_field_name,
55 boost::shared_ptr<MatrixDouble> mDptr,
56 boost::shared_ptr<VectorDouble> coeff_expansion_ptr)
57 : AssemblyDomainEleOp(row_field_name, col_field_name,
58 DomainEleOp::OPROWCOL),
59 mDPtr(mDptr), coeffExpansionPtr(coeff_expansion_ptr) {
60 sYmm = false;
61}
62
63MoFEMErrorCode
64OpKCauchyThermoElasticity::iNtegrate(EntitiesFieldData::EntData &row_data,
65 EntitiesFieldData::EntData &col_data) {
67
68 auto &locMat = AssemblyDomainEleOp::locMat;
69
70 const auto nb_integration_pts = row_data.getN().size1();
71 const auto nb_row_base_functions = row_data.getN().size2();
72 auto t_w = getFTensor0IntegrationWeight();
73 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
74 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mDPtr);
75
81
83 t_coeff_exp(i, j) = 0;
84 for (auto d = 0; d != SPACE_DIM; ++d) {
85 t_coeff_exp(d, d) = (*coeffExpansionPtr)[d];
86 }
87
88 t_eigen_strain(i, j) = (t_D(i, j, k, l) * t_coeff_exp(k, l));
89
90 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
91
92 double alpha = getMeasure() * t_w;
93 auto rr = 0;
94 for (; rr != AssemblyDomainEleOp::nbRows / SPACE_DIM; ++rr) {
95 auto t_mat = getFTensor1FromMat<SPACE_DIM, 1>(locMat, rr * SPACE_DIM);
96 auto t_col_base = col_data.getFTensor0N(gg, 0);
97 for (auto cc = 0; cc != AssemblyDomainEleOp::nbCols; cc++) {
98
99 t_mat(i) -=
100 (t_row_diff_base(j) * t_eigen_strain(i, j)) * (t_col_base * alpha);
101
102 ++t_mat;
103 ++t_col_base;
104 }
105
106 ++t_row_diff_base;
107 }
108 for (; rr != nb_row_base_functions; ++rr)
109 ++t_row_diff_base;
110
111 ++t_w;
112 }
113
115}
116
118 const std::string field_name, boost::shared_ptr<MatrixDouble> strain_ptr,
119 boost::shared_ptr<VectorDouble> temp_ptr,
120 boost::shared_ptr<MatrixDouble> m_D_ptr,
121 boost::shared_ptr<VectorDouble> coeff_expansion_ptr,
122 boost::shared_ptr<MatrixDouble> stress_ptr)
123 : DomainEleOp(field_name, DomainEleOp::OPROW), strainPtr(strain_ptr),
124 tempPtr(temp_ptr), mDPtr(m_D_ptr), coeffExpansionPtr(coeff_expansion_ptr),
125 stressPtr(stress_ptr) {
126 // Operator is only executed for vertices
127 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
128}
129
131 boost::shared_ptr<MatrixDouble> strain_ptr,
132 boost::shared_ptr<VectorDouble> temp_ptr,
133 boost::shared_ptr<MatrixDouble> m_D_ptr,
134 boost::shared_ptr<VectorDouble> coeff_expansion_ptr,
135 boost::shared_ptr<double> ref_temp_ptr,
136 boost::shared_ptr<MatrixDouble> stress_ptr)
137 : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), strainPtr(strain_ptr),
138 tempPtr(temp_ptr), mDPtr(m_D_ptr), coeffExpansionPtr(coeff_expansion_ptr),
139 refTempPtr(ref_temp_ptr), stressPtr(stress_ptr) {}
140
141//! [Calculate stress]
142MoFEMErrorCode OpStressThermal::doWork(int side, EntityType type,
143 EntData &data) {
145 const auto nb_gauss_pts = getGaussPts().size2();
146 stressPtr->resize((SPACE_DIM * (SPACE_DIM + 1)) / 2, nb_gauss_pts);
147
148 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mDPtr);
149 auto t_strain = getFTensor2SymmetricFromMat<SPACE_DIM>(*strainPtr);
150 auto t_stress = getFTensor2SymmetricFromMat<SPACE_DIM>(*stressPtr);
151 auto t_temp = getFTensor0FromVec(*tempPtr);
152
157
159 t_coeff_exp(i, j) = 0;
160 for (auto d = 0; d != SPACE_DIM; ++d) {
161 t_coeff_exp(d, d) = (*coeffExpansionPtr)[d];
162 }
163
164 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
165 t_stress(i, j) =
166 t_D(i, j, k, l) *
167 (t_strain(k, l) - t_coeff_exp(k, l) * (t_temp - (*refTempPtr)));
168 ++t_strain;
169 ++t_stress;
170 ++t_temp;
171 }
173}
174//! [Calculate stress]
175
176struct SetTargetTemperature;
177
178} // namespace ThermoElasticOps
179
180//! [Target temperature]
181
182template <AssemblyType A, IntegrationType I, typename OpBase>
183struct MoFEM::OpFluxRhsImpl<ThermoElasticOps::SetTargetTemperature, 1, 1, A, I,
184 OpBase>;
185
186template <AssemblyType A, IntegrationType I, typename OpBase>
187struct MoFEM::OpFluxLhsImpl<ThermoElasticOps::SetTargetTemperature, 1, 1, A, I,
188 OpBase>;
189
190template <AssemblyType A, IntegrationType I, typename OpBase>
192
193 MoFEM::OpFluxRhsImpl<ThermoElasticOps::SetTargetTemperature, 1, 1, A, I,
194 OpBase>,
195 A, I, OpBase
196
197 > {
198
200
202
203 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
204 MoFEM::Interface &m_field, const std::string field_name,
205 boost::shared_ptr<VectorDouble> temp_ptr, std::string block_name, Sev sev
206
207 ) {
209
210 using OP_SOURCE = typename FormsIntegrators<OpBase>::template Assembly<
211 A>::template LinearForm<I>::template OpSource<1, 1>;
212 using OP_TEMP = typename FormsIntegrators<OpBase>::template Assembly<
213 A>::template LinearForm<I>::template OpBaseTimesScalar<1>;
214
215 auto add_op = [&](auto &&meshset_vec_ptr) {
217 for (auto m : meshset_vec_ptr) {
218 std::vector<double> block_data;
219 m->getAttributes(block_data);
220 if (block_data.size() < 2) {
221 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
222 "Expected two parameters");
223 }
224 double target_temperature = block_data[0];
225 double beta =
226 block_data[1]; // Set temperature parameter [ W/K * (1/m^3)]
227 auto ents_ptr = boost::make_shared<Range>();
228 CHKERR m_field.get_moab().get_entities_by_handle(m->meshset,
229 *(ents_ptr), true);
230
231 MOFEM_TAG_AND_LOG("WORLD", sev, "SetTargetTemperature")
232 << "Add " << *m << " target temperature " << target_temperature
233 << " penalty " << beta;
234
235 pipeline.push_back(new OP_SOURCE(
237 [target_temperature, beta](double, double, double) {
238 return target_temperature * beta;
239 },
240 ents_ptr));
241 pipeline.push_back(new OP_TEMP(
242 field_name, temp_ptr,
243 [beta](double, double, double) { return -beta; }, ents_ptr));
244 }
246 };
247
248 CHKERR add_op(
249
250 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
251
252 (boost::format("%s(.*)") % block_name).str()
253
254 ))
255
256 );
257
258 MOFEM_LOG_CHANNEL("WORLD");
259
261 }
262};
263
264template <AssemblyType A, IntegrationType I, typename OpBase>
265struct AddFluxToLhsPipelineImpl<
266
267 OpFluxLhsImpl<ThermoElasticOps::SetTargetTemperature, 1, 1, A, I, OpBase>,
268 A, I, OpBase
269
270 > {
271
273
274 static MoFEMErrorCode add(
275
276 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
277 MoFEM::Interface &m_field, const std::string field_name,
278 boost::shared_ptr<VectorDouble> temp_ptr, std::string block_name, Sev sev
279
280 ) {
282
283 using OP_MASS = typename FormsIntegrators<OpBase>::template Assembly<
285
286 auto add_op = [&](auto &&meshset_vec_ptr) {
288 for (auto m : meshset_vec_ptr) {
289 std::vector<double> block_data;
290 m->getAttributes(block_data);
291 if (block_data.size() != 2) {
292 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
293 "Expected two parameters");
294 }
295 double beta =
296 block_data[1]; // Set temperature parameter [ W/K * (1/m^3)]
297 auto ents_ptr = boost::make_shared<Range>();
298 CHKERR m_field.get_moab().get_entities_by_handle(m->meshset,
299 *(ents_ptr), true);
300
301 MOFEM_TAG_AND_LOG("WORLD", sev, "SetTargetTemperature")
302 << "Add " << *m << " penalty " << beta;
303
304 pipeline.push_back(new OP_MASS(
306 [beta](double, double, double) { return -beta; }, ents_ptr));
307 }
309 };
310
311 CHKERR add_op(
312
313 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
314
315 (boost::format("%s(.*)") % block_name).str()
316
317 ))
318
319 );
320
321 MOFEM_LOG_CHANNEL("WORLD");
322
324 }
325};
326
327//! [Target temperature]
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
constexpr int SPACE_DIM
EntitiesFieldData::EntData EntData
Data on entities.
@ NOSPACE
Definition definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define DEPRECATED
Definition definitions.h:17
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
SeverityLevel
Severity levels.
#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
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
constexpr IntegrationType I
constexpr auto field_name
FTensor::Index< 'm', 3 > m
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
Deprecated interface functions.
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
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