v0.15.0
Loading...
Searching...
No Matches
ADOLCPlasticity.hpp
Go to the documentation of this file.
1/** \file ADOLCPlasticity.hpp
2 * \ingroup adoc_plasticity
3 * \example ADOLCPlasticity.hpp
4 *
5 * \brief Operators and data structures for small strain plasticity
6 *
7 * \defgroup adoc_plasticity ADOL-C plasticity
8 * \ingroup user_modules
9 * \defgroup user_modules User modules
10 *
11**/
12
13#ifndef __ADOLCPLASTICITY_HPP_
14#define __ADOLCPLASTICITY_HPP_
15
16#ifndef WITH_ADOL_C
17#error "MoFEM need to be compiled with ADOL-C"
18#endif
19
20/**
21 * \ingroup adoc_plasticity
22 *
23 */
24namespace ADOLCPlasticity {
25
27
28/**
29 * \brief Op convert Vight strain vector to strain tensor
30*/
32 return FTensor::Dg<double, 3, 6>{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
33 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0,
34 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
35 0.0, 0.0, 0.5, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0};
36};
37
38/**
39 * \brief Op convert strain tensor to Vight strain vector
40*/
42 return FTensor::Dg<double, 3, 6>{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
43 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
44 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
45 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0};
46};
47
48/**
49 * \brief Op convert Vight stress vector to stress tensor
50*/
52 return FTensor::Dg<double, 3, 6>{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
53 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
54 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
55 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0};
56};
57
58/** \brief common data used by volume elements
59 * \ingroup nonlinear_elastic_elem
60 */
61struct CommonData : boost::enable_shared_from_this<CommonData> {
62
63 MatrixDouble activeVariablesW;
64 MatrixDouble gradientW;
65
66 inline auto getFTensor1StressAtGaussPts(int gg = 0) {
68 &gradientW(gg, 7),
69 &gradientW(gg, 8),
70 &gradientW(gg, 9),
71 &gradientW(gg, 10),
72 &gradientW(gg, 11),
73 static_cast<int>(gradientW.size2())};
74 };
75
76 inline auto getFTensor1PlasticStrainAtGaussPts(int gg = 0) {
78 &activeVariablesW(gg, 0),
79 &activeVariablesW(gg, 1),
80 &activeVariablesW(gg, 2),
81 &activeVariablesW(gg, 3),
82 &activeVariablesW(gg, 4),
83 &activeVariablesW(gg, 5),
84 static_cast<int>(activeVariablesW.size2())};
85 };
86
87 MatrixDouble gradAtGaussPts;
89
90 inline auto getPlasticStrain(int gg = 0) {
91 return getVectorAdaptor(&(activeVariablesW(gg, 0)), 6);
92 }
93
94 inline auto getInternalVariables(int gg = 0) {
95 return getVectorAdaptor(&(activeVariablesW(gg, 12)),
96 activeVariablesW.size2() - 12);
97 }
98
99 vector<double> deltaGamma; //< Lagrange plastic multiplier
104
105 inline auto getGradAtGaussPtsPtr() {
106 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &gradAtGaussPts);
107 }
108
109 inline auto getMatTangentPtr() {
110 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &materialTangentOperator);
111 }
112
113 bool bBar = true;
114
117 PetscBool b_bar = bBar ? PETSC_TRUE : PETSC_FALSE;
118 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, PETSC_NULLPTR, "-b_bar", &b_bar,
119 PETSC_NULLPTR);
120 bBar = b_bar;
122 }
123
124 boost::shared_ptr<MatrixDouble> getStressMatrixPtr() {
125 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &stressMatrix);
126 }
127 boost::shared_ptr<MatrixDouble> getPlasticStrainMatrixPtr() {
128 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
130 }
131
132 MatrixDouble stressMatrix;
134
136
137 CHK_THROW_MESSAGE(getDefaultMaterialParameters(), "get parameters failed");
138 }
139};
140
141/** \brief Closest point projection algorithm
142* \ingroup small_strain_plasticity
143
144*/
146
148 boost::function<int(int, int, int)> integrationRule = [](int, int, int p) {
149 return 2 * (p - 1);
150 };
151
152 VectorDouble internalVariables0;
153 VectorDouble plasticStrain0;
154
155 inline VectorAdaptor getPlasticStrain() {
156 return getVectorAdaptor(&(activeVariablesW[0]), 6);
157 }
158 inline VectorAdaptor getTotalStrain() {
159 return getVectorAdaptor(&(activeVariablesW[6]), 6);
160 }
161 inline VectorAdaptor getInternalVariables() {
162 return getVectorAdaptor(&(activeVariablesW[12]), nbInternalVariables);
163 }
164 inline VectorAdaptor getActiveVariablesYH() {
165 return getVectorAdaptor(&(gradientW[6]), 6 + nbInternalVariables);
166 }
167 inline VectorAdaptor getStress() {
168 return getVectorAdaptor(&(gradientW[6]), 6);
169 }
170 inline VectorAdaptor getInternalFluxes() {
171 return getVectorAdaptor(&(gradientW[12]), nbInternalVariables);
172 }
173
174 enum TypesTags { W = 0, Y, H, LAST_TAPE }; //< W - energy, Y - yield, H - flow
175 std::array<int, LAST_TAPE> tapesTags; //< Tapes nmbers
176
178
179 VectorAdaptor activeVariablesW;
180 VectorAdaptor gradientW;
181
182 double w;
183 double y;
184 double h;
185
186 double deltaGamma; // Increment of plastic multiplier
187 MatrixDouble Ep, Cp, Cep;
188 ublas::symmetric_matrix<double, ublas::lower> C, D;
189
190 PetscBool implHessianW;
191
192 /**
193 * \brief Record strain energy
194 */
195 MoFEMErrorCode recordW();
196
197 /**
198 * \brief Record yield function
199 */
200 MoFEMErrorCode recordY();
201
202 /**
203 * \brief Record flow potential
204 */
205 MoFEMErrorCode recordH();
206
207 MatrixDouble hessianW; //< Hessian of energy
208 MoFEMErrorCode playW();
209 MoFEMErrorCode playW_NoHessian();
210
211 VectorDouble gradientY; //< Gradient of yield function
212 MoFEMErrorCode playY();
213 MoFEMErrorCode playY_NoGradient();
214
215 VectorDouble gradientH; //< Gradient of flow potential
216 MatrixDouble hessianH; //< Hessian of flow potential
217 MoFEMErrorCode playH();
218 MoFEMErrorCode playH_NoHessian();
219
220 MoFEMErrorCode createMatAVecR(); //< For integration point SNES solver
221
222 MoFEMErrorCode evaluatePotentials(); //< Evaluate potentials
223
224 MoFEMErrorCode
225 playPotentials(); //< Play potentials from recorded ADOl-C tapes
226
227 MoFEMErrorCode
228 playPotentials_NoHessian(); //< Play potentials from recorded ADOl-C tapes
229
230 MoFEMErrorCode calculateR(Vec R); //< Calculate residual
231
232 MoFEMErrorCode calculateA(); //< Calculate tangent matrix
233
234 /**
235 * \brief Function executed by nested SENES at evaluationg residual
236 */
237 friend MoFEMErrorCode ADOLCPlasticityRes(SNES, Vec, Vec, void *ctx);
238
239 /**
240 * \brief Function executed by nested SENES at evaluationg tangent matrix
241 */
242 friend MoFEMErrorCode ADOLCPlasticityJac(SNES, Vec, Mat, Mat, void *ctx);
243
244 /**
245 * \brief Create nested snes
246 */
247 MoFEMErrorCode snesCreate();
248
249 /**
250 * \brief Solve nonlinear system of equations to find stress, internal
251 * fluxes, and Lagrange plastic multiplier
252 */
253 MoFEMErrorCode solveClosestProjection();
254
255 /**
256 * \brief Calculate consistent tangent matrix
257 */
258 MoFEMErrorCode consistentTangent();
259
260 /**
261 * \brief Record tapes
262 */
263 MoFEMErrorCode recordTapes();
264
265 /**
266 * \brief Get model parameters from blocks
267 */
268 virtual MoFEMErrorCode
270 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
271 std::string block_name, Sev sev) {
272 return 0;
273 }
274
275 /**
276 * \brief Set parameters for ADOL-C of constitutive relations
277 *
278 * \param tag [in] - tag of the tape
279 * \param recalculate_elastic_tangent [out] - if setParam set it to true,
280 * tangent matrix for elastic domain should be recalculated
281 */
282 virtual MoFEMErrorCode setParams(short tag,
283 bool &recalculate_elastic_tangent) {
284 return 0;
285 }
286
287 /**
288 * \brief Set Hemholtz energy
289 */
290 virtual MoFEMErrorCode freeHemholtzEnergy() = 0;
291
292 /**
293 * \brief Set yield function
294 */
295 virtual MoFEMErrorCode yieldFunction() = 0;
296
297 /**
298 * \brief Set flow potential
299 */
300 virtual MoFEMErrorCode flowPotential() = 0;
301
302 virtual MoFEMErrorCode codedHessianW(vector<double *>) {
303 return 0;
304 }
305
306 // protected:
307
308 MatrixDouble partialWStrainPlasticStrain; //< Partial derivative of free energy
309 // with respect to plastic strain
310 VectorAdaptor partialYSigma; //< Partial derivative of yield function with
311 // respect to stress
312 VectorAdaptor partialYFlux; //< Partial derivative of yield function with
313 // respect to internal fluxes
314 VectorAdaptor partialHSigma; //< Partial derivative of flow potential with
315 // respect to stress
316 VectorAdaptor partialHFlux; //< Partial derivative of flow potential with
317 // respect to internal fluxes
318
319 /// Second partial derivative of flow potential with respect to stresses
320 /// and internal
321 ublas::symmetric_matrix<double, ublas::lower> partial2HSigma;
322 /// Second partial derivative of flow potential with respect to internal
323 /// fluxes
324 ublas::symmetric_matrix<double, ublas::lower> partial2HFlux;
325 /// Mixed decond partial derivative of flow potential with respect to stresses
326 /// and internal fluxes
327 MatrixDouble partial2HSigmaFlux;
328
329 SmartPetscObj<Mat> A; //< Nested SNES tangent matrix
330 SmartPetscObj<Vec> R; //< Nested SNES residual
331 SmartPetscObj<Vec> Chi; //< Nested SNES unknown vector
332 SmartPetscObj<Vec> dChi; //< Nested SNES increment vector
333 ublas::matrix<double, ublas::column_major> dataA;
334 SmartPetscObj<SNES> sNes; //< Nested SNES solver
335
336 ublas::vector<adouble> a_sTrain;
337 ublas::vector<adouble> a_plasticStrain;
338 ublas::vector<adouble> a_internalVariables;
339 ublas::vector<adouble> a_sTress, a_internalFluxes;
343};
344
345/**
346 * \brief Internal SNES function used at integration points to calulate stress
347*/
348MoFEMErrorCode ADOLCPlasticityRes(SNES snes, Vec chi, Vec r, void *ctx);
349
350/**
351 * \brief Internal SNES function used at integration points to calulate
352 * tangent matrix
353*/
354MoFEMErrorCode ADOLCPlasticityJac(SNES, Vec, Mat, Mat, void *ctx);
355
356/**
357 * \brief Get opreator to calulate stress
358*/
359template <int DIM, StrainType STRAIN>
360ForcesAndSourcesCore::UserDataOperator *getRawPtrOpCalculateStress(
361 MoFEM::Interface &m_field, boost::shared_ptr<CommonData> common_data_ptr,
362 boost::shared_ptr<ClosestPointProjection> cp_ptr, bool calc_lhs);
363
364template <>
365ForcesAndSourcesCore::UserDataOperator *getRawPtrOpCalculateStress<3, LARGE_STRAIN>(
366 MoFEM::Interface &m_field, boost::shared_ptr<CommonData> common_data_ptr,
367 boost::shared_ptr<ClosestPointProjection> cp_ptr, bool calc_lhs);
368
369template <>
370ForcesAndSourcesCore::UserDataOperator *getRawPtrOpCalculateStress<3, SMALL_STRAIN>(
371 MoFEM::Interface &m_field, boost::shared_ptr<CommonData> common_data_ptr,
372 boost::shared_ptr<ClosestPointProjection> cp_ptr, bool calc_lhs);
373
374template <>
375ForcesAndSourcesCore::UserDataOperator *getRawPtrOpCalculateStress<2, LARGE_STRAIN>(
376 MoFEM::Interface &m_field, boost::shared_ptr<CommonData> common_data_ptr,
377 boost::shared_ptr<ClosestPointProjection> cp_ptr, bool calc_lhs);
378
379template <>
380ForcesAndSourcesCore::UserDataOperator *getRawPtrOpCalculateStress<2, SMALL_STRAIN>(
381 MoFEM::Interface &m_field, boost::shared_ptr<CommonData> common_data_ptr,
382 boost::shared_ptr<ClosestPointProjection> cp_ptr, bool calc_lhs);
383
384/**
385 * \brief Assemble right hand side
386 *
387 * \param field_name [in] - name of field
388 * \param common_data_ptr [in] - shared pointer to common data
389 *
390 * \ingroup small_strain_plasticity
391 */
392template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
394
395/**
396 * \brief Assemble left hand side
397 *
398 * \param field_name [in] - name of field
399 * \param common_data_ptr [in] - shared pointer to common data
400 *
401 * \ingroup small_strain_plasticity
402 */
403template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
405
406template <typename DomainEleOp> struct ADOLCPlasticityIntegrators {
407
408 template <AssemblyType A> struct Assembly {
409
411 typename FormsIntegrators<DomainEleOp>::template Assembly<A>::OpBase;
412
413 template <int DIM, IntegrationType I>
415
416 template <int DIM, IntegrationType I>
418 };
419};
420
421using Pip = boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator>;
422
423/**
424 * @brief Assemble the left hand side, i.e. tangent matrix
425 * @ingroup adoc_plasticity
426 *
427 * @tparam DIM dimension of the problem
428 * @tparam A assembly type
429 * @tparam I integration type
430 * @tparam DomainEleOp operator type
431 * @param m_field
432 * @param field_name
433 * @param pip
434 * @param block_name esh block name caring material parameters
435 * @param common_data_ptr
436 * @param cp_ptr
437 * @return MoFEMErrorCode
438 */
439template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
440MoFEMErrorCode opFactoryDomainRhs(
441 MoFEM::Interface &m_field, string field_name, Pip &pip,
442 std::string block_name, boost::shared_ptr<CommonData> common_data_ptr,
443 boost::shared_ptr<ClosestPointProjection> cp_ptr, Sev sev = Sev::inform) {
446 CHKERR cp_ptr->addMatBlockOps(m_field, pip, block_name, sev);
447 pip.push_back(
448 getRawPtrOpCalculateStress<DIM, SMALL_STRAIN>(m_field, common_data_ptr, cp_ptr, false));
449 pip.push_back(new typename P::template Assembly<A>::template OpRhs<DIM, I>(
450 field_name, common_data_ptr));
452}
453/**
454 * @brief Assemble the left hand side, i.e. tangent matrix
455 * @ingroup adoc_plasticity
456 *
457 * @tparam DIM dimension of the problem
458 * @tparam A assembly type
459 * @tparam I integration type
460 * @tparam DomainEleOp operator type
461 * @param m_field
462 * @param field_name
463 * @param pip
464 * @param block_name esh block name caring material parameters
465 * @param common_data_ptr
466 * @param cp_ptr
467 * @return MoFEMErrorCode
468 */
469template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
470MoFEMErrorCode
472 std::string block_name,
473 boost::shared_ptr<CommonData> common_data_ptr,
474 boost::shared_ptr<ClosestPointProjection> cp_ptr) {
477 CHKERR cp_ptr->addMatBlockOps(m_field, pip, block_name, Sev::noisy);
478 pip.push_back(
479 getRawPtrOpCalculateStress<DIM, SMALL_STRAIN>(m_field, common_data_ptr, cp_ptr, true));
480 pip.push_back(new typename P::template Assembly<A>::template OpLhs<DIM, I>(
481 field_name, common_data_ptr));
483}
484
485/**
486 * @brief Push operators to update history variables
487 * @ingroup adoc_plasticity
488 *
489 * @tparam DIM dimension of the problem
490 * @param m_field core interface
491 * @param pip
492 * @param block_name mesh block name caring material parameters
493 * @param common_data_ptr
494 * @param cp_ptr
495 * @return MoFEMErrorCode
496 */
497template <int DIM>
498MoFEMErrorCode
500 std::string block_name,
501 boost::shared_ptr<CommonData> common_data_ptr,
502 boost::shared_ptr<ClosestPointProjection> cp_ptr);
503
504/**
505 * @copydoc ADOLCPlasticity::opFactoryDomainUpdate
506 */
507template <>
508MoFEMErrorCode
510 std::string block_name,
511 boost::shared_ptr<CommonData> common_data_ptr,
512 boost::shared_ptr<ClosestPointProjection> cp_ptr);
513
514/**
515 * @copydoc ADOLCPlasticity::opFactoryDomainUpdate
516 */
517template <>
518MoFEMErrorCode
520 std::string block_name,
521 boost::shared_ptr<CommonData> common_data_ptr,
522 boost::shared_ptr<ClosestPointProjection> cp_ptr);
523
524
525
526/**
527 * \brief Update internal fluxes (update history variables)
528 * \ingroup adoc_plasticity
529 */
530struct TSUpdate {
531 virtual MoFEMErrorCode postProcess(TS ts) = 0;
532};
533boost::shared_ptr<TSUpdate> createTSUpdate(std::string fe_name,
534 boost::shared_ptr<FEMethod> fe_ptr);
535
536/**
537 * \brief Calculate tensorial base functions. Apply bBar method when needed
538 */
539struct MakeB {
540
542 getFTensor2SymmetricDiffBase(DataForcesAndSourcesCore::EntData &data,
543 MatrixDouble &storage, const bool b_bar,
544 const int nb_integration_pts, double *w_ptr,
546
548 getFTensor2SymmetricDiffBase(DataForcesAndSourcesCore::EntData &data,
549 MatrixDouble &storage, const bool b_bar,
550 const int nb_integration_pts, double *w_ptr,
552};
553
554
555template <int DIM, typename AssemblyDomainEleOp>
557 OpRhsImpl(std::string field_name,
558 boost::shared_ptr<CommonData> common_data_ptr);
559 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
560
561private:
562 boost::shared_ptr<CommonData> commonDataPtr;
563 MatrixDouble baseStorage; ///< Store tensorial base functions
564};
565
566template <int DIM, typename AssemblyDomainEleOp>
569 boost::shared_ptr<CommonData> common_data_ptr);
570 MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data,
571 DataForcesAndSourcesCore::EntData &col_data);
572
573protected:
574 boost::shared_ptr<CommonData> commonDataPtr;
575 MatrixDouble baseRowStorage; ///< Store tensorial base functions
576 MatrixDouble baseColStorage; ///< Store tensorial base functions
577};
578
579template <int DIM, typename AssemblyDomainEleOp>
581 string field_name, boost::shared_ptr<CommonData> common_data_ptr)
583 commonDataPtr(common_data_ptr) {}
584
585//! [OpRhsImpl integrate]
586template <int DIM, typename AssemblyDomainEleOp>
588 EntitiesFieldData::EntData &data) {
590 using OP = AssemblyDomainEleOp;
591 FTensor::Index<'i', 3> i;
592 FTensor::Index<'j', 3> j;
593 double *w_ptr = &(OP::getGaussPts()(DIM, 0));
594 // Calulate tensorial base functions. Apply bBar method when needed
596 data, baseStorage, commonDataPtr->bBar, OP::getGaussPts().size2(), w_ptr,
598
599 auto t_stress = getFTensor2SymmetricFromMat<3>(commonDataPtr->stressMatrix);
600
601 const auto vol = OP::getMeasure();
602 auto t_w = OP::getFTensor0IntegrationWeight();
603 for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
604 double alpha = vol * t_w;
605 for (int bb = 0; bb != OP::nbRows; ++bb) {
606 OP::locF[bb] += alpha * (t_stress(i, j) * t_diff(i, j));
607 ++t_diff;
608 }
609 ++t_w;
610 ++t_stress;
611 }
613}
614//! [OpRhsImpl integrate]
615
616template <int DIM, typename AssemblyDomainEleOp>
618 string field_name, boost::shared_ptr<CommonData> common_data_ptr)
620 AssemblyDomainEleOp::OPROWCOL),
621 commonDataPtr(common_data_ptr) {
622 this->sYmm = false; // It has to be false for not-associtive plasticity
623}
624
625template <int DIM, typename AssemblyDomainEleOp>
627 DataForcesAndSourcesCore::EntData &row_data,
628 DataForcesAndSourcesCore::EntData &col_data) {
630 using OP = AssemblyDomainEleOp;
631
632 double *w_ptr = &(OP::getGaussPts()(DIM, 0));
633 auto t_diff_row = MakeB::getFTensor2SymmetricDiffBase(
634 row_data, baseRowStorage, commonDataPtr->bBar, OP::getGaussPts().size2(),
635 w_ptr, FTensor::Number<DIM>());
636 auto t_diff_col = MakeB::getFTensor2SymmetricDiffBase(
637 col_data, baseColStorage, commonDataPtr->bBar, OP::getGaussPts().size2(),
638 w_ptr, FTensor::Number<DIM>());
639
640 auto get_ftensor2_symmetric = [&](auto &storage, const auto gg,
641 const auto rr) {
643 &storage(gg, 6 * rr + 0), &storage(gg, 6 * rr + 1),
644 &storage(gg, 6 * rr + 2), &storage(gg, 6 * rr + 3),
645 &storage(gg, 6 * rr + 4), &storage(gg, 6 * rr + 5)};
646 };
647
648 FTensor::Index<'i', 3> i;
649 FTensor::Index<'j', 3> j;
650 FTensor::Index<'k', 3> k;
651 FTensor::Index<'l', 3> l;
652
653 auto t_Cep =
654 getFTensor4DdgFromMat<3, 3, 1>(commonDataPtr->materialTangentOperator);
655 const auto vol = OP::getMeasure();
656 auto t_w = OP::getFTensor0IntegrationWeight();
657 for (int gg = 0; gg != OP::nbIntegrationPts; ++gg) {
658 const double alpha = vol * t_w;
659 ++t_w;
660
661 for (auto rr = 0; rr != OP::nbRows; ++rr) {
663 t_rowCep(k, l) = t_Cep(i, j, k, l) * t_diff_row(i, j);
664 auto t_diff_col = get_ftensor2_symmetric(baseColStorage, gg, 0);
665 for (auto cc = 0; cc != OP::nbCols; ++cc) {
666 OP::locMat(rr, cc) += alpha * (t_rowCep(k, l) * t_diff_col(k, l));
667 ++t_diff_col;
668 }
669 ++t_diff_row;
670 }
671
672 ++t_Cep;
673 }
675}
676
677} // namespace ADOLCPlasticity
678
679#endif //__ADOLCPLASTICITY_HPP_
680
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
MoFEMErrorCode opFactoryDomainLhs(MoFEM::Interface &m_field, string field_name, Pip &pip, std::string block_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr)
Assemble the left hand side, i.e. tangent matrix.
MoFEMErrorCode opFactoryDomainUpdate(MoFEM::Interface &m_field, Pip &pip, std::string block_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr)
Push operators to update history variables.
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, string field_name, Pip &pip, std::string block_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr, Sev sev=Sev::inform)
Assemble the left hand side, i.e. tangent matrix.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
ForcesAndSourcesCore::UserDataOperator * getRawPtrOpCalculateStress< 3, SMALL_STRAIN >(MoFEM::Interface &m_field, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr, bool calc_lhs)
ForcesAndSourcesCore::UserDataOperator * getRawPtrOpCalculateStress< 2, SMALL_STRAIN >(MoFEM::Interface &m_field, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr, bool calc_lhs)
FTensor::Dg< double, 3, 6 > strain_to_voight_op()
Op convert strain tensor to Vight strain vector.
FTensor::Dg< double, 3, 6 > voight_to_stress_op()
Op convert Vight stress vector to stress tensor.
MoFEMErrorCode opFactoryDomainUpdate< 2 >(MoFEM::Interface &m_field, Pip &pip, std::string block_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr)
Push operators to update history variables.
MoFEMErrorCode ADOLCPlasticityJac(SNES, Vec, Mat, Mat, void *ctx)
Internal SNES function used at integration points to calulate tangent matrix.
FTensor::Dg< double, 3, 6 > voight_to_strain_op()
Op convert Vight strain vector to strain tensor.
MoFEMErrorCode ADOLCPlasticityRes(SNES snes, Vec chi, Vec r, void *ctx)
Internal SNES function used at integration points to calulate stress.
ForcesAndSourcesCore::UserDataOperator * getRawPtrOpCalculateStress(MoFEM::Interface &m_field, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr, bool calc_lhs)
Get opreator to calulate stress.
boost::shared_ptr< TSUpdate > createTSUpdate(std::string fe_name, boost::shared_ptr< FEMethod > fe_ptr)
ForcesAndSourcesCore::UserDataOperator * getRawPtrOpCalculateStress< 2, LARGE_STRAIN >(MoFEM::Interface &m_field, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr, bool calc_lhs)
ForcesAndSourcesCore::UserDataOperator * getRawPtrOpCalculateStress< 3, LARGE_STRAIN >(MoFEM::Interface &m_field, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr, bool calc_lhs)
boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > Pip
MoFEMErrorCode opFactoryDomainUpdate< 3 >(MoFEM::Interface &m_field, Pip &pip, std::string block_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr)
Push operators to update history variables.
constexpr auto field_name
Closest point projection algorithm.
std::array< int, LAST_TAPE > tapesTags
MoFEMErrorCode solveClosestProjection()
Solve nonlinear system of equations to find stress, internal fluxes, and Lagrange plastic multiplier.
MoFEMErrorCode recordH()
Record flow potential.
friend MoFEMErrorCode ADOLCPlasticityRes(SNES, Vec, Vec, void *ctx)
Function executed by nested SENES at evaluationg residual.
MoFEMErrorCode recordY()
Record yield function.
boost::function< int(int, int, int)> integrationRule
virtual MoFEMErrorCode flowPotential()=0
Set flow potential.
ublas::symmetric_matrix< double, ublas::lower > partial2HFlux
MoFEMErrorCode recordW()
Record strain energy.
virtual MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string block_name, Sev sev)
Get model parameters from blocks.
ublas::symmetric_matrix< double, ublas::lower > D
ublas::symmetric_matrix< double, ublas::lower > C
virtual MoFEMErrorCode setParams(short tag, bool &recalculate_elastic_tangent)
Set parameters for ADOL-C of constitutive relations.
virtual MoFEMErrorCode codedHessianW(vector< double * >)
virtual MoFEMErrorCode yieldFunction()=0
Set yield function.
MoFEMErrorCode consistentTangent()
Calculate consistent tangent matrix.
ublas::symmetric_matrix< double, ublas::lower > partial2HSigma
friend MoFEMErrorCode ADOLCPlasticityJac(SNES, Vec, Mat, Mat, void *ctx)
Function executed by nested SENES at evaluationg tangent matrix.
ublas::matrix< double, ublas::column_major > dataA
virtual MoFEMErrorCode freeHemholtzEnergy()=0
Set Hemholtz energy.
MoFEMErrorCode snesCreate()
Create nested snes.
common data used by volume elements
MoFEMErrorCode getDefaultMaterialParameters()
auto getFTensor1StressAtGaussPts(int gg=0)
auto getFTensor1PlasticStrainAtGaussPts(int gg=0)
boost::shared_ptr< MatrixDouble > getStressMatrixPtr()
boost::shared_ptr< MatrixDouble > getPlasticStrainMatrixPtr()
Calculate tensorial base functions. Apply bBar method when needed.
static FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 6 >, 3 > getFTensor2SymmetricDiffBase(DataForcesAndSourcesCore::EntData &data, MatrixDouble &storage, const bool b_bar, const int nb_integration_pts, double *w_ptr, FTensor::Number< 2 >)
MatrixDouble baseRowStorage
Store tensorial base functions.
OpLhsImpl(std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
MatrixDouble baseColStorage
Store tensorial base functions.
Assemble left hand side.
MatrixDouble baseStorage
Store tensorial base functions.
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Assemble right hand side.
Update internal fluxes (update history variables)
virtual MoFEMErrorCode postProcess(TS ts)=0
Deprecated interface functions.
FormsIntegrators< DomainEleOp >::Assembly< A >::OpBase AssemblyDomainEleOp