v0.14.0
PlasticOps.hpp
Go to the documentation of this file.
1 
2 
3 /** \file PlasticOps.hpp
4  * \example PlasticOps.hpp
5 
6 \f[
7 \left\{
8 \begin{array}{ll}
9 \frac{\partial \sigma_{ij}}{\partial x_j} - b_i = 0 & \forall x \in \Omega \\
10 \varepsilon_{ij} = \frac{1}{2}\left( \frac{\partial u_i}{\partial x_j} +
11 \frac{\partial u_j}{\partial x_i} \right)\\
12 \sigma_{ij} = D_{ijkl}\left(\varepsilon_{kl}-\varepsilon^p_{kl}\right) \\
13 \dot{\varepsilon}^p_{kl} - \dot{\tau} \left( \left. \frac{\partial f}{\partial
14 \sigma_{kl}} \right|_{(\sigma,\tau) } \right) = 0 \\
15 f(\sigma, \tau) \leq 0,\; \dot{\tau} \geq 0,\;\dot{\tau}f(\sigma, \tau)=0\\
16 u_i = \overline{u}_i & \forall x \in \partial\Omega_u \\
17 \sigma_{ij}n_j = \overline{t}_i & \forall x \in \partial\Omega_\sigma \\
18 \Omega_u \cup \Omega_\sigma = \Omega \\
19 \Omega_u \cap \Omega_\sigma = \emptyset
20 \end{array}
21 \right.
22 \f]
23 
24 \f[
25 \left\{
26 \begin{array}{ll}
27 \left(\frac{\partial \delta u_i}{\partial x_j},\sigma_{ij}\right)_\Omega-(\delta
28 u_i,b_i)_\Omega -(\delta u_i,\overline{t}_i)_{\partial\Omega_\sigma}=0 & \forall
29 \delta u_i \in H^1(\Omega)\\ \left(\delta\varepsilon^p_{kl} ,D_{ijkl}\left(
30 \dot{\varepsilon}^p_{kl} - \dot{\tau} A_{kl} \right)\right) = 0
31 & \forall \delta\varepsilon^p_{ij} \in L^2(\Omega) \cap \mathcal{S} \\
32 \left(\delta\tau,c_n\dot{\tau} - \frac{1}{2}\left\{c_n \dot{\tau} +
33 (f(\pmb\sigma,\tau) - \sigma_y) +
34 \| c_n \dot{\tau} + (f(\pmb\sigma,\tau) - \sigma_y) \|\right\}\right) = 0 &
35 \forall \delta\tau \in L^2(\Omega) \end{array} \right.
36 \f]
37 
38 */
39 
40 namespace PlasticOps {
41 
42 //! [Common data]
43 struct CommonData : public boost::enable_shared_from_this<CommonData> {
44 
45  enum ParamsIndexes {
46  YOUNG_MODULUS,
47  POISSON_RATIO,
48  SIGMA_Y,
49  H,
50  VIS_H,
51  QINF,
52  BISO,
53  C1_k,
54  LAST_PARAM
55  };
56 
57  using BlockParams = std::array<double, LAST_PARAM>;
58  BlockParams blockParams;
59 
60  inline auto getParamsPtr() {
61  return boost::shared_ptr<BlockParams>(shared_from_this(), &blockParams);
62  };
63 
64  //! [Common data set externally]
65  boost::shared_ptr<MatrixDouble> mDPtr;
66  boost::shared_ptr<MatrixDouble> mGradPtr;
67  boost::shared_ptr<MatrixDouble> mStrainPtr;
68  boost::shared_ptr<MatrixDouble> mStressPtr;
69  //! [Common data set externally]
70 
71  VectorDouble plasticSurface;
72  MatrixDouble plasticFlow;
73  VectorDouble plasticTau;
74  VectorDouble plasticTauDot;
75  MatrixDouble plasticStrain;
76  MatrixDouble plasticStrainDot;
77 
78  VectorDouble resC;
79  VectorDouble resCdTau;
80  MatrixDouble resCdStrain;
81  MatrixDouble resCdPlasticStrain;
82  MatrixDouble resFlow;
83  MatrixDouble resFlowDtau;
84  MatrixDouble resFlowDstrain;
85  MatrixDouble resFlowDstrainDot;
86 
87  inline auto getPlasticSurfacePtr() {
88  return boost::shared_ptr<VectorDouble>(shared_from_this(), &plasticSurface);
89  }
90  inline auto getPlasticTauPtr() {
91  return boost::shared_ptr<VectorDouble>(shared_from_this(), &plasticTau);
92  }
93  inline auto getPlasticTauDotPtr() {
94  return boost::shared_ptr<VectorDouble>(shared_from_this(), &plasticTauDot);
95  }
96  inline auto getPlasticStrainPtr() {
97  return boost::shared_ptr<MatrixDouble>(shared_from_this(), &plasticStrain);
98  }
99  inline auto getPlasticStrainDotPtr() {
100  return boost::shared_ptr<MatrixDouble>(shared_from_this(),
101  &plasticStrainDot);
102  }
103  inline auto getPlasticFlowPtr() {
104  return boost::shared_ptr<MatrixDouble>(shared_from_this(), &plasticFlow);
105  }
106 
107  static std::array<int, 5> activityData;
108 };
109 
110 std::array<int, 5> CommonData::activityData = {0, 0, 0, 0, 0};
111 
112 //! [Common data]
113 
114 
115 FTensor::Index<'I', 3> I;
116 FTensor::Index<'J', 3> J;
117 FTensor::Index<'M', 3> M;
118 FTensor::Index<'N', 3> N;
119 
120 template <int DIM, IntegrationType I, typename DomainEleOp>
121 struct OpCalculatePlasticSurfaceImpl;
122 
123 template <int DIM, IntegrationType I, typename DomainEleOp>
124 struct OpCalculatePlasticityImpl;
125 
126 template <int DIM, IntegrationType I, typename DomainEleOp>
127 struct OpPlasticStressImpl;
128 
129 template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
130 struct OpCalculatePlasticFlowRhsImpl;
131 
132 template <IntegrationType I, typename AssemblyDomainEleOp>
133 struct OpCalculateConstraintsRhsImpl;
134 
135 template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
136 struct OpCalculatePlasticInternalForceLhs_dEPImpl;
137 
138 template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
139 struct OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl;
140 
141 template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
142 struct OpCalculatePlasticFlowLhs_dUImpl;
143 
144 template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
145 struct OpCalculatePlasticFlowLhs_LogStrain_dUImpl;
146 
147 template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
148 struct OpCalculatePlasticFlowLhs_dEPImpl;
149 
150 template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
151 struct OpCalculatePlasticFlowLhs_dTAUImpl;
152 
153 template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
154 struct OpCalculateConstraintsLhs_dUImpl;
155 
156 template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
157 struct OpCalculateConstraintsLhs_LogStrain_dUImpl;
158 
159 template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
160 struct OpCalculateConstraintsLhs_dEPImpl;
161 
162 template <IntegrationType I, typename AssemblyDomainEleOp>
163 struct OpCalculateConstraintsLhs_dTAUImpl;
164 
165 template <typename DomainEleOp> struct PlasticityIntegrators {
166  template <int DIM, IntegrationType I>
167  using OpCalculatePlasticSurface =
168  OpCalculatePlasticSurfaceImpl<DIM, I, DomainEleOp>;
169 
170  template <int DIM, IntegrationType I>
171  using OpCalculatePlasticity = OpCalculatePlasticityImpl<DIM, I, DomainEleOp>;
172 
173  template <int DIM, IntegrationType I>
174  using OpPlasticStress = OpPlasticStressImpl<DIM, I, DomainEleOp>;
175 
176  template <AssemblyType A> struct Assembly {
177 
178  using AssemblyDomainEleOp =
179  typename FormsIntegrators<DomainEleOp>::template Assembly<A>::OpBase;
180 
181  template <int DIM, IntegrationType I>
182  using OpCalculatePlasticFlowRhs =
183  OpCalculatePlasticFlowRhsImpl<DIM, I, AssemblyDomainEleOp>;
184 
185  template <IntegrationType I>
186  using OpCalculateConstraintsRhs =
187  OpCalculateConstraintsRhsImpl<I, AssemblyDomainEleOp>;
188 
189  template <int DIM, IntegrationType I>
190  using OpCalculatePlasticInternalForceLhs_dEP =
191  OpCalculatePlasticInternalForceLhs_dEPImpl<DIM, I, AssemblyDomainEleOp>;
192 
193  template <int DIM, IntegrationType I>
194  using OpCalculatePlasticInternalForceLhs_LogStrain_dEP =
195  OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl<
196  DIM, I, AssemblyDomainEleOp>;
197 
198  template <int DIM, IntegrationType I>
199  using OpCalculatePlasticFlowLhs_dU =
200  OpCalculatePlasticFlowLhs_dUImpl<DIM, I, AssemblyDomainEleOp>;
201 
202  template <int DIM, IntegrationType I>
203  using OpCalculatePlasticFlowLhs_LogStrain_dU =
204  OpCalculatePlasticFlowLhs_LogStrain_dUImpl<DIM, I, AssemblyDomainEleOp>;
205 
206  template <int DIM, IntegrationType I>
207  using OpCalculatePlasticFlowLhs_dEP =
208  OpCalculatePlasticFlowLhs_dEPImpl<DIM, I, AssemblyDomainEleOp>;
209 
210  template <int DIM, IntegrationType I>
211  using OpCalculatePlasticFlowLhs_dTAU =
212  OpCalculatePlasticFlowLhs_dTAUImpl<DIM, I, AssemblyDomainEleOp>;
213 
214  template <int DIM, IntegrationType I>
215  using OpCalculateConstraintsLhs_dU =
216  OpCalculateConstraintsLhs_dUImpl<DIM, I, AssemblyDomainEleOp>;
217 
218  template <int DIM, IntegrationType I>
219  using OpCalculateConstraintsLhs_LogStrain_dU =
220  OpCalculateConstraintsLhs_LogStrain_dUImpl<DIM, I, AssemblyDomainEleOp>;
221 
222  template <int DIM, IntegrationType I>
223  using OpCalculateConstraintsLhs_dEP =
224  OpCalculateConstraintsLhs_dEPImpl<DIM, I, AssemblyDomainEleOp>;
225 
226  template <IntegrationType I>
227  using OpCalculateConstraintsLhs_dTAU =
228  OpCalculateConstraintsLhs_dTAUImpl<I, AssemblyDomainEleOp>;
229 
230  };
231 };
232 
233 }; // namespace PlasticOps
234 
235 #include <PlasticOpsGeneric.hpp>
236 #include <PlasticOpsSmallStrains.hpp>
237 #include <PlasticOpsLargeStrains.hpp>
238 #include <PlasticOpsMonitor.hpp>
239 
240 namespace PlasticOps {
241 
242 using Pip = boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator>;
243 using CommonPlasticPtr = boost::shared_ptr<PlasticOps::CommonData>;
244 using CommonHenkyPtr = boost::shared_ptr<HenckyOps::CommonData>;
245 
246 template <int DIM>
247 MoFEMErrorCode
248 addMatBlockOps(MoFEM::Interface &m_field, std::string block_name, Pip &pip,
249  boost::shared_ptr<MatrixDouble> mat_D_Ptr,
250  boost::shared_ptr<CommonData::BlockParams> mat_params_ptr,
251  double scale_value, Sev sev) {
252  MoFEMFunctionBegin;
253 
254  struct OpMatBlocks : public DomainEleOp {
255  OpMatBlocks(boost::shared_ptr<MatrixDouble> m_D_ptr,
256  boost::shared_ptr<CommonData::BlockParams> mat_params_ptr,
257  double scale_value, MoFEM::Interface &m_field, Sev sev,
258  std::vector<const CubitMeshSets *> meshset_vec_ptr)
259  : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), matDPtr(m_D_ptr),
260  matParamsPtr(mat_params_ptr), scaleVal(scale_value) {
261  CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
262  "Can not get data from block");
263  }
264 
265  MoFEMErrorCode doWork(int side, EntityType type,
266  EntitiesFieldData::EntData &data) {
267  MoFEMFunctionBegin;
268 
269  auto getK = [](auto &p) {
270  auto young_modulus = p[CommonData::YOUNG_MODULUS];
271  auto poisson_ratio = p[CommonData::POISSON_RATIO];
272  return young_modulus / (3 * (1 - 2 * poisson_ratio));
273  };
274 
275  auto getG = [](auto &p) {
276  auto young_modulus = p[CommonData::YOUNG_MODULUS];
277  auto poisson_ratio = p[CommonData::POISSON_RATIO];
278  return young_modulus / (2 * (1 + poisson_ratio));
279  };
280 
281  auto scale_fun = [this](auto &p) {
282  p[CommonData::YOUNG_MODULUS] *= scaleVal;
283  p[CommonData::SIGMA_Y] *= scaleVal;
284  p[CommonData::H] *= scaleVal;
285  p[CommonData::VIS_H] *= scaleVal;
286  p[CommonData::QINF] *= scaleVal;
287  p[CommonData::C1_k] *= scaleVal;
288  };
289 
290  for (auto &b : blockData) {
291  if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
292  *matParamsPtr = b.bParams;
293  scale_fun(*matParamsPtr);
294  CHKERR getMatDPtr(matDPtr, getK(*matParamsPtr), getG(*matParamsPtr));
295  MoFEMFunctionReturnHot(0);
296  }
297  }
298 
299  (*matParamsPtr) = {young_modulus, poisson_ratio, sigmaY, H,
300  visH, Qinf, b_iso, C1_k};
301  scale_fun(*matParamsPtr);
302  CHKERR getMatDPtr(matDPtr, getK(*matParamsPtr), getG(*matParamsPtr));
303 
304  MoFEMFunctionReturn(0);
305  }
306 
307  private:
308  boost::shared_ptr<MatrixDouble> matDPtr;
309  boost::shared_ptr<CommonData::BlockParams> matParamsPtr;
310  const double scaleVal;
311 
312  struct BlockData {
313  std::array<double, CommonData::LAST_PARAM> bParams;
314  Range blockEnts;
315  };
316  std::vector<BlockData> blockData;
317 
318  /**
319  * @brief Extract block data from meshsets
320  *
321  * @param m_field
322  * @param meshset_vec_ptr
323  * @param sev
324  * @return MoFEMErrorCode
325  */
326  MoFEMErrorCode
327  extractBlockData(MoFEM::Interface &m_field,
328  std::vector<const CubitMeshSets *> meshset_vec_ptr,
329  Sev sev) {
330  MoFEMFunctionBegin;
331 
332  for (auto m : meshset_vec_ptr) {
333  MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
334  std::vector<double> block_data;
335  CHKERR m->getAttributes(block_data);
336  if (block_data.size() != CommonData::LAST_PARAM) {
337  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
338  "Wrong number of block attribute");
339  }
340  auto get_block_ents = [&]() {
341  Range ents;
342  CHKERR m_field.get_moab().get_entities_by_handle(m->meshset, ents,
343  true);
344  return ents;
345  };
346 
347  CommonData::BlockParams block_params;
348  for (auto i = 0; i != CommonData::LAST_PARAM; ++i) {
349  block_params[i] = block_data[i];
350  }
351 
352  MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
353  << "E = " << block_params[CommonData::YOUNG_MODULUS]
354  << " nu = " << block_params[CommonData::POISSON_RATIO];
355  MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
356  << std::endl
357  << "sigma_y = " << block_params[CommonData::SIGMA_Y] << std::endl
358  << "h = " << block_params[CommonData::H] << std::endl
359  << "vis_h = " << block_params[CommonData::VIS_H] << std::endl
360  << "qinf = " << block_params[CommonData::QINF] << std::endl
361  << "biso = " << block_params[CommonData::BISO] << std::endl
362  << "C1_k = " << block_params[CommonData::C1_k] << std::endl;
363 
364  blockData.push_back({block_params, get_block_ents()});
365  }
366  MOFEM_LOG_CHANNEL("WORLD");
367  MoFEMFunctionReturn(0);
368  }
369 
370  /**
371  * @brief Get elasticity tensor
372  *
373  * Calculate elasticity tensor for given material parameters
374  *
375  * @param mat_D_ptr
376  * @param bulk_modulus_K
377  * @param shear_modulus_G
378  * @return MoFEMErrorCode
379  *
380  */
381  MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
382  double bulk_modulus_K, double shear_modulus_G) {
383  MoFEMFunctionBegin;
384  //! [Calculate elasticity tensor]
385  auto set_material_stiffness = [&]() {
386  FTensor::Index<'i', DIM> i;
387  FTensor::Index<'j', DIM> j;
388  FTensor::Index<'k', DIM> k;
389  FTensor::Index<'l', DIM> l;
390  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
391  double A = (DIM == 2)
392  ? 2 * shear_modulus_G /
393  (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
394  : 1;
395  auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mat_D_ptr);
396  t_D(i, j, k, l) =
397  2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
398  A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
399  t_kd(k, l);
400  };
401  //! [Calculate elasticity tensor]
402  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
403  mat_D_ptr->resize(size_symm * size_symm, 1);
404  set_material_stiffness();
405  MoFEMFunctionReturn(0);
406  }
407  };
408 
409  // push operator to calculate material stiffness matrix for each block
410  pip.push_back(new OpMatBlocks(
411  mat_D_Ptr, mat_params_ptr, scale_value, m_field, sev,
412 
413  // Get blockset using regular expression
414  m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
415 
416  (boost::format("%s(.*)") % block_name).str()
417 
418  ))
419 
420  ));
421 
422  MoFEMFunctionReturn(0);
423 }
424 
425 template <int DIM, IntegrationType I, typename DomainEleOp>
426 auto createCommonPlasticOps(
427  MoFEM::Interface &m_field, std::string block_name,
428  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
429  std::string u, std::string ep, std::string tau, double scale, Sev sev) {
430 
431  using P = PlasticityIntegrators<DomainEleOp>;
432 
433  auto common_plastic_ptr = boost::make_shared<PlasticOps::CommonData>();
434  common_plastic_ptr = boost::make_shared<PlasticOps::CommonData>();
435 
436  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
437  auto make_d_mat = []() {
438  return boost::make_shared<MatrixDouble>(size_symm * size_symm, 1);
439  };
440 
441  common_plastic_ptr->mDPtr = make_d_mat();
442  common_plastic_ptr->mGradPtr = boost::make_shared<MatrixDouble>();
443  common_plastic_ptr->mStrainPtr = boost::make_shared<MatrixDouble>();
444  common_plastic_ptr->mStressPtr = boost::make_shared<MatrixDouble>();
445 
446  auto m_D_ptr = common_plastic_ptr->mDPtr;
447 
448  CHK_THROW_MESSAGE(addMatBlockOps<DIM>(m_field, block_name, pip, m_D_ptr,
449  common_plastic_ptr->getParamsPtr(),
450  scale, sev),
451  "add mat block ops");
452 
453  pip.push_back(new OpCalculateScalarFieldValues(
454  tau, common_plastic_ptr->getPlasticTauPtr()));
455  pip.push_back(new OpCalculateTensor2SymmetricFieldValues<DIM>(
456  ep, common_plastic_ptr->getPlasticStrainPtr()));
457  pip.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(
458  u, common_plastic_ptr->mGradPtr));
459 
460  CommonHenkyPtr common_henky_ptr;
461 
462  if (is_large_strains) {
463  common_henky_ptr = boost::make_shared<HenckyOps::CommonData>();
464  common_henky_ptr->matGradPtr = common_plastic_ptr->mGradPtr;
465  common_henky_ptr->matDPtr = common_plastic_ptr->mDPtr;
466  common_henky_ptr->matLogCPlastic =
467  common_plastic_ptr->getPlasticStrainPtr();
468  common_plastic_ptr->mStrainPtr = common_henky_ptr->getMatLogC();
469  common_plastic_ptr->mStressPtr = common_henky_ptr->getMatHenckyStress();
470 
471  using H = HenckyOps::HenkyIntegrators<DomainEleOp>;
472 
473  pip.push_back(new typename H::template OpCalculateEigenVals<DIM, I>(
474  u, common_henky_ptr));
475  pip.push_back(
476  new typename H::template OpCalculateLogC<DIM, I>(u, common_henky_ptr));
477  pip.push_back(new typename H::template OpCalculateLogC_dC<DIM, I>(
478  u, common_henky_ptr));
479  pip.push_back(new
480  typename H::template OpCalculateHenckyPlasticStress<DIM, I>(
481  u, common_henky_ptr, m_D_ptr));
482  pip.push_back(new typename H::template OpCalculatePiolaStress<DIM, I>(
483  u, common_henky_ptr));
484  } else {
485 
486  pip.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
487  common_plastic_ptr->mGradPtr, common_plastic_ptr->mStrainPtr));
488  pip.push_back(new typename P::template OpPlasticStress<DIM, I>(
489  common_plastic_ptr, m_D_ptr));
490  }
491 
492  pip.push_back(new typename P::template OpCalculatePlasticSurface<DIM, I>(
493  u, common_plastic_ptr));
494 
495  return std::make_tuple(common_plastic_ptr, common_henky_ptr);
496 }
497 
498 template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
499 MoFEMErrorCode
500 opFactoryDomainRhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip,
501  std::string u, std::string ep, std::string tau) {
502  MoFEMFunctionBegin;
503 
504  using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
505  A>::template LinearForm<I>;
506  using OpInternalForceCauchy =
507  typename B::template OpGradTimesSymTensor<1, DIM, DIM>;
508  using OpInternalForcePiola =
509  typename B::template OpGradTimesTensor<1, DIM, DIM>;
510 
511  using P = PlasticityIntegrators<DomainEleOp>;
512 
513  auto [common_plastic_ptr, common_henky_ptr] =
514  createCommonPlasticOps<DIM, I, DomainEleOp>(m_field, block_name, pip, u,
515  ep, tau, scale, Sev::inform);
516 
517  auto m_D_ptr = common_plastic_ptr->mDPtr;
518 
519  pip.push_back(new OpCalculateTensor2SymmetricFieldValuesDot<DIM>(
520  ep, common_plastic_ptr->getPlasticStrainDotPtr()));
521  pip.push_back(new OpCalculateScalarFieldValuesDot(
522  tau, common_plastic_ptr->getPlasticTauDotPtr()));
523  pip.push_back(new typename P::template OpCalculatePlasticity<DIM, I>(
524  u, common_plastic_ptr, m_D_ptr));
525 
526  // Calculate internal forces
527  if (common_henky_ptr) {
528  pip.push_back(new OpInternalForcePiola(
529  u, common_henky_ptr->getMatFirstPiolaStress()));
530  } else {
531  pip.push_back(new OpInternalForceCauchy(u, common_plastic_ptr->mStressPtr));
532  }
533 
534  pip.push_back(
535  new
536  typename P::template Assembly<A>::template OpCalculateConstraintsRhs<I>(
537  tau, common_plastic_ptr, m_D_ptr));
538  pip.push_back(
539  new typename P::template Assembly<A>::template OpCalculatePlasticFlowRhs<
540  DIM, I>(ep, common_plastic_ptr, m_D_ptr));
541 
542  MoFEMFunctionReturn(0);
543 }
544 
545 template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
546 MoFEMErrorCode
547 opFactoryDomainLhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip,
548  std::string u, std::string ep, std::string tau) {
549  MoFEMFunctionBegin;
550 
551  using namespace HenckyOps;
552 
553  using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
554  A>::template BiLinearForm<I>;
555  using OpKPiola = typename B::template OpGradTensorGrad<1, DIM, DIM, 1>;
556  using OpKCauchy = typename B::template OpGradSymTensorGrad<1, DIM, DIM, 0>;
557 
558  using P = PlasticityIntegrators<DomainEleOp>;
559 
560  auto [common_plastic_ptr, common_henky_ptr] =
561  createCommonPlasticOps<DIM, I, DomainEleOp>(m_field, block_name, pip, u,
562  ep, tau, scale, Sev::verbose);
563 
564  auto m_D_ptr = common_plastic_ptr->mDPtr;
565 
566  pip.push_back(new OpCalculateTensor2SymmetricFieldValuesDot<DIM>(
567  ep, common_plastic_ptr->getPlasticStrainDotPtr()));
568  pip.push_back(new OpCalculateScalarFieldValuesDot(
569  tau, common_plastic_ptr->getPlasticTauDotPtr()));
570  pip.push_back(new typename P::template OpCalculatePlasticity<DIM, I>(
571  u, common_plastic_ptr, m_D_ptr));
572 
573  if (common_henky_ptr) {
574  using H = HenckyOps::HenkyIntegrators<DomainEleOp>;
575  pip.push_back(new typename H::template OpHenckyTangent<DIM, I>(
576  u, common_henky_ptr, m_D_ptr));
577  pip.push_back(new OpKPiola(u, u, common_henky_ptr->getMatTangent()));
578  pip.push_back(
579  new typename P::template Assembly<A>::
580  template OpCalculatePlasticInternalForceLhs_LogStrain_dEP<DIM, I>(
581  u, ep, common_plastic_ptr, common_henky_ptr, m_D_ptr));
582  } else {
583  pip.push_back(new OpKCauchy(u, u, m_D_ptr));
584  pip.push_back(new typename P::template Assembly<A>::
585  template OpCalculatePlasticInternalForceLhs_dEP<DIM, I>(
586  u, ep, common_plastic_ptr, m_D_ptr));
587  }
588 
589  if (common_henky_ptr) {
590  pip.push_back(
591  new typename P::template Assembly<A>::
592  template OpCalculateConstraintsLhs_LogStrain_dU<DIM, I>(
593  tau, u, common_plastic_ptr, common_henky_ptr, m_D_ptr));
594  pip.push_back(
595  new typename P::template Assembly<A>::
596  template OpCalculatePlasticFlowLhs_LogStrain_dU<DIM, I>(
597  ep, u, common_plastic_ptr, common_henky_ptr, m_D_ptr));
598  } else {
599  pip.push_back(
600  new
601  typename P::template Assembly<A>::template OpCalculateConstraintsLhs_dU<
602  DIM, I>(tau, u, common_plastic_ptr, m_D_ptr));
603  pip.push_back(
604  new
605  typename P::template Assembly<A>::template OpCalculatePlasticFlowLhs_dU<
606  DIM, I>(ep, u, common_plastic_ptr, m_D_ptr));
607  }
608 
609  pip.push_back(
610  new
611  typename P::template Assembly<A>::template OpCalculatePlasticFlowLhs_dEP<
612  DIM, I>(ep, ep, common_plastic_ptr, m_D_ptr));
613  pip.push_back(
614  new
615  typename P::template Assembly<A>::template OpCalculatePlasticFlowLhs_dTAU<
616  DIM, I>(ep, tau, common_plastic_ptr, m_D_ptr));
617  pip.push_back(
618  new
619  typename P::template Assembly<A>::template OpCalculateConstraintsLhs_dEP<
620  DIM, I>(tau, ep, common_plastic_ptr, m_D_ptr));
621  pip.push_back(
622  new
623  typename P::template Assembly<A>::template OpCalculateConstraintsLhs_dTAU<
624  I>(tau, tau, common_plastic_ptr));
625 
626  MoFEMFunctionReturn(0);
627 }
628 
629 template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
630 MoFEMErrorCode opFactoryDomainReactions(MoFEM::Interface &m_field,
631  std::string block_name, Pip &pip,
632  std::string u, std::string ep,
633  std::string tau) {
634  MoFEMFunctionBegin;
635 
636  using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
637  A>::template LinearForm<I>;
638  using OpInternalForceCauchy =
639  typename B::template OpGradTimesSymTensor<1, DIM, DIM>;
640  using OpInternalForcePiola =
641  typename B::template OpGradTimesTensor<1, DIM, DIM>;
642 
643  auto [common_plastic_ptr, common_henky_ptr] =
644  createCommonPlasticOps<DIM, I, DomainEleOp>(m_field, block_name, pip, u,
645  ep, tau, 1., Sev::inform);
646 
647  // Calculate internal forces
648  if (common_henky_ptr) {
649  pip.push_back(new OpInternalForcePiola(
650  u, common_henky_ptr->getMatFirstPiolaStress()));
651  } else {
652  pip.push_back(new OpInternalForceCauchy(u, common_plastic_ptr->mStressPtr));
653  }
654 
655  MoFEMFunctionReturn(0);
656 }
657 
658 } // namespace PlasticOps
PlasticOps::OpCalculateConstraintsRhsImpl
Definition: PlasticOps.hpp:133
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
PlasticOps::CommonData::BISO
@ BISO
Definition: PlasticOps.hpp:52
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
PlasticOps::OpCalculateConstraintsLhs_dTAUImpl
Definition: PlasticOps.hpp:163
PlasticOps::CommonData::resCdStrain
MatrixDouble resCdStrain
Definition: PlasticOps.hpp:80
PlasticOps::J
FTensor::Index< 'J', 3 > J
Definition: PlasticOps.hpp:116
PlasticOps::CommonData::BlockParams
std::array< double, LAST_PARAM > BlockParams
Definition: PlasticOps.hpp:57
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:172
sigmaY
double sigmaY
Yield stress.
Definition: plastic.cpp:174
PlasticOps::CommonData::LAST_PARAM
@ LAST_PARAM
Definition: PlasticOps.hpp:54
PlasticOps::opFactoryDomainRhs
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
Definition: PlasticOps.hpp:500
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:596
PlasticOps::CommonData::mStressPtr
boost::shared_ptr< MatrixDouble > mStressPtr
Definition: PlasticOps.hpp:68
PlasticOps::CommonData::C1_k
@ C1_k
Definition: PlasticOps.hpp:53
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
PlasticOps::CommonData::resC
VectorDouble resC
Definition: PlasticOps.hpp:78
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
PlasticOps::CommonHenkyPtr
boost::shared_ptr< HenckyOps::CommonData > CommonHenkyPtr
Definition: PlasticOps.hpp:244
PlasticOps::CommonData::getPlasticStrainPtr
auto getPlasticStrainPtr()
Definition: PlasticOps.hpp:96
MoFEM::OpCalculateScalarFieldValuesDot
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
Definition: UserDataOperators.hpp:273
OpKPiola
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: seepage.cpp:64
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
bulk_modulus_K
double bulk_modulus_K
Definition: dynamic_first_order_con_law.cpp:96
PlasticOps::CommonData::plasticTau
VectorDouble plasticTau
Definition: PlasticOps.hpp:73
HenckyOps
Definition: HenckyOps.hpp:11
PlasticOps::CommonPlasticPtr
boost::shared_ptr< PlasticOps::CommonData > CommonPlasticPtr
Definition: PlasticOps.hpp:243
PlasticOps::createCommonPlasticOps
auto createCommonPlasticOps(MoFEM::Interface &m_field, std::string block_name, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string u, std::string ep, std::string tau, double scale, Sev sev)
Definition: PlasticOps.hpp:426
PlasticOps::opFactoryDomainLhs
MoFEMErrorCode opFactoryDomainLhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
Definition: PlasticOps.hpp:547
PlasticOps::PlasticityIntegrators::Assembly
Definition: PlasticOps.hpp:176
PlasticOps::CommonData::getPlasticFlowPtr
auto getPlasticFlowPtr()
Definition: PlasticOps.hpp:103
PlasticOps::CommonData::getParamsPtr
auto getParamsPtr()
Definition: PlasticOps.hpp:60
PlasticOpsGeneric.hpp
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
PlasticOps::OpCalculateConstraintsLhs_LogStrain_dUImpl
Definition: PlasticOps.hpp:157
PlasticOps::N
FTensor::Index< 'N', 3 > N
Definition: PlasticOps.hpp:118
PlasticOps::CommonData::getPlasticSurfacePtr
auto getPlasticSurfacePtr()
Definition: PlasticOps.hpp:87
PlasticOps::CommonData::mGradPtr
boost::shared_ptr< MatrixDouble > mGradPtr
Definition: PlasticOps.hpp:66
PlasticOps::Pip
boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > Pip
Definition: PlasticOps.hpp:242
MoFEM::Sev
MoFEM::LogManager::SeverityLevel Sev
Definition: CoreTemplates.hpp:17
PlasticOps::M
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:117
PlasticOps::CommonData::H
@ H
Definition: PlasticOps.hpp:49
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
ContactOps::scale
double scale
Definition: EshelbianContact.hpp:22
PlasticOps::CommonData::getPlasticStrainDotPtr
auto getPlasticStrainDotPtr()
Definition: PlasticOps.hpp:99
PlasticOps::CommonData::blockParams
BlockParams blockParams
Definition: PlasticOps.hpp:58
PlasticOps::OpCalculatePlasticInternalForceLhs_LogStrain_dEPImpl
Definition: PlasticOps.hpp:139
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
DomainEleOp
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
Definition: child_and_parent.cpp:36
OpInternalForcePiola
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: seepage.cpp:66
BlockData
Definition: gel_analysis.cpp:74
PlasticOps::OpCalculatePlasticityImpl
Definition: PlasticOps.hpp:124
AssemblyDomainEleOp
PlasticOps::CommonData::activityData
static std::array< int, 5 > activityData
Definition: PlasticOps.hpp:107
convert.type
type
Definition: convert.py:64
PlasticOps::CommonData::mDPtr
boost::shared_ptr< MatrixDouble > mDPtr
[Common data set externally]
Definition: PlasticOps.hpp:62
EshelbianPlasticity::P
@ P
Definition: EshelbianContact.cpp:193
C1_k
double C1_k
Kinematic hardening.
Definition: plastic.cpp:180
PlasticOps::CommonData::ParamsIndexes
ParamsIndexes
Definition: PlasticOps.hpp:45
PlasticOps::opFactoryDomainReactions
MoFEMErrorCode opFactoryDomainReactions(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
Definition: PlasticOps.hpp:630
PlasticOps::CommonData::mStrainPtr
boost::shared_ptr< MatrixDouble > mStrainPtr
Definition: PlasticOps.hpp:67
PlasticOps::CommonData::plasticTauDot
VectorDouble plasticTauDot
Definition: PlasticOps.hpp:74
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:173
PlasticOps::OpCalculatePlasticFlowLhs_dUImpl
Definition: PlasticOps.hpp:142
OpKCauchy
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Linear elastic problem]
Definition: thermo_elastic.cpp:34
HenckyOps::HenkyIntegrators
Definition: HenckyOps.hpp:569
PlasticOps::CommonData::getPlasticTauPtr
auto getPlasticTauPtr()
Definition: PlasticOps.hpp:90
PlasticOps::CommonData::plasticStrain
MatrixDouble plasticStrain
Definition: PlasticOps.hpp:75
visH
double visH
Viscous hardening.
Definition: plastic.cpp:176
size_symm
constexpr auto size_symm
Definition: plastic.cpp:42
PlasticOps::CommonData::YOUNG_MODULUS
@ YOUNG_MODULUS
Definition: PlasticOps.hpp:46
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
OpInternalForceCauchy
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: thermo_elastic.cpp:38
OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
Definition: operators_tests.cpp:48
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
FTensor::Index
Definition: Index.hpp:23
PlasticOps::CommonData::resFlowDstrainDot
MatrixDouble resFlowDstrainDot
Definition: PlasticOps.hpp:85
PlasticOps::CommonData::QINF
@ QINF
Definition: PlasticOps.hpp:51
PlasticOps::OpCalculatePlasticFlowLhs_dTAUImpl
Definition: PlasticOps.hpp:151
b_iso
double b_iso
Saturation exponent.
Definition: plastic.cpp:179
PlasticOps::CommonData::resCdPlasticStrain
MatrixDouble resCdPlasticStrain
Definition: PlasticOps.hpp:81
Range
DomainEleOp
PlasticOps::OpCalculatePlasticFlowLhs_LogStrain_dUImpl
Definition: PlasticOps.hpp:145
PlasticOps
Definition: PlasticNaturalBCs.hpp:13
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
shear_modulus_G
double shear_modulus_G
Definition: dynamic_first_order_con_law.cpp:97
PlasticOpsLargeStrains.hpp
is_large_strains
PetscBool is_large_strains
Large strains.
Definition: plastic.cpp:167
PlasticOps::OpCalculatePlasticFlowLhs_dEPImpl
Definition: PlasticOps.hpp:148
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
PlasticOps::OpCalculateConstraintsLhs_dEPImpl
Definition: PlasticOps.hpp:160
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
Qinf
double Qinf
Saturation yield stress.
Definition: plastic.cpp:178
PlasticOps::PlasticityIntegrators
Definition: PlasticOps.hpp:165
PlasticOps::CommonData::resFlowDtau
MatrixDouble resFlowDtau
Definition: PlasticOps.hpp:83
PlasticOps::CommonData::resFlow
MatrixDouble resFlow
Definition: PlasticOps.hpp:82
PlasticOps::OpCalculatePlasticSurfaceImpl
Definition: PlasticOps.hpp:121
PlasticOpsSmallStrains.hpp
PlasticOps::CommonData::plasticSurface
VectorDouble plasticSurface
[Common data set externally]
Definition: PlasticOps.hpp:71
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
PlasticOps::CommonData::plasticStrainDot
MatrixDouble plasticStrainDot
Definition: PlasticOps.hpp:76
PlasticOps::CommonData::VIS_H
@ VIS_H
Definition: PlasticOps.hpp:50
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
PlasticOps::CommonData::getPlasticTauDotPtr
auto getPlasticTauDotPtr()
Definition: PlasticOps.hpp:93
PlasticOps::I
FTensor::Index< 'I', 3 > I
[Common data]
Definition: PlasticOps.hpp:115
PlasticOps::OpPlasticStressImpl
Definition: PlasticOps.hpp:127
PlasticOps::CommonData::resFlowDstrain
MatrixDouble resFlowDstrain
Definition: PlasticOps.hpp:84
PlasticOps::CommonData::POISSON_RATIO
@ POISSON_RATIO
Definition: PlasticOps.hpp:47
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:416
PlasticOps::OpCalculatePlasticFlowRhsImpl
Definition: PlasticOps.hpp:130
H
double H
Hardening.
Definition: plastic.cpp:175
PlasticOps::OpCalculatePlasticInternalForceLhs_dEPImpl
Definition: PlasticOps.hpp:136
OpCalculateVectorFieldGradient
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
PlasticOps::CommonData::plasticFlow
MatrixDouble plasticFlow
Definition: PlasticOps.hpp:72
PlasticOpsMonitor.hpp
PlasticOps::CommonData::resCdTau
VectorDouble resCdTau
Definition: PlasticOps.hpp:79
PlasticOps::CommonData::SIGMA_Y
@ SIGMA_Y
Definition: PlasticOps.hpp:48
PlasticOps::OpCalculateConstraintsLhs_dUImpl
Definition: PlasticOps.hpp:154
PlasticOps::CommonData
[Common data]
Definition: PlasticOps.hpp:43
PlasticOps::addMatBlockOps
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, std::string block_name, Pip &pip, boost::shared_ptr< MatrixDouble > mat_D_Ptr, boost::shared_ptr< CommonData::BlockParams > mat_params_ptr, double scale_value, Sev sev)
Definition: PlasticOps.hpp:248
Generated by Doxygen 1.8.17 and hosted at University of Glasgow