v0.14.0
PlasticOpsLargeStrains.hpp
Go to the documentation of this file.
1 /** \file PlasticOpsLargeStrains.hpp
2  * \example PlasticOpsLargeStrains.hpp
3  */
4 
5 namespace PlasticOps {
6 
7 template <int DIM, typename AssemblyDomainEleOp>
10  : public AssemblyDomainEleOp {
12  const std::string row_field_name, const std::string col_field_name,
13  boost::shared_ptr<CommonData> common_data_ptr,
14  boost::shared_ptr<HenckyOps::CommonData> common_henky_data_ptr,
15  boost::shared_ptr<MatrixDouble> m_D_ptr);
16  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
17  EntitiesFieldData::EntData &col_data);
18 
19 private:
20  boost::shared_ptr<CommonData> commonDataPtr;
21  boost::shared_ptr<HenckyOps::CommonData> commonHenckyDataPtr;
22  boost::shared_ptr<MatrixDouble> mDPtr;
25 };
26 
27 template <int DIM, typename AssemblyDomainEleOp>
30  OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl(
31  const std::string row_field_name, const std::string col_field_name,
32  boost::shared_ptr<CommonData> common_data_ptr,
33  boost::shared_ptr<HenckyOps::CommonData> common_henky_data_ptr,
34  boost::shared_ptr<MatrixDouble> m_D_ptr)
35  : AssemblyDomainEleOp(row_field_name, col_field_name,
36  AssemblyDomainEleOp::OPROWCOL),
37  commonDataPtr(common_data_ptr),
38  commonHenckyDataPtr(common_henky_data_ptr), mDPtr(m_D_ptr) {
39  AssemblyDomainEleOp::sYmm = false;
40 }
41 
42 template <int DIM, typename AssemblyDomainEleOp>
44  DIM, GAUSS,
46  EntitiesFieldData::EntData &col_data) {
48 
55 
56  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
58 
59  auto &locMat = AssemblyDomainEleOp::locMat;
60 
61  const size_t nb_integration_pts = row_data.getN().size1();
62  const size_t nb_row_base_functions = row_data.getN().size2();
63 
64  if (AssemblyDomainEleOp::rowType == MBVERTEX &&
66  resDiff.resize(DIM * DIM * size_symm, nb_integration_pts, false);
67  auto t_res_diff = getFTensor3FromMat<DIM, DIM, size_symm>(resDiff);
68  auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mDPtr);
69  auto t_logC_dC =
70  getFTensor4DdgFromMat<DIM, DIM>(commonHenckyDataPtr->matLogCdC);
71  auto t_grad =
72  getFTensor2FromMat<DIM, DIM>(*(commonHenckyDataPtr->matGradPtr));
73  auto t_L = symm_L_tensor(FTensor::Number<DIM>());
74  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
75  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
77  t_F(i, j) = t_grad(i, j) + t_kd(i, j);
79  t_DLogC_dC(i, j, k, l) = t_D(m, n, k, l) * t_logC_dC(m, n, i, j);
81  t_FDLogC_dC(i, j, k, l) = t_F(i, m) * t_DLogC_dC(m, j, k, l);
83  t_res_diff(i, j, L) = t_FDLogC_dC(i, j, k, l) * t_L(k, l, L);
84  ++t_logC_dC;
85  ++t_grad;
86  ++t_res_diff;
87  }
88  }
89 
90  auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
91  auto t_row_diff_base = row_data.getFTensor1DiffN<DIM>();
92  auto t_res_diff = getFTensor3FromMat<DIM, DIM, size_symm>(resDiff);
93 
94  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
95  double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
96 
97  size_t rr = 0;
98  for (; rr != AssemblyDomainEleOp::nbRows / DIM; ++rr) {
99 
100  auto t_mat =
102 
104  t_tmp(i, L) = (t_res_diff(i, j, L) * (alpha * t_row_diff_base(j)));
105 
106  auto t_col_base = col_data.getFTensor0N(gg, 0);
107  for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols / size_symm; ++cc) {
108 
109  t_mat(i, L) -= (t_col_base * t_tmp(i, L));
110 
111  ++t_mat;
112  ++t_col_base;
113  }
114 
115  ++t_row_diff_base;
116  }
117 
118  for (; rr < nb_row_base_functions; ++rr)
119  ++t_row_diff_base;
120 
121  ++t_w;
122  ++t_res_diff;
123  }
124 
126 }
127 
128 template <int DIM, typename AssemblyDomainEleOp>
131  : public AssemblyDomainEleOp {
133  const std::string row_field_name, const std::string col_field_name,
134  boost::shared_ptr<CommonData> common_data_ptr,
135  boost::shared_ptr<HenckyOps::CommonData> comman_henky_data_ptr,
136  boost::shared_ptr<MatrixDouble> m_D_ptr);
137  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
138  EntitiesFieldData::EntData &col_data);
139 
140 private:
141  boost::shared_ptr<CommonData> commonDataPtr;
142  boost::shared_ptr<HenckyOps::CommonData> commonHenckyDataPtr;
143  boost::shared_ptr<MatrixDouble> mDPtr;
145 };
146 
147 template <int DIM, typename AssemblyDomainEleOp>
150  const std::string row_field_name, const std::string col_field_name,
151  boost::shared_ptr<CommonData> common_data_ptr,
152  boost::shared_ptr<HenckyOps::CommonData> comman_henky_data_ptr,
153  boost::shared_ptr<MatrixDouble> m_D_ptr)
154  : AssemblyDomainEleOp(row_field_name, col_field_name,
155  AssemblyDomainEleOp::OPROWCOL),
156  commonDataPtr(common_data_ptr),
157  commonHenckyDataPtr(comman_henky_data_ptr), mDPtr(m_D_ptr) {
158  AssemblyDomainEleOp::sYmm = false;
159 }
160 
161 template <int DIM, typename AssemblyDomainEleOp>
165  EntitiesFieldData::EntData &col_data) {
167 
174  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
176 
177  auto &locMat = AssemblyDomainEleOp::locMat;
178 
179  const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
180  const auto nb_row_base_functions = row_data.getN().size2();
181 
182  if (AssemblyDomainEleOp::colType == MBVERTEX &&
184 
185  resDiff.resize(size_symm * DIM * DIM, nb_integration_pts, false);
186  auto t_res_diff = getFTensor3FromMat<size_symm, DIM, DIM>(resDiff);
187 
188  auto t_res_flow_dstrain =
189  getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrain);
190  auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->mGradPtr));
191  auto t_logC_dC =
192  getFTensor4DdgFromMat<DIM, DIM>(commonHenckyDataPtr->matLogCdC);
193 
194  auto next = [&]() {
195  ++t_res_flow_dstrain;
196  ++t_grad;
197  ++t_logC_dC;
198  ++t_res_diff;
199  };
200 
202  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
204  t_diff_grad(i, j, k, l) = t_kd(i, k) * t_kd(j, l);
205 
206  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
207  FTensor::Ddg<double, DIM, DIM> t_diff_ls_dlogC_dC;
210 
211  t_diff_ls_dlogC_dC(i, j, k, l) =
212  (t_res_flow_dstrain(i, j, m, n)) * (t_logC_dC(m, n, k, l) / 2);
213 
214  t_F(i, j) = t_grad(i, j) + t_kd(i, j);
215  t_dC_dF(i, j, k, l) = (t_kd(i, l) * t_F(k, j)) + (t_kd(j, l) * t_F(k, i));
216 
218  t_res_diff(L, i, j) =
219  (t_L(m, n, L) * t_diff_ls_dlogC_dC(m, n, k, l)) * t_dC_dF(k, l, i, j);
220  next();
221  }
222  }
223 
224  auto t_res_diff = getFTensor3FromMat<size_symm, DIM, DIM>(resDiff);
225 
226  auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
227  auto t_row_base = row_data.getFTensor0N();
228  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
229  double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
230  ++t_w;
231 
232  size_t rr = 0;
233  for (; rr != AssemblyDomainEleOp::nbRows / size_symm; ++rr) {
234  const auto row_base = alpha * t_row_base;
235  auto t_mat =
237  auto t_col_diff_base = col_data.getFTensor1DiffN<DIM>(gg, 0);
238  for (auto cc = 0; cc != AssemblyDomainEleOp::nbCols / DIM; ++cc) {
239  t_mat(L, l) += row_base * (t_res_diff(L, l, k) * t_col_diff_base(k));
240  ++t_mat;
241  ++t_col_diff_base;
242  }
243  ++t_row_base;
244  }
245 
246  for (; rr < nb_row_base_functions; ++rr)
247  ++t_row_base;
248 
249  ++t_res_diff;
250  }
251 
253 }
254 
255 template <int DIM, typename AssemblyDomainEleOp>
257  public AssemblyDomainEleOp {
259  const std::string row_field_name, const std::string col_field_name,
260  boost::shared_ptr<CommonData> common_data_ptr,
261  boost::shared_ptr<HenckyOps::CommonData> comman_henky_data_ptr,
262  boost::shared_ptr<MatrixDouble> m_D_ptr);
263  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
264  EntitiesFieldData::EntData &col_data);
265 
266 private:
267  boost::shared_ptr<CommonData> commonDataPtr;
268  boost::shared_ptr<HenckyOps::CommonData> commonHenckyDataPtr;
269  boost::shared_ptr<MatrixDouble> mDPtr;
271 };
272 
273 template <int DIM, typename AssemblyDomainEleOp>
276  const std::string row_field_name, const std::string col_field_name,
277  boost::shared_ptr<CommonData> common_data_ptr,
278  boost::shared_ptr<HenckyOps::CommonData> comman_henky_data_ptr,
279  boost::shared_ptr<MatrixDouble> m_D_ptr)
280  : AssemblyDomainEleOp(row_field_name, col_field_name,
281  DomainEleOp::OPROWCOL),
282  commonDataPtr(common_data_ptr),
283  commonHenckyDataPtr(comman_henky_data_ptr), mDPtr(m_D_ptr) {
284  AssemblyDomainEleOp::sYmm = false;
285 }
286 
287 template <int DIM, typename AssemblyDomainEleOp>
291  EntitiesFieldData::EntData &col_data) {
293 
298 
299  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
300 
302  false);
304 
305  const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
306  const auto nb_row_base_functions = row_data.getN().size2();
307 
308  if (AssemblyDomainEleOp::colType == MBVERTEX &&
310 
311  resDiff.resize(DIM * DIM, nb_integration_pts, false);
312  auto t_res_diff = getFTensor2FromMat<DIM, DIM>(resDiff);
313  auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->mGradPtr));
314  auto t_logC_dC =
315  getFTensor4DdgFromMat<DIM, DIM>(commonHenckyDataPtr->matLogCdC);
316  auto t_c_dstrain =
317  getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resCdStrain);
318 
319  auto next = [&]() {
320  ++t_grad;
321  ++t_logC_dC;
322  ++t_c_dstrain;
323  ++t_res_diff;
324  };
325 
326  auto t_diff_grad_symmetrise = diff_symmetrize(FTensor::Number<DIM>());
327 
328  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
332 
333  t_diff_ls_dlog_c(k, l) =
334  (t_c_dstrain(i, j)) * (t_logC_dC(i, j, k, l) / 2);
335  t_F(i, j) = t_grad(i, j) + t_kd(i, j);
336  t_dC_dF(i, j, k, l) = (t_kd(i, l) * t_F(k, j)) + (t_kd(j, l) * t_F(k, i));
337 
339  t_res_diff(i, j) = (t_diff_ls_dlog_c(k, l) * t_dC_dF(k, l, i, j));
340  next();
341  }
342  }
343 
344  auto t_res_diff = getFTensor2FromMat<DIM, DIM>(resDiff);
345 
346  auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
347  auto t_row_base = row_data.getFTensor0N();
348  for (auto gg = 0; gg != nb_integration_pts; ++gg) {
349  double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
350  ++t_w;
351 
352  auto t_mat =
353  getFTensor1FromPtr<DIM>(AssemblyDomainEleOp::locMat.data().data());
354  size_t rr = 0;
355  for (; rr != AssemblyDomainEleOp::nbRows; ++rr) {
356  const auto row_base = alpha * t_row_base;
357  auto t_col_diff_base = col_data.getFTensor1DiffN<DIM>(gg, 0);
358  for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols / DIM; cc++) {
359  t_mat(i) += row_base * (t_res_diff(i, j) * t_col_diff_base(j));
360  ++t_mat;
361  ++t_col_diff_base;
362  }
363  ++t_row_base;
364  }
365  for (; rr != nb_row_base_functions; ++rr)
366  ++t_row_base;
367 
368  ++t_res_diff;
369  }
370 
372 }
373 
374 }; // namespace PlasticOps
PlasticOps::OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl< DIM, GAUSS, AssemblyDomainEleOp >::mDPtr
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOpsLargeStrains.hpp:22
MoFEM::OpBaseImpl::locMat
MatrixDouble locMat
local entity block matrix
Definition: FormsIntegrators.hpp:249
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
PlasticOps::OpCalculateConstraintsLhs_LogStrain_dUImpl< DIM, GAUSS, AssemblyDomainEleOp >::commonHenckyDataPtr
boost::shared_ptr< HenckyOps::CommonData > commonHenckyDataPtr
Definition: PlasticOpsLargeStrains.hpp:268
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
PlasticOps::OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl< DIM, GAUSS, AssemblyDomainEleOp >::resDiff
MatrixDouble resDiff
Definition: PlasticOpsLargeStrains.hpp:24
PlasticOps::OpCalculateConstraintsLhs_LogStrain_dUImpl
Definition: PlasticOps.hpp:157
PlasticOps::OpCalculatePlasticFlowLhs_LogStrain_dUImpl< DIM, GAUSS, AssemblyDomainEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOpsLargeStrains.hpp:141
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
FTensor::Number
Definition: Number.hpp:11
FTensor::Tensor3
Definition: Tensor3_value.hpp:12
PlasticOps::OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl
Definition: PlasticOps.hpp:139
PlasticOps::OpCalculatePlasticFlowLhs_LogStrain_dUImpl< DIM, GAUSS, AssemblyDomainEleOp >::resDiff
MatrixDouble resDiff
Definition: PlasticOpsLargeStrains.hpp:144
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
MoFEM::OpBaseImpl::colSide
int colSide
column side number
Definition: FormsIntegrators.hpp:242
PlasticOps::OpCalculatePlasticFlowLhs_LogStrain_dUImpl< DIM, GAUSS, AssemblyDomainEleOp >::mDPtr
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOpsLargeStrains.hpp:143
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
PlasticOps::get_mat_tensor_sym_dvector
static auto get_mat_tensor_sym_dvector(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
[Lambda functions]
Definition: PlasticOpsGeneric.hpp:366
FTensor::Tensor4
Definition: Tensor4_value.hpp:18
size_symm
constexpr auto size_symm
Definition: plastic.cpp:42
PlasticOps::OpCalculatePlasticFlowLhs_LogStrain_dUImpl< DIM, GAUSS, AssemblyDomainEleOp >::commonHenckyDataPtr
boost::shared_ptr< HenckyOps::CommonData > commonHenckyDataPtr
Definition: PlasticOpsLargeStrains.hpp:142
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
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
FTensor::Index< 'i', DIM >
convert.n
n
Definition: convert.py:82
MoFEM::OpBaseImpl::rowSide
int rowSide
row side number
Definition: FormsIntegrators.hpp:241
MoFEM::OpBaseImpl::nbCols
int nbCols
number if dof on column
Definition: FormsIntegrators.hpp:237
PlasticOps::OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl< DIM, GAUSS, AssemblyDomainEleOp >::commonHenckyDataPtr
boost::shared_ptr< HenckyOps::CommonData > commonHenckyDataPtr
Definition: PlasticOpsLargeStrains.hpp:21
DomainEleOp
PlasticOps::OpCalculatePlasticFlowLhs_LogStrain_dUImpl
Definition: PlasticOps.hpp:145
PlasticOps
Definition: PlasticNaturalBCs.hpp:13
MoFEM::OpBaseImpl::nbRows
int nbRows
number of dofs on rows
Definition: FormsIntegrators.hpp:236
PlasticOps::symm_L_tensor
auto symm_L_tensor(FTensor::Number< DIM >)
Definition: PlasticOpsGeneric.hpp:21
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
FTensor::Ddg
Definition: Ddg_value.hpp:7
PlasticOps::OpCalculateConstraintsLhs_LogStrain_dUImpl< DIM, GAUSS, AssemblyDomainEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOpsLargeStrains.hpp:267
PlasticOps::OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl< DIM, GAUSS, AssemblyDomainEleOp >::locMat
MatrixDouble locMat
Definition: PlasticOpsLargeStrains.hpp:23
PlasticOps::OpCalculateConstraintsLhs_LogStrain_dUImpl< DIM, GAUSS, AssemblyDomainEleOp >::resDiff
MatrixDouble resDiff
Definition: PlasticOpsLargeStrains.hpp:270
PlasticOps::diff_symmetrize
auto diff_symmetrize(FTensor::Number< DIM >)
Definition: PlasticOpsGeneric.hpp:43
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
PlasticOps::get_mat_vector_dtensor_sym
static FTensor::Tensor2< FTensor::PackPtr< double *, 3 >, 2, 3 > get_mat_vector_dtensor_sym(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
Definition: PlasticOpsSmallStrains.hpp:38
MoFEM::OpBaseImpl::rowType
EntityType rowType
row type
Definition: FormsIntegrators.hpp:243
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
PlasticOps::OpCalculateConstraintsLhs_LogStrain_dUImpl< DIM, GAUSS, AssemblyDomainEleOp >::mDPtr
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOpsLargeStrains.hpp:269
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
PlasticOps::OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl< DIM, GAUSS, AssemblyDomainEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOpsLargeStrains.hpp:20
MoFEM::OpBaseImpl::colType
EntityType colType
column type
Definition: FormsIntegrators.hpp:244
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21