v0.14.0
ThermoElasticOps.hpp
Go to the documentation of this file.
1 
2 
3 /** \file ThermoElasticOps.hpp
4  * \example ThermoElasticOps.hpp
5  */
6 
7 namespace ThermoElasticOps {
8 
9 struct OpKCauchyThermoElasticity : public AssemblyDomainEleOp {
10  OpKCauchyThermoElasticity(
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 
18 private:
19  boost::shared_ptr<MatrixDouble> mDPtr;
20  boost::shared_ptr<VectorDouble> coeffExpansionPtr;
21 };
22 
23 struct OpStressThermal : public DomainEleOp {
24  /**
25  * @deprecated do not use this constructor
26  */
27  DEPRECATED
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<MatrixDouble> stress_ptr);
40 
41  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
42 
43 private:
44  boost::shared_ptr<MatrixDouble> strainPtr;
45  boost::shared_ptr<VectorDouble> tempPtr;
46  boost::shared_ptr<MatrixDouble> mDPtr;
47  boost::shared_ptr<VectorDouble> coeffExpansionPtr;
48  boost::shared_ptr<MatrixDouble> stressPtr;
49 };
50 
51 OpKCauchyThermoElasticity::OpKCauchyThermoElasticity(
52  const std::string row_field_name, const std::string col_field_name,
53  boost::shared_ptr<MatrixDouble> mDptr,
54  boost::shared_ptr<VectorDouble> coeff_expansion_ptr)
55  : AssemblyDomainEleOp(row_field_name, col_field_name,
56  DomainEleOp::OPROWCOL),
57  mDPtr(mDptr), coeffExpansionPtr(coeff_expansion_ptr) {
58  sYmm = false;
59 }
60 
61 MoFEMErrorCode
62 OpKCauchyThermoElasticity::iNtegrate(EntitiesFieldData::EntData &row_data,
63  EntitiesFieldData::EntData &col_data) {
64  MoFEMFunctionBegin;
65 
66  auto &locMat = AssemblyDomainEleOp::locMat;
67 
68  const auto nb_integration_pts = row_data.getN().size1();
69  const auto nb_row_base_functions = row_data.getN().size2();
70  auto t_w = getFTensor0IntegrationWeight();
71  auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
72  auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mDPtr);
73 
74  FTensor::Index<'i', SPACE_DIM> i;
75  FTensor::Index<'j', SPACE_DIM> j;
76  FTensor::Index<'k', SPACE_DIM> k;
77  FTensor::Index<'l', SPACE_DIM> l;
78  FTensor::Tensor2_symmetric<double, SPACE_DIM> t_eigen_strain;
79 
80  auto t_coeff_exp = FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM>();
81  t_coeff_exp(i, j) = 0;
82  for (auto d = 0; d != SPACE_DIM; ++d) {
83  t_coeff_exp(d, d) = (*coeffExpansionPtr)[d];
84  }
85 
86  t_eigen_strain(i, j) = (t_D(i, j, k, l) * t_coeff_exp(k, l));
87 
88  for (auto gg = 0; gg != nb_integration_pts; ++gg) {
89 
90  double alpha = getMeasure() * t_w;
91  auto rr = 0;
92  for (; rr != AssemblyDomainEleOp::nbRows / SPACE_DIM; ++rr) {
93  auto t_mat = getFTensor1FromMat<SPACE_DIM, 1>(locMat, rr * SPACE_DIM);
94  auto t_col_base = col_data.getFTensor0N(gg, 0);
95  for (auto cc = 0; cc != AssemblyDomainEleOp::nbCols; cc++) {
96 
97  t_mat(i) -=
98  (t_row_diff_base(j) * t_eigen_strain(i, j)) * (t_col_base * alpha);
99 
100  ++t_mat;
101  ++t_col_base;
102  }
103 
104  ++t_row_diff_base;
105  }
106  for (; rr != nb_row_base_functions; ++rr)
107  ++t_row_diff_base;
108 
109  ++t_w;
110  }
111 
112  MoFEMFunctionReturn(0);
113 }
114 
115 OpStressThermal::OpStressThermal(
116  const std::string field_name, boost::shared_ptr<MatrixDouble> strain_ptr,
117  boost::shared_ptr<VectorDouble> temp_ptr,
118  boost::shared_ptr<MatrixDouble> m_D_ptr,
119  boost::shared_ptr<VectorDouble> coeff_expansion_ptr,
120  boost::shared_ptr<MatrixDouble> stress_ptr)
121  : DomainEleOp(field_name, DomainEleOp::OPROW), strainPtr(strain_ptr),
122  tempPtr(temp_ptr), mDPtr(m_D_ptr), coeffExpansionPtr(coeff_expansion_ptr),
123  stressPtr(stress_ptr) {
124  // Operator is only executed for vertices
125  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
126 }
127 
128 OpStressThermal::OpStressThermal(
129  boost::shared_ptr<MatrixDouble> strain_ptr,
130  boost::shared_ptr<VectorDouble> temp_ptr,
131  boost::shared_ptr<MatrixDouble> m_D_ptr,
132  boost::shared_ptr<VectorDouble> coeff_expansion_ptr,
133  boost::shared_ptr<MatrixDouble> stress_ptr)
134  : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), strainPtr(strain_ptr),
135  tempPtr(temp_ptr), mDPtr(m_D_ptr), coeffExpansionPtr(coeff_expansion_ptr),
136  stressPtr(stress_ptr) {}
137 
138 //! [Calculate stress]
139 MoFEMErrorCode OpStressThermal::doWork(int side, EntityType type,
140  EntData &data) {
141  MoFEMFunctionBegin;
142  const auto nb_gauss_pts = getGaussPts().size2();
143  stressPtr->resize((SPACE_DIM * (SPACE_DIM + 1)) / 2, nb_gauss_pts);
144 
145  auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mDPtr);
146  auto t_strain = getFTensor2SymmetricFromMat<SPACE_DIM>(*strainPtr);
147  auto t_stress = getFTensor2SymmetricFromMat<SPACE_DIM>(*stressPtr);
148  auto t_temp = getFTensor0FromVec(*tempPtr);
149 
150  FTensor::Index<'i', SPACE_DIM> i;
151  FTensor::Index<'j', SPACE_DIM> j;
152  FTensor::Index<'k', SPACE_DIM> k;
153  FTensor::Index<'l', SPACE_DIM> l;
154 
155  auto t_coeff_exp = FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM>();
156  t_coeff_exp(i, j) = 0;
157  for (auto d = 0; d != SPACE_DIM; ++d) {
158  t_coeff_exp(d, d) = (*coeffExpansionPtr)[d];
159  }
160 
161  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
162  t_stress(i, j) = t_D(i, j, k, l) *
163  (t_strain(k, l) - t_coeff_exp(k, l) * (t_temp - ref_temp));
164  ++t_strain;
165  ++t_stress;
166  ++t_temp;
167  }
168  MoFEMFunctionReturn(0);
169 }
170 //! [Calculate stress]
171 
172 struct SetTargetTemperature;
173 
174 } // namespace ThermoElasticOps
175 
176 //! [Target temperature]
177 
178 template <AssemblyType A, IntegrationType I, typename OpBase>
179 struct MoFEM::OpFluxRhsImpl<ThermoElasticOps::SetTargetTemperature, 1, 1, A, I,
180  OpBase>;
181 
182 template <AssemblyType A, IntegrationType I, typename OpBase>
183 struct MoFEM::OpFluxLhsImpl<ThermoElasticOps::SetTargetTemperature, 1, 1, A, I,
184  OpBase>;
185 
186 template <AssemblyType A, IntegrationType I, typename OpBase>
187 struct MoFEM::AddFluxToRhsPipelineImpl<
188 
189  MoFEM::OpFluxRhsImpl<ThermoElasticOps::SetTargetTemperature, 1, 1, A, I,
190  OpBase>,
191  A, I, OpBase
192 
193  > {
194 
195  AddFluxToRhsPipelineImpl() = delete;
196 
197  static MoFEMErrorCode add(
198 
199  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
200  MoFEM::Interface &m_field, const std::string field_name,
201  boost::shared_ptr<VectorDouble> temp_ptr, std::string block_name, Sev sev
202 
203  ) {
204  MoFEMFunctionBegin;
205 
206  using OP_SOURCE = typename FormsIntegrators<OpBase>::template Assembly<
207  A>::template LinearForm<I>::template OpSource<1, 1>;
208  using OP_TEMP = typename FormsIntegrators<OpBase>::template Assembly<
209  A>::template LinearForm<I>::template OpBaseTimesScalar<1>;
210 
211  auto add_op = [&](auto &&meshset_vec_ptr) {
212  MoFEMFunctionBegin;
213  for (auto m : meshset_vec_ptr) {
214  std::vector<double> block_data;
215  m->getAttributes(block_data);
216  if (block_data.size() < 2) {
217  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
218  "Expected two parameters");
219  }
220  double target_temperature = block_data[0];
221  double beta =
222  block_data[1]; // Set temperature parameter [ W/K * (1/m^3)]
223  auto ents_ptr = boost::make_shared<Range>();
224  CHKERR m_field.get_moab().get_entities_by_handle(m->meshset,
225  *(ents_ptr), true);
226 
227  MOFEM_TAG_AND_LOG("WORLD", sev, "SetTargetTemperature")
228  << "Add " << *m << " target temperature " << target_temperature
229  << " penalty " << beta;
230 
231  pipeline.push_back(new OP_SOURCE(
232  field_name,
233  [target_temperature, beta](double, double, double) {
234  return target_temperature * beta;
235  },
236  ents_ptr));
237  pipeline.push_back(new OP_TEMP(
238  field_name, temp_ptr,
239  [beta](double, double, double) { return -beta; }, ents_ptr));
240  }
241  MoFEMFunctionReturn(0);
242  };
243 
244  CHKERR add_op(
245 
246  m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
247 
248  (boost::format("%s(.*)") % block_name).str()
249 
250  ))
251 
252  );
253 
254  MOFEM_LOG_CHANNEL("WORLD");
255 
256  MoFEMFunctionReturn(0);
257  }
258 };
259 
260 template <AssemblyType A, IntegrationType I, typename OpBase>
261 struct AddFluxToLhsPipelineImpl<
262 
263  OpFluxLhsImpl<ThermoElasticOps::SetTargetTemperature, 1, 1, A, I, OpBase>,
264  A, I, OpBase
265 
266  > {
267 
268  AddFluxToLhsPipelineImpl() = delete;
269 
270  static MoFEMErrorCode add(
271 
272  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
273  MoFEM::Interface &m_field, const std::string field_name,
274  boost::shared_ptr<VectorDouble> temp_ptr, std::string block_name, Sev sev
275 
276  ) {
277  MoFEMFunctionBegin;
278 
279  using OP_MASS = typename FormsIntegrators<OpBase>::template Assembly<
280  A>::template BiLinearForm<I>::template OpMass<1, 1>;
281 
282  auto add_op = [&](auto &&meshset_vec_ptr) {
283  MoFEMFunctionBegin;
284  for (auto m : meshset_vec_ptr) {
285  std::vector<double> block_data;
286  m->getAttributes(block_data);
287  if (block_data.size() != 2) {
288  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
289  "Expected two parameters");
290  }
291  double beta =
292  block_data[1]; // Set temperature parameter [ W/K * (1/m^3)]
293  auto ents_ptr = boost::make_shared<Range>();
294  CHKERR m_field.get_moab().get_entities_by_handle(m->meshset,
295  *(ents_ptr), true);
296 
297  MOFEM_TAG_AND_LOG("WORLD", sev, "SetTargetTemperature")
298  << "Add " << *m << " penalty " << beta;
299 
300  pipeline.push_back(new OP_MASS(
301  field_name, field_name,
302  [beta](double, double, double) { return -beta; }, ents_ptr));
303  }
304  MoFEMFunctionReturn(0);
305  };
306 
307  CHKERR add_op(
308 
309  m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
310 
311  (boost::format("%s(.*)") % block_name).str()
312 
313  ))
314 
315  );
316 
317  MOFEM_LOG_CHANNEL("WORLD");
318 
319  MoFEMFunctionReturn(0);
320  }
321 };
322 
323 //! [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:46
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::OpKCauchyThermoElasticity
Definition: ThermoElasticOps.hpp:9
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
ref_temp
double ref_temp
Definition: thermo_elastic.cpp:125
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< VectorDouble > coeff_expansion_ptr, boost::shared_ptr< MatrixDouble > stress_ptr)
Definition: ThermoElasticOps.hpp:115
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
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:44
MoFEM::Sev
MoFEM::LogManager::SeverityLevel Sev
Definition: CoreTemplates.hpp:17
ThermoElasticOps::OpStressThermal
Definition: ThermoElasticOps.hpp:23
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
AddFluxToLhsPipelineImpl< OpFluxLhsImpl< ThermoElasticOps::SetTargetTemperature, 1, 1, A, I, OpBase >, A, I, OpBase >::add
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)
Definition: ThermoElasticOps.hpp:270
ThermoElasticOps::OpKCauchyThermoElasticity::mDPtr
boost::shared_ptr< MatrixDouble > mDPtr
Definition: ThermoElasticOps.hpp:19
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
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
MoFEM::AddFluxToRhsPipelineImpl< MoFEM::OpFluxRhsImpl< ThermoElasticOps::SetTargetTemperature, 1, 1, A, I, OpBase >, A, I, OpBase >::add
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)
Definition: ThermoElasticOps.hpp:197
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:139
ThermoElasticOps::OpStressThermal::coeffExpansionPtr
boost::shared_ptr< VectorDouble > coeffExpansionPtr
Definition: ThermoElasticOps.hpp:47
MoFEM::LogManager::SeverityLevel
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
MoFEM::AddFluxToRhsPipelineImpl
Definition: Natural.hpp:46
ThermoElasticOps
Definition: ThermalConvection.hpp:12
ThermoElasticOps::OpKCauchyThermoElasticity::coeffExpansionPtr
boost::shared_ptr< VectorDouble > coeffExpansionPtr
Definition: ThermoElasticOps.hpp:20
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
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::OpKCauchyThermoElasticity::OpKCauchyThermoElasticity
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)
Definition: ThermoElasticOps.hpp:51
ThermoElasticOps::OpStressThermal::stressPtr
boost::shared_ptr< MatrixDouble > stressPtr
Definition: ThermoElasticOps.hpp:48
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
ThermoElasticOps::OpKCauchyThermoElasticity::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ThermoElasticOps.hpp:62
MoFEM::FormsIntegrators
Integrator forms.
Definition: FormsIntegrators.hpp:306
ThermoElasticOps::OpStressThermal::tempPtr
boost::shared_ptr< VectorDouble > tempPtr
Definition: ThermoElasticOps.hpp:45
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
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:578
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
Generated by Doxygen 1.8.17 and hosted at University of Glasgow