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