v0.13.1
NonLinearElasticElement.hpp
Go to the documentation of this file.
1/** \file NonLinearElasticElement.hpp
2 * \ingroup nonlinear_elastic_elem
3 * \brief Operators and data structures for non-linear elastic analysis
4 *
5 * Implementation of nonlinear elastic element.
6 */
7
8
9
10#ifndef __NONLINEAR_ELASTIC_HPP
11#define __NONLINEAR_ELASTIC_HPP
12
13#ifndef WITH_ADOL_C
14#error "MoFEM need to be compiled with ADOL-C"
15#endif
16
17/** \brief structure grouping operators and data used for calculation of
18 * nonlinear elastic element \ingroup nonlinear_elastic_elem
19 *
20 * In order to assemble matrices and right hand vectors, the loops over
21 * elements, entities over that elements and finally loop over integration
22 * points are executed.
23 *
24 * Following implementation separate those three categories of loops and to each
25 * loop attach operator.
26 *
27 */
29
30 /// \brief definition of volume element
32
33 Mat A;
35
37
39 virtual ~MyVolumeFE() = default;
40
41 /** \brief it is used to calculate nb. of Gauss integration points
42 *
43 * for more details pleas look
44 * Reference:
45 *
46 * Albert Nijenhuis, Herbert Wilf,
47 * Combinatorial Algorithms for Computers and Calculators,
48 * Second Edition,
49 * Academic Press, 1978,
50 * ISBN: 0-12-519260-6,
51 * LC: QA164.N54.
52 *
53 * More details about algorithm
54 * http://people.sc.fsu.edu/~jburkardt/cpp_src/gm_rule/gm_rule.html
55 **/
56 int getRule(int order);
57
58 SmartPetscObj<Vec> V;
59 double eNergy;
60
63 };
64
65 MyVolumeFE feRhs; ///< calculate right hand side for tetrahedral elements
66 MyVolumeFE &getLoopFeRhs() { return feRhs; } ///< get rhs volume element
67 MyVolumeFE feLhs; //< calculate left hand side for tetrahedral elements
68 MyVolumeFE &getLoopFeLhs() { return feLhs; } ///< get lhs volume element
69
70 MyVolumeFE feEnergy; ///< calculate elastic energy
71 MyVolumeFE &getLoopFeEnergy() { return feEnergy; } ///< get energy fe
72
74 short int tAg;
75
76 NonlinearElasticElement(MoFEM::Interface &m_field, short int tag);
77 virtual ~NonlinearElasticElement() = default;
78
79 template <typename TYPE> struct FunctionsToCalculatePiolaKirchhoffI;
80
81 /** \brief data for calculation het conductivity and heat capacity elements
82 * \ingroup nonlinear_elastic_elem
83 */
84 struct BlockData {
85 int iD;
86 double E;
87 double PoissonRatio;
88 // Eberlein Fibres stiffness properties
89 double k1, k2;
90 Range tEts; ///< constrains elements in block set
91 boost::shared_ptr<FunctionsToCalculatePiolaKirchhoffI<adouble>>
93 boost::shared_ptr<FunctionsToCalculatePiolaKirchhoffI<double>>
97 };
98
99 std::map<int, BlockData>
100 setOfBlocks; ///< maps block set id with appropriate BlockData
101
102 /** \brief common data used by volume elements
103 * \ingroup nonlinear_elastic_elem
104 */
105 struct CommonData {
106
107 std::map<std::string, std::vector<VectorDouble>> dataAtGaussPts;
108 std::map<std::string, std::vector<MatrixDouble>> gradAtGaussPts;
111 std::vector<MatrixDouble3by3> sTress;
112 std::vector<MatrixDouble>
113 jacStress; ///< this is simply material tangent operator
114
115 // This part can be used to calculate stress directly from potential
116
117 std::vector<double> eNergy;
118 std::vector<VectorDouble> jacEnergy;
119 std::vector<VectorDouble> hessianEnergy;
120 };
122
126
127 /** \brief Implementation of elastic (non-linear) St. Kirchhoff equation
128 * \ingroup nonlinear_elastic_elem
129 */
130 template <typename TYPE> struct FunctionsToCalculatePiolaKirchhoffI {
131
135
142 }
143
145
146 double lambda, mu;
147 MatrixBoundedArray<TYPE, 9> F, C, E, S, invF, P, sIGma, h, H, invH,
151
153
154 int gG; ///< Gauss point number
155 CommonData *commonDataPtr; ///< common data shared between entities (f.e.
156 ///< field values at Gauss pts.)
158 *opPtr; ///< pointer to finite element tetrahedral operator
159
162 t_C(i, j) = t_F(k, i) * t_F(k, j);
164 }
165
168 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
169 t_E(i, j) = 0.5 * (t_C(i, j) - t_kd(i, j));
171 }
172
173 // St. Venant–Kirchhoff Material
176 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
177 t_S(i, j) = (2 * mu) * t_E(i, j) + (lambda * t_E(k, k)) * t_kd(i, j);
179 }
180
181 /** \brief Function overload to implement user material
182 *
183
184 * Calculation of Piola Kirchhoff I is implemented by user. Tangent matrix
185 * user implemented physical equation is calculated using automatic
186 * differentiation.
187
188 * \f$\mathbf{S} =
189 \lambda\textrm{tr}[\mathbf{E}]\mathbf{I}+2\mu\mathbf{E}\f$
190
191 * Notes: <br>
192 * Number of actual Gauss point is accessed from variable gG. <br>
193 * Access to operator data structures is available by variable opPtr. <br>
194 * Access to common data is by commonDataPtr. <br>
195
196 * \param block_data used to give access to material parameters
197 * \param fe_ptr pointer to element data structures
198
199 For details look to: <br>
200 NONLINEAR CONTINUUM MECHANICS FOR FINITE ELEMENT ANALYSIS, Javier Bonet,
201 Richard D. Wood
202
203 */
205 const BlockData block_data,
206 boost::shared_ptr<const NumeredEntFiniteElement> fe_ptr) {
208 lambda = LAMBDA(block_data.E, block_data.PoissonRatio);
209 mu = MU(block_data.E, block_data.PoissonRatio);
213 t_P(i, j) = t_F(i, k) * t_S(k, j);
215 }
216
217 /** \brief Function overload to implement user material
218 *
219
220 * Calculation of Piola Kirchhoff I is implemented by user. Tangent matrix
221 * user implemented physical equation is calculated using automatic
222 * differentiation.
223
224 * \f$\mathbf{S} =
225 \lambda\textrm{tr}[\mathbf{E}]\mathbf{I}+2\mu\mathbf{E}\f$
226
227 * Notes: <br>
228 * Number of actual Gauss point is accessed from variable gG. <br>
229 * Access to operator data structures is available by variable opPtr. <br>
230 * Access to common data is by commonDataPtr. <br>
231
232 * \param block_data used to give access to material parameters
233 * \param fe_ptr pointer to element data structures
234
235 For details look to: <br>
236 NONLINEAR CONTINUUM MECHANICS FOR FINITE ELEMENT ANALYSIS, Javier Bonet,
237 Richard D. Wood
238
239 */
241 const BlockData block_data,
242 boost::shared_ptr<const NumeredEntFiniteElement> fe_ptr) {
244 sigmaCauchy.resize(3, 3);
245 t_sigmaCauchy(i, j) = t_P(i, k) * t_F(j, k);
248 }
249
250 /**
251 * \brief add additional active variables
252 *
253 * \note This member function if used should be implement by template member
254 * function Specialization, different implementation needed for TYPE=double
255 * or TYPE=adouble
256 *
257 * More complex physical models depend on gradient of defamation and some
258 * additional variables. For example can depend on temperature. This
259 * function adds additional independent variables to the model.
260 *
261 * @param nb_active_variables number of active variables
262 * @return error code
263 */
264 virtual MoFEMErrorCode setUserActiveVariables(int &nb_active_variables) {
267 }
268
269 /**
270 * \brief Add additional independent variables
271 * More complex physical models depend on gradient of defamation and some
272 * additional variables. For example can depend on temperature. This
273 * function adds additional independent variables to the model.
274 *
275 * /note First 9 elements are reserved for gradient of deformation.
276 * @param activeVariables vector of deepened variables, values after index
277 * 9 should be add.
278 *
279 * @return error code
280 */
281 virtual MoFEMErrorCode
285 }
286
287 /** \brief Calculate elastic energy density
288 *
289 * \f[\Psi =
290 * \frac{1}{2}\lambda(\textrm{tr}[\mathbf{E}])^2+\mu\mathbf{E}:\mathbf{E}\f]
291 */
293 const BlockData block_data,
294 boost::shared_ptr<const NumeredEntFiniteElement> fe_ptr) {
296 lambda = LAMBDA(block_data.E, block_data.PoissonRatio);
297 mu = MU(block_data.E, block_data.PoissonRatio);
300 TYPE trace = 0;
301 eNergy = 0;
302 for (int ii = 0; ii < 3; ii++) {
303 trace += E(ii, ii);
304 for (int jj = 0; jj < 3; jj++) {
305 TYPE e = E(ii, jj);
306 eNergy += mu * e * e;
307 }
308 }
309 eNergy += 0.5 * lambda * trace * trace;
311 }
312
313 /** \brief Calculate Eshelby stress
314 */
316 const BlockData block_data,
317 boost::shared_ptr<const NumeredEntFiniteElement> fe_ptr) {
319 CHKERR calculateP_PiolaKirchhoffI(block_data, fe_ptr);
320 CHKERR calculateElasticEnergy(block_data, fe_ptr);
321 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
322 t_sIGma(i, j) = t_kd(i, j) * eNergy - t_F(k, i) * t_P(k, j);
324 }
325
326 /** \brief Do operations when pre-process
327 */
329 std::map<std::string, std::vector<VectorDouble>> &field_map,
330 std::map<std::string, std::vector<MatrixDouble>> &grad_map) {
333 }
334
335 protected:
336 inline static auto resizeAndSet(MatrixBoundedArray<TYPE, 9> &m) {
337 m.resize(3, 3, false);
339 };
340 };
341
344
345 std::vector<VectorDouble> &valuesAtGaussPts;
346 std::vector<MatrixDouble> &gradientAtGaussPts;
348
349 OpGetDataAtGaussPts(const std::string field_name,
350 std::vector<VectorDouble> &values_at_gauss_pts,
351 std::vector<MatrixDouble> &gradient_at_gauss_pts);
352
353 /** \brief operator calculating deformation gradient
354 *
355 * temperature gradient is calculated multiplying derivatives of shape
356 * functions by degrees of freedom
357 */
360 };
361
363 OpGetCommonDataAtGaussPts(const std::string field_name,
364 CommonData &common_data);
365 };
366
367 /**
368 * \brief Operator performs automatic differentiation.
369 */
372
373 BlockData &dAta; ///< Structure keeping data about problem, like material
374 ///< parameters
375 CommonData &commonData; ///< Structure keeping data abut this particular
376 ///< element, e.g. gradient of deformation at
377 ///< integration points
378 int tAg; //,lastId; ///< ADOL-C tag used for recording operations
379 int adlocReturnValue; ///< return value from ADOL-C, if non-zero that is
380 ///< error.
381 bool jAcobian; ///< if true Jacobian is calculated
382 bool fUnction; ///< if true stress i calculated
383 bool aLe; ///< true if arbitrary Lagrangian-Eulerian formulation
384 bool fieldDisp; ///< true if field of displacements is given, usually
385 ///< spatial positions are given.
386
387 /**
388 \brief Construct operator to calculate Piola-Kirchhoff stress or its
389 derivatives over gradient deformation
390
391 \param field_name approximation field name of spatial positions or
392 displacements \param data reference to block data (what is Young modulus,
393 Poisson ratio or what elements are part of the block) \param tag adol-c
394 unique tag of the tape \param jacobian if true derivative of Piola Stress
395 is calculated otherwise just stress is calculated \param field_disp if
396 true approximation field keeps displacements not spatial positions
397
398 */
400 BlockData &data, CommonData &common_data,
401 int tag, bool jacobian, bool ale,
402 bool field_disp);
403
406
407 std::vector<MatrixDouble> *ptrh;
408 std::vector<MatrixDouble> *ptrH;
409
413
414 /**
415 * \brief Calculate Paola-Kirchhoff I stress
416 * @return error code
417 */
418 virtual MoFEMErrorCode calculateStress(const int gg);
419
420 /**
421 * \brief Record ADOL-C tape
422 * @return error code
423 */
424 virtual MoFEMErrorCode recordTag(const int gg);
425
426 /**
427 * \brief Play ADOL-C tape
428 * @return error code
429 */
430 virtual MoFEMErrorCode playTag(const int gg);
431
432 /**
433 * \brief Cgeck if tape is recorded for given integration point
434 * @param gg integration point
435 * @return true if tag is recorded
436 */
437 virtual bool recordTagForIntegrationPoint(const int gg) {
438 // return true;
439 if (gg == 0)
440 return true;
441 return false;
442 }
443
444 /**
445 * \brief Calculate stress or jacobian at gauss points
446 *
447 * @param row_side
448 * @param row_type
449 * @param row_data
450 * @return error code
451 */
452 MoFEMErrorCode doWork(int row_side, EntityType row_type,
454 };
455
456 /**
457 * \brief Calculate explicit derivative of free energy
458 */
461
464
465 int tAg; ///< tape tag
466 bool gRadient; ///< if set true gradient of energy is calculated
467 bool hEssian; ///< if set true hessian of energy is calculated
468 bool aLe; ///< true if arbitrary Lagrangian-Eulerian formulation
469 bool fieldDisp; ///< true if displacements instead spatial positions used
470
471 OpJacobianEnergy(const std::string field_name, ///< field name for spatial
472 ///< positions or
473 ///< displacements
474 BlockData &data, CommonData &common_data, int tag,
475 bool gradient, bool hessian, bool ale, bool field_disp);
476
479
480 std::vector<MatrixDouble> *ptrh;
481 std::vector<MatrixDouble> *ptrH;
482
486
487 /**
488 * \brief Check if tape is recorded for given integration point
489 * @param gg integration point
490 * @return true if tag is recorded
491 */
492 virtual bool recordTagForIntegrationPoint(const int gg) {
493 if (gg == 0)
494 return true;
495 return false;
496 }
497
498 /**
499 * \brief Calculate Paola-Kirchhoff I stress
500 * @return error code
501 */
502 virtual MoFEMErrorCode calculateEnergy(const int gg);
503
504 /**
505 * \brief Record ADOL-C tape
506 * @return error code
507 */
508 virtual MoFEMErrorCode recordTag(const int gg);
509
510 /**
511 * \brief Play ADOL-C tape
512 * @return error code
513 */
514 virtual MoFEMErrorCode playTag(const int gg);
515
516 MoFEMErrorCode doWork(int row_side, EntityType row_type,
518 };
519
522
526 bool aLe;
527
528 ublas::vector<int> iNdices;
529 OpRhsPiolaKirchhoff(const std::string field_name, BlockData &data,
530 CommonData &common_data);
531
533 MoFEMErrorCode doWork(int row_side, EntityType row_type,
535
536 virtual MoFEMErrorCode aSemble(int row_side, EntityType row_type,
538 };
539
540 struct OpEnergy
542
545 SmartPetscObj<Vec> ghostVec;
547
548 OpEnergy(const std::string field_name, BlockData &data,
549 CommonData &common_data, SmartPetscObj<Vec> ghost_vec,
550 bool field_disp);
551
552 MoFEMErrorCode doWork(int row_side, EntityType row_type,
554 };
555
558
561 int tAg;
562 bool aLe;
563
564 ublas::vector<int> rowIndices;
565 ublas::vector<int> colIndices;
566
567 OpLhsPiolaKirchhoff_dx(const std::string vel_field,
568 const std::string field_name, BlockData &data,
569 CommonData &common_data);
570
572
573 /**
574 \brief Directive of Piola Kirchhoff stress over spatial DOFs
575
576 This project derivative \f$\frac{\partial P}{\partial F}\f$, that is
577 \f[
578 \frac{\partial P}{\partial x_\textrm{DOF}} = \frac{\partial P}{\partial
579 F}\frac{\partial F}{\partial x_\textrm{DOF}}, \f] where second therm
580 \f$\frac{\partial F}{\partial x_\textrm{DOF}}\f$ is derivative of shape
581 function
582
583 */
585 int gg);
586
587 virtual MoFEMErrorCode aSemble(int row_side, int col_side,
588 EntityType row_type, EntityType col_type,
591
592 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
593 EntityType col_type,
596 };
597
599
600 OpLhsPiolaKirchhoff_dX(const std::string vel_field,
601 const std::string field_name, BlockData &data,
602 CommonData &common_data);
603
604 /// \brief Derivative of Piola Kirchhoff stress over material DOFs
606
607 MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type,
608 EntityType col_type,
611 };
612
614
615 OpJacobianEshelbyStress(const std::string field_name, BlockData &data,
616 CommonData &common_data, int tag, bool jacobian,
617 bool ale);
618
619 MoFEMErrorCode calculateStress(const int gg);
620 };
621
623
624 OpRhsEshelbyStress(const std::string field_name, BlockData &data,
625 CommonData &common_data);
626 };
627
628 /**
629 * \deprecated name with spelling mistake
630 */
632
634
635 OpLhsEshelby_dx(const std::string vel_field, const std::string field_name,
636 BlockData &data, CommonData &common_data);
637
639 };
640
642
643 OpLhsEshelby_dX(const std::string vel_field, const std::string field_name,
644 BlockData &data, CommonData &common_data);
645
647 };
648
651 materialDoublePtr,
653 materialAdoublePtr);
654
656 const std::string element_name,
657 const std::string spatial_position_field_name,
658 const std::string material_position_field_name = "MESH_NODE_POSITIONS",
659 const bool ale = false);
660
661 /** \brief Set operators to calculate left hand tangent matrix and right hand
662 * residual
663 *
664 * \param fun class needed to calculate Piola Kirchhoff I Stress tensor
665 * \param spatial_position_field_name name of approximation field
666 * \param material_position_field_name name of field to define geometry
667 * \param ale true if arbitrary Lagrangian Eulerian formulation
668 * \param field_disp true if approximation field represents displacements
669 * otherwise it is field of spatial positions
670 */
672 const std::string spatial_position_field_name,
673 const std::string material_position_field_name = "MESH_NODE_POSITIONS",
674 const bool ale = false, const bool field_disp = false);
675};
676
677#endif //__NONLINEAR_ELASTIC_HPP
678
679/**
680 * \defgroup nonlinear_elastic_elem NonLinear Elastic Element
681 * \ingroup user_modules
682 * \defgroup user_modules User modules
683 **/
EntitiesFieldData::EntData EntData
Kronecker Delta class.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define DEPRECATED
Definition: definitions.h:17
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
#define MU(E, NU)
Definition: fem_tools.h:23
#define LAMBDA(E, NU)
Definition: fem_tools.h:22
FTensor::Index< 'm', SPACE_DIM > m
constexpr auto t_kd
@ TYPE
Definition: inflate.h:32
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr)
Definition: Templates.hpp:978
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1199
double trace(FTensor::Tensor2_symmetric< T, SPACE_DIM > &t_stress)
Definition: PlasticOps.hpp:337
constexpr auto field_name
Deprecated interface functions.
FTensor::Tensor2< double *, 3, 3 > & getJac()
get element Jacobian
data for calculation heat conductivity and heat capacity elements
boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< adouble > > materialAdoublePtr
boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< double > > materialDoublePtr
Range tEts
constrains elements in block set
common data used by volume elements
std::vector< MatrixDouble > jacStress
this is simply material tangent operator
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
std::vector< MatrixDouble3by3 > sTress
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
Implementation of elastic (non-linear) St. Kirchhoff equation.
FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 > t_h
virtual MoFEMErrorCode calculatesIGma_EshelbyStress(const BlockData block_data, boost::shared_ptr< const NumeredEntFiniteElement > fe_ptr)
Calculate Eshelby stress.
FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 > t_invH
FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 > t_C
virtual MoFEMErrorCode calculateP_PiolaKirchhoffI(const BlockData block_data, boost::shared_ptr< const NumeredEntFiniteElement > fe_ptr)
Function overload to implement user material.
FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 > t_sIGma
FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 > t_F
FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 > t_sigmaCauchy
FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 > t_P
virtual MoFEMErrorCode calculateElasticEnergy(const BlockData block_data, boost::shared_ptr< const NumeredEntFiniteElement > fe_ptr)
Calculate elastic energy density.
FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 > t_H
FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 > t_E
FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 > t_invF
FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 > t_S
virtual MoFEMErrorCode setUserActiveVariables(int &nb_active_variables)
add additional active variables
virtual MoFEMErrorCode getDataOnPostProcessor(std::map< std::string, std::vector< VectorDouble > > &field_map, std::map< std::string, std::vector< MatrixDouble > > &grad_map)
Do operations when pre-process.
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator * opPtr
pointer to finite element tetrahedral operator
static auto resizeAndSet(MatrixBoundedArray< TYPE, 9 > &m)
virtual MoFEMErrorCode calculateCauchyStress(const BlockData block_data, boost::shared_ptr< const NumeredEntFiniteElement > fe_ptr)
Function overload to implement user material.
virtual MoFEMErrorCode setUserActiveVariables(VectorDouble &activeVariables)
Add additional independent variables More complex physical models depend on gradient of defamation an...
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode postProcess()
function is run at the end of loop
int getRule(int order)
it is used to calculate nb. of Gauss integration points
OpEnergy(const std::string field_name, BlockData &data, CommonData &common_data, SmartPetscObj< Vec > ghost_vec, bool field_disp)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpGetCommonDataAtGaussPts(const std::string field_name, CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
operator calculating deformation gradient
OpGetDataAtGaussPts(const std::string field_name, std::vector< VectorDouble > &values_at_gauss_pts, std::vector< MatrixDouble > &gradient_at_gauss_pts)
Calculate explicit derivative of free energy.
bool gRadient
if set true gradient of energy is calculated
bool fieldDisp
true if displacements instead spatial positions used
virtual bool recordTagForIntegrationPoint(const int gg)
Check if tape is recorded for given integration point.
virtual MoFEMErrorCode playTag(const int gg)
Play ADOL-C tape.
bool hEssian
if set true hessian of energy is calculated
OpJacobianEnergy(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool gradient, bool hessian, bool ale, bool field_disp)
virtual MoFEMErrorCode recordTag(const int gg)
Record ADOL-C tape.
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
virtual MoFEMErrorCode calculateEnergy(const int gg)
Calculate Paola-Kirchhoff I stress.
bool aLe
true if arbitrary Lagrangian-Eulerian formulation
MoFEMErrorCode calculateStress(const int gg)
Calculate Paola-Kirchhoff I stress.
OpJacobianEshelbyStress(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian, bool ale)
virtual MoFEMErrorCode recordTag(const int gg)
Record ADOL-C tape.
bool aLe
true if arbitrary Lagrangian-Eulerian formulation
virtual MoFEMErrorCode playTag(const int gg)
Play ADOL-C tape.
OpJacobianPiolaKirchhoffStress(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian, bool ale, bool field_disp)
Construct operator to calculate Piola-Kirchhoff stress or its derivatives over gradient deformation.
virtual bool recordTagForIntegrationPoint(const int gg)
Cgeck if tape is recorded for given integration point.
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Calculate stress or jacobian at gauss points.
virtual MoFEMErrorCode calculateStress(const int gg)
Calculate Paola-Kirchhoff I stress.
OpLhsEshelby_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
OpLhsEshelby_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
OpLhsPiolaKirchhoff_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
virtual MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpLhsPiolaKirchhoff_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
OpRhsEshelbyStress(const std::string field_name, BlockData &data, CommonData &common_data)
virtual MoFEMErrorCode aSemble(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpRhsPiolaKirchhoff(const std::string field_name, BlockData &data, CommonData &common_data)
structure grouping operators and data used for calculation of nonlinear elastic element
MyVolumeFE & getLoopFeLhs()
get lhs volume element
MyVolumeFE feRhs
calculate right hand side for tetrahedral elements
NonlinearElasticElement(MoFEM::Interface &m_field, short int tag)
MoFEMErrorCode addElement(const std::string element_name, const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false)
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
FTensor::Index< 'k', 3 > k
DEPRECATED typedef OpRhsEshelbyStress OpRhsEshelbyStrees
FTensor::Index< 'j', 3 > j
MoFEMErrorCode setBlocks(boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< double > > materialDoublePtr, boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< adouble > > materialAdoublePtr)
MyVolumeFE & getLoopFeRhs()
get rhs volume element
MyVolumeFE & getLoopFeEnergy()
get energy fe
virtual ~NonlinearElasticElement()=default
MyVolumeFE feEnergy
calculate elastic energy
FTensor::Index< 'i', 3 > i
MoFEMErrorCode setOperators(const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false, const bool field_disp=false)
Set operators to calculate left hand tangent matrix and right hand residual.