v0.10.0
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  * This file is part of MoFEM.
10  * MoFEM is free software: you can redistribute it and/or modify it under
11  * the terms of the GNU Lesser General Public License as published by the
12  * Free Software Foundation, either version 3 of the License, or (at your
13  * option) any later version.
14  *
15  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
16  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18  * License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
22 
23 #ifndef __NONLINEAR_ELASTIC_HPP
24 #define __NONLINEAR_ELASTIC_HPP
25 
26 #ifndef WITH_ADOL_C
27 #error "MoFEM need to be compiled with ADOL-C"
28 #endif
29 
30 /** \brief structure grouping operators and data used for calculation of
31  * nonlinear elastic element \ingroup nonlinear_elastic_elem
32  *
33  * In order to assemble matrices and right hand vectors, the loops over
34  * elements, entities over that elements and finally loop over integration
35  * points are executed.
36  *
37  * Following implementation separate those three categories of loops and to each
38  * loop attach operator.
39  *
40  */
42 
43  /// \brief definition of volume element
45 
46  Mat A;
47  Vec F;
48 
49  int addToRule;
50 
51  MyVolumeFE(MoFEM::Interface &m_field);
52  virtual ~MyVolumeFE();
53 
54  /** \brief it is used to calculate nb. of Gauss integration points
55  *
56  * for more details pleas look
57  * Reference:
58  *
59  * Albert Nijenhuis, Herbert Wilf,
60  * Combinatorial Algorithms for Computers and Calculators,
61  * Second Edition,
62  * Academic Press, 1978,
63  * ISBN: 0-12-519260-6,
64  * LC: QA164.N54.
65  *
66  * More details about algorithm
67  * http://people.sc.fsu.edu/~jburkardt/cpp_src/gm_rule/gm_rule.html
68  **/
69  int getRule(int order);
70 
71  Vec V;
72  double eNergy;
73 
76  };
77 
78  MyVolumeFE feRhs; ///< calculate right hand side for tetrahedral elements
79  MyVolumeFE &getLoopFeRhs() { return feRhs; } ///< get rhs volume element
80  MyVolumeFE feLhs; //< calculate left hand side for tetrahedral elements
81  MyVolumeFE &getLoopFeLhs() { return feLhs; } ///< get lhs volume element
82 
83  MyVolumeFE feEnergy; ///< calculate elastic energy
84  MyVolumeFE &getLoopFeEnergy() { return feEnergy; } ///< get energy fe
85 
87  short int tAg;
88 
89  NonlinearElasticElement(MoFEM::Interface &m_field, short int tag);
90  virtual ~NonlinearElasticElement() = default;
91 
92  template <typename TYPE> struct FunctionsToCalculatePiolaKirchhoffI;
93 
94  /** \brief data for calculation het conductivity and heat capacity elements
95  * \ingroup nonlinear_elastic_elem
96  */
97  struct BlockData {
98  int iD;
99  double E;
100  double PoissonRatio;
101  // Eberlein Fibres stiffness properties
102  double k1, k2;
103  Range tEts; ///< constrains elements in block set
104  boost::shared_ptr<FunctionsToCalculatePiolaKirchhoffI<adouble>>
106  boost::shared_ptr<FunctionsToCalculatePiolaKirchhoffI<double>>
110  };
111 
112  std::map<int, BlockData>
113  setOfBlocks; ///< maps block set id with appropriate BlockData
114 
115  /** \brief common data used by volume elements
116  * \ingroup nonlinear_elastic_elem
117  */
118  struct CommonData {
119 
120  std::map<std::string, std::vector<VectorDouble>> dataAtGaussPts;
121  std::map<std::string, std::vector<MatrixDouble>> gradAtGaussPts;
124  std::vector<MatrixDouble3by3> sTress;
125  std::vector<MatrixDouble>
126  jacStress; ///< this is simply material tangent operator
127 
128  // This part can be used to calculate stress directly from potential
129 
130  std::vector<double> eNergy;
131  std::vector<VectorDouble> jacEnergy;
132  std::vector<VectorDouble> hessianEnergy;
133  };
135 
136  /** \brief Implementation of elastic (non-linear) St. Kirchhoff equation
137  * \ingroup nonlinear_elastic_elem
138  */
139  template <typename TYPE> struct FunctionsToCalculatePiolaKirchhoffI {
140 
142 
143  /** \brief Calculate determinant of 3x3 matrix
144  */
146  ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>> &a,
147  TYPE &det) {
149  // a11a22a33
150  //+a21a32a13
151  //+a31a12a23
152  //-a11a32a23
153  //-a31a22a13
154  //-a21a12a33
155  // http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html
156  // http://mathworld.wolfram.com/MatrixInverse.html
157  det = a(0, 0) * a(1, 1) * a(2, 2) + a(1, 0) * a(2, 1) * a(0, 2) +
158  a(2, 0) * a(0, 1) * a(1, 2) - a(0, 0) * a(2, 1) * a(1, 2) -
159  a(2, 0) * a(1, 1) * a(0, 2) - a(1, 0) * a(0, 1) * a(2, 2);
161  }
162 
163  /** \brief Calculate inverse of 3x3 matrix
164  */
166  TYPE det,
167  ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>> &a,
168  ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>>
169  &inv_a) {
171  //
172  inv_a.resize(3, 3);
173  // http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html
174  // http://mathworld.wolfram.com/MatrixInverse.html
175  inv_a(0, 0) = a(1, 1) * a(2, 2) - a(1, 2) * a(2, 1);
176  inv_a(0, 1) = a(0, 2) * a(2, 1) - a(0, 1) * a(2, 2);
177  inv_a(0, 2) = a(0, 1) * a(1, 2) - a(0, 2) * a(1, 1);
178  inv_a(1, 0) = a(1, 2) * a(2, 0) - a(1, 0) * a(2, 2);
179  inv_a(1, 1) = a(0, 0) * a(2, 2) - a(0, 2) * a(2, 0);
180  inv_a(1, 2) = a(0, 2) * a(1, 0) - a(0, 0) * a(1, 2);
181  inv_a(2, 0) = a(1, 0) * a(2, 1) - a(1, 1) * a(2, 0);
182  inv_a(2, 1) = a(0, 1) * a(2, 0) - a(0, 0) * a(2, 1);
183  inv_a(2, 2) = a(0, 0) * a(1, 1) - a(0, 1) * a(1, 0);
184  inv_a /= det;
186  }
187 
188  double lambda, mu;
189  ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>> F, C,
190  E, S, invF, P, SiGma, h, H, invH;
191  TYPE J, eNergy, detH;
192 
193  int gG; ///< Gauss point number
194  CommonData *commonDataPtr; ///< common data shared between entities (f.e.
195  ///< field values at Gauss pts.)
197  *opPtr; ///< pointer to finite element tetrahedral operator
198 
201  C.resize(3, 3);
202  noalias(C) = prod(trans(F), F);
204  }
205 
208  E.resize(3, 3);
209  noalias(E) = C;
210  for (int dd = 0; dd < 3; dd++) {
211  E(dd, dd) -= 1;
212  }
213  E *= 0.5;
215  }
216 
217  // St. Venant–Kirchhoff Material
220  TYPE trE = 0;
221  for (int dd = 0; dd < 3; dd++) {
222  trE += E(dd, dd);
223  }
224  S.resize(3, 3);
225  S.clear();
226  for (int dd = 0; dd < 3; dd++) {
227  S(dd, dd) = trE * lambda;
228  }
229  S += 2 * mu * E;
231  }
232 
233  /** \brief Function overload to implement user material
234  *
235 
236  * Calculation of Piola Kirchhoff I is implemented by user. Tangent matrix
237  * user implemented physical equation is calculated using automatic
238  * differentiation.
239 
240  * \f$\mathbf{S} =
241  \lambda\textrm{tr}[\mathbf{E}]\mathbf{I}+2\mu\mathbf{E}\f$
242 
243  * Notes: <br>
244  * Number of actual Gauss point is accessed from variable gG. <br>
245  * Access to operator data structures is available by variable opPtr. <br>
246  * Access to common data is by commonDataPtr. <br>
247 
248  * \param block_data used to give access to material parameters
249  * \param fe_ptr pointer to element data structures
250 
251  For details look to: <br>
252  NONLINEAR CONTINUUM MECHANICS FOR FINITE ELEMENT ANALYSIS, Javier Bonet,
253  Richard D. Wood
254 
255  */
257  const BlockData block_data,
258  boost::shared_ptr<const NumeredEntFiniteElement> fe_ptr) {
260  lambda = LAMBDA(block_data.E, block_data.PoissonRatio);
261  mu = MU(block_data.E, block_data.PoissonRatio);
265  P.resize(3, 3);
266  noalias(P) = prod(F, S);
268  }
269 
270  /**
271  * \brief add additional active variables
272  *
273  * \note This member function if used should be implement by template member
274  * function Specialization, different implementation needed for TYPE=double
275  * or TYPE=adouble
276  *
277  * More complex physical models depend on gradient of defamation and some
278  * additional variables. For example can depend on temperature. This
279  * function adds additional independent variables to the model.
280  *
281  * @param nb_active_variables number of active variables
282  * @return error code
283  */
284  virtual MoFEMErrorCode setUserActiveVariables(int &nb_active_variables) {
287  }
288 
289  /**
290  * \brief Add additional independent variables
291  * More complex physical models depend on gradient of defamation and some
292  * additional variables. For example can depend on temperature. This
293  * function adds additional independent variables to the model.
294  *
295  * /note First 9 elements are reserved for gradient of deformation.
296  * @param activeVariables vector of deepened variables, values after index
297  * 9 should be add.
298  *
299  * @return error code
300  */
301  virtual MoFEMErrorCode
305  }
306 
307  /** \brief Calculate elastic energy density
308  *
309  * \f[\Psi =
310  * \frac{1}{2}\lambda(\textrm{tr}[\mathbf{E}])^2+\mu\mathbf{E}:\mathbf{E}\f]
311  */
313  const BlockData block_data,
314  boost::shared_ptr<const NumeredEntFiniteElement> fe_ptr) {
316  lambda = LAMBDA(block_data.E, block_data.PoissonRatio);
317  mu = MU(block_data.E, block_data.PoissonRatio);
320  TYPE trace = 0;
321  eNergy = 0;
322  for (int ii = 0; ii < 3; ii++) {
323  trace += E(ii, ii);
324  for (int jj = 0; jj < 3; jj++) {
325  TYPE e = E(ii, jj);
326  eNergy += mu * e * e;
327  }
328  }
329  eNergy += 0.5 * lambda * trace * trace;
331  }
332 
333  /** \brief Calculate Eshelby stress
334  */
336  const BlockData block_data,
337  boost::shared_ptr<const NumeredEntFiniteElement> fe_ptr) {
339  CHKERR calculateP_PiolaKirchhoffI(block_data, fe_ptr);
340  CHKERR calculateElasticEnergy(block_data, fe_ptr);
341  SiGma.resize(3, 3, false);
342  noalias(SiGma) = -prod(trans(F), P);
343  for (int dd = 0; dd < 3; dd++) {
344  SiGma(dd, dd) += eNergy;
345  }
347  }
348 
349  /** \brief Do operations when pre-process
350  */
352  std::map<std::string, std::vector<VectorDouble>> &field_map,
353  std::map<std::string, std::vector<MatrixDouble>> &grad_map) {
356  }
357  };
358 
361 
362  std::vector<VectorDouble> &valuesAtGaussPts;
363  std::vector<MatrixDouble> &gradientAtGaussPts;
364  const EntityType zeroAtType;
365 
366  OpGetDataAtGaussPts(const std::string field_name,
367  std::vector<VectorDouble> &values_at_gauss_pts,
368  std::vector<MatrixDouble> &gradient_at_gauss_pts);
369 
370  /** \brief operator calculating deformation gradient
371  *
372  * temperature gradient is calculated multiplying derivatives of shape
373  * functions by degrees of freedom
374  */
375  MoFEMErrorCode doWork(int side, EntityType type,
377  };
378 
380  OpGetCommonDataAtGaussPts(const std::string field_name,
381  CommonData &common_data);
382  };
383 
384  /**
385  * \brief Operator performs automatic differentiation.
386  */
389 
390  BlockData &dAta; ///< Structure keeping data about problem, like material
391  ///< parameters
392  CommonData &commonData; ///< Structure keeping data abut this particular
393  ///< element, e.g. gradient of deformation at
394  ///< integration points
395  int tAg; //,lastId; ///< ADOL-C tag used for recording operations
396  int adlocReturnValue; ///< return value from ADOL-C, if non-zero that is
397  ///< error.
398  bool jAcobian; ///< if true Jacobian is calculated
399  bool fUnction; ///< if true stress i calculated
400  bool aLe; ///< true if arbitrary Lagrangian-Eulerian formulation
401  bool fieldDisp; ///< true if field of displacements is given, usually
402  ///< spatial positions are given.
403 
404  /**
405  \brief Construct operator to calculate Piola-Kirchhoff stress or its
406  derivatives over gradient deformation
407 
408  \param field_name approximation field name of spatial positions or
409  displacements \param data reference to block data (what is Young modulus,
410  Poisson ratio or what elements are part of the block) \param tag adol-c
411  unique tag of the tape \param jacobian if true derivative of Piola Stress
412  is calculated otherwise just stress is calculated \param field_disp if
413  true approximation field keeps displacements not spatial positions
414 
415  */
416  OpJacobianPiolaKirchhoffStress(const std::string field_name,
417  BlockData &data, CommonData &common_data,
418  int tag, bool jacobian, bool ale,
419  bool field_disp);
420 
423 
424  std::vector<MatrixDouble> *ptrh;
425  std::vector<MatrixDouble> *ptrH;
426 
427  /**
428  * \brief Calculate Paola-Kirchhoff I stress
429  * @return error code
430  */
431  virtual MoFEMErrorCode calculateStress(const int gg);
432 
433  /**
434  * \brief Record ADOL-C tape
435  * @return error code
436  */
437  virtual MoFEMErrorCode recordTag(const int gg);
438 
439  /**
440  * \brief Play ADOL-C tape
441  * @return error code
442  */
443  virtual MoFEMErrorCode playTag(const int gg);
444 
445  /**
446  * \brief Cgeck if tape is recorded for given integration point
447  * @param gg integration point
448  * @return true if tag is recorded
449  */
450  virtual bool recordTagForIntegrationPoint(const int gg) {
451  // return true;
452  if (gg == 0)
453  return true;
454  return false;
455  }
456 
457  /**
458  * \brief Calculate stress or jacobian at gauss points
459  *
460  * @param row_side
461  * @param row_type
462  * @param row_data
463  * @return error code
464  */
465  MoFEMErrorCode doWork(int row_side, EntityType row_type,
467  };
468 
469  /**
470  * \brief Calculate explicit derivative of free energy
471  */
474 
477 
478  int tAg; ///< tape tag
479  bool gRadient; ///< if set true gradient of energy is calculated
480  bool hEssian; ///< if set true hessian of energy is calculated
481  bool aLe; ///< true if arbitrary Lagrangian-Eulerian formulation
482  bool fieldDisp; ///< true if displacements instead spatial positions used
483 
484  OpJacobianEnergy(const std::string field_name, ///< field name for spatial
485  ///< positions or
486  ///< displacements
487  BlockData &data, CommonData &common_data, int tag,
488  bool gradient, bool hessian, bool ale, bool field_disp);
489 
492 
493  std::vector<MatrixDouble> *ptrh;
494  std::vector<MatrixDouble> *ptrH;
495 
496  /**
497  * \brief Check if tape is recorded for given integration point
498  * @param gg integration point
499  * @return true if tag is recorded
500  */
501  virtual bool recordTagForIntegrationPoint(const int gg) {
502  if (gg == 0)
503  return true;
504  return false;
505  }
506 
507  /**
508  * \brief Calculate Paola-Kirchhoff I stress
509  * @return error code
510  */
511  virtual MoFEMErrorCode calculateEnergy(const int gg);
512 
513  /**
514  * \brief Record ADOL-C tape
515  * @return error code
516  */
517  virtual MoFEMErrorCode recordTag(const int gg);
518 
519  /**
520  * \brief Play ADOL-C tape
521  * @return error code
522  */
523  virtual MoFEMErrorCode playTag(const int gg);
524 
525  MoFEMErrorCode doWork(int row_side, EntityType row_type,
527  };
528 
531 
534  bool fieldDisp;
535  bool aLe;
536 
537  ublas::vector<int> iNdices;
538  OpRhsPiolaKirchhoff(const std::string field_name, BlockData &data,
539  CommonData &common_data);
540 
542  MoFEMErrorCode doWork(int row_side, EntityType row_type,
544 
545  virtual MoFEMErrorCode aSemble(int row_side, EntityType row_type,
547  };
548 
549  struct OpEnergy
551 
554  Vec ghostVec;
555  bool fieldDisp;
556 
557  OpEnergy(const std::string field_name, BlockData &data,
558  CommonData &common_data, Vec ghost_vec, bool field_disp);
559  ~OpEnergy();
560 
561  MoFEMErrorCode doWork(int row_side, EntityType row_type,
563  };
564 
567 
570  int tAg;
571  bool aLe;
572 
573  ublas::vector<int> rowIndices;
574  ublas::vector<int> colIndices;
575 
576  OpLhsPiolaKirchhoff_dx(const std::string vel_field,
577  const std::string field_name, BlockData &data,
578  CommonData &common_data);
579 
581 
582  /**
583  \brief Directive of Piola Kirchhoff stress over spatial DOFs
584 
585  This project derivative \f$\frac{\partial P}{\partial F}\f$, that is
586  \f[
587  \frac{\partial P}{\partial x_\textrm{DOF}} = \frac{\partial P}{\partial
588  F}\frac{\partial F}{\partial x_\textrm{DOF}}, \f] where second therm
589  \f$\frac{\partial F}{\partial x_\textrm{DOF}}\f$ is derivative of shape
590  function
591 
592  */
594  int gg);
595 
596  virtual MoFEMErrorCode aSemble(int row_side, int col_side,
597  EntityType row_type, EntityType col_type,
600 
601  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
602  EntityType col_type,
605  };
606 
608 
609  OpLhsPiolaKirchhoff_dX(const std::string vel_field,
610  const std::string field_name, BlockData &data,
611  CommonData &common_data);
612 
613  /// \brief Derivative of Piola Kirchhoff stress over material DOFs
615 
616  MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type,
617  EntityType col_type,
620  };
621 
623 
624  OpJacobianEshelbyStress(const std::string field_name, BlockData &data,
625  CommonData &common_data, int tag, bool jacobian,
626  bool ale);
627 
628  MoFEMErrorCode calculateStress(const int gg);
629  };
630 
632 
633  OpRhsEshelbyStress(const std::string field_name, BlockData &data,
634  CommonData &common_data);
635  };
636 
637  /**
638  * \deprecated name with spelling mistake
639  */
641 
643 
644  OpLhsEshelby_dx(const std::string vel_field, const std::string field_name,
645  BlockData &data, CommonData &common_data);
646 
648  };
649 
651 
652  OpLhsEshelby_dX(const std::string vel_field, const std::string field_name,
653  BlockData &data, CommonData &common_data);
654 
656  };
657 
660  materialDoublePtr,
662  materialAdoublePtr);
663 
665  const std::string element_name,
666  const std::string spatial_position_field_name,
667  const std::string material_position_field_name = "MESH_NODE_POSITIONS",
668  const bool ale = false);
669 
670  /** \brief Set operators to calculate left hand tangent matrix and right hand
671  * residual
672  *
673  * \param fun class needed to calculate Piola Kirchhoff I Stress tensor
674  * \param spatial_position_field_name name of approximation field
675  * \param material_position_field_name name of field to define geometry
676  * \param ale true if arbitrary Lagrangian Eulerian formulation
677  * \param field_disp true if approximation field represents displacements
678  * otherwise it is field of spatial positions
679  */
681  const std::string spatial_position_field_name,
682  const std::string material_position_field_name = "MESH_NODE_POSITIONS",
683  const bool ale = false, const bool field_disp = false);
684 };
685 
686 #endif //__NONLINEAR_ELASTIC_HPP
687 
688 /**
689  * \defgroup nonlinear_elastic_elem NonLinear Elastic Element
690  * \ingroup user_modules
691  * \defgroup user_modules User modules
692  **/
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.
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > H
MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
#define LAMBDA(E, NU)
Definition: fem_tools.h:32
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > S
Deprecated interface functions.
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.
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > F
boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< double > > materialDoublePtr
boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< adouble > > materialAdoublePtr
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.
OpJacobianEnergy(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool gradient, bool hessian, bool ale, bool field_disp)
FTensor::Tensor2< double *, 3, 3 > & getJac()
get element Jacobian
std::vector< MatrixDouble3by3 > sTress
OpGetCommonDataAtGaussPts(const std::string field_name, CommonData &common_data)
bool hEssian
if set true hessian of energy is calculated
MoFEMErrorCode dEterminant(ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 >> &a, TYPE &det)
Calculate determinant of 3x3 matrix.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
Implementation of elastic (non-linear) St. Kirchhoff equation.
OpLhsEshelby_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
virtual MoFEMErrorCode calculateStress(const int gg)
Calculate Paola-Kirchhoff I stress.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
virtual bool recordTagForIntegrationPoint(const int gg)
Check if tape is recorded for given integration point.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
MoFEMErrorCode calculateStress(const int gg)
Calculate Paola-Kirchhoff I stress.
virtual MoFEMErrorCode calculateSiGma_EshelbyStress(const BlockData block_data, boost::shared_ptr< const NumeredEntFiniteElement > fe_ptr)
Calculate Eshelby stress.
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > P
Calculate explicit derivative of free energy.
OpJacobianEshelbyStress(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian, bool ale)
MyVolumeFE feEnergy
calculate elastic energy
virtual MoFEMErrorCode calculateEnergy(const int gg)
Calculate Paola-Kirchhoff I stress.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
MyVolumeFE feRhs
calculate right hand side for tetrahedral elements
MyVolumeFE & getLoopFeRhs()
get rhs volume element
MoFEMErrorCode postProcess()
function is run at the end of loop
MyVolumeFE & getLoopFeLhs()
get lhs volume element
virtual MoFEMErrorCode playTag(const int gg)
Play ADOL-C tape.
MoFEMErrorCode setBlocks(boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< double >> materialDoublePtr, boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< adouble >> materialAdoublePtr)
bool gRadient
if set true gradient of energy is calculated
Range tEts
constrains elements in block set
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
OpLhsEshelby_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
MyVolumeFE & getLoopFeEnergy()
get energy fe
OpLhsPiolaKirchhoff_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
OpLhsPiolaKirchhoff_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
common data used by volume elements
NonlinearElasticElement(MoFEM::Interface &m_field, short int tag)
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Calculate stress or jacobian at gauss points.
OpRhsPiolaKirchhoff(const std::string field_name, BlockData &data, CommonData &common_data)
OpEnergy(const std::string field_name, BlockData &data, CommonData &common_data, Vec ghost_vec, bool field_disp)
DataForcesAndSourcesCore::EntData EntData
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > SiGma
MoFEMErrorCode preProcess()
function is run at the beginning of loop
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
virtual ~NonlinearElasticElement()=default
virtual bool recordTagForIntegrationPoint(const int gg)
Cgeck if tape is recorded for given integration point.
DEPRECATED typedef OpRhsEshelbyStress OpRhsEshelbyStrees
virtual MoFEMErrorCode calculateP_PiolaKirchhoffI(const BlockData block_data, boost::shared_ptr< const NumeredEntFiniteElement > fe_ptr)
Function overload to implement user material.
double trace(FTensor::Tensor2_symmetric< T, 2 > &t_stress)
Definition: PlasticOps.hpp:240
std::vector< MatrixDouble > jacStress
this is simply material tangent operator
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator * opPtr
pointer to finite element tetrahedral operator
virtual MoFEMErrorCode setUserActiveVariables(int &nb_active_variables)
add additional active variables
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
int getRule(int order)
it is used to calculate nb. of Gauss integration points
#define CHKERR
Inline error check.
Definition: definitions.h:604
#define MU(E, NU)
Definition: fem_tools.h:33
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > h
constexpr int order
virtual MoFEMErrorCode calculateElasticEnergy(const BlockData block_data, boost::shared_ptr< const NumeredEntFiniteElement > fe_ptr)
Calculate elastic energy density.
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
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)
#define DEPRECATED
Definition: definitions.h:29
virtual MoFEMErrorCode recordTag(const int gg)
Record ADOL-C tape.
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > C
bool aLe
true if arbitrary Lagrangian-Eulerian formulation
MoFEMErrorCode iNvert(TYPE det, ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 >> &a, ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 >> &inv_a)
Calculate inverse of 3x3 matrix.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
operator calculating deformation gradient
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > invH
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
bool aLe
true if arbitrary Lagrangian-Eulerian formulation
virtual MoFEMErrorCode playTag(const int gg)
Play ADOL-C tape.
virtual MoFEMErrorCode recordTag(const int gg)
Record ADOL-C tape.
OpRhsEshelbyStress(const std::string field_name, BlockData &data, CommonData &common_data)
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
OpGetDataAtGaussPts(const std::string field_name, std::vector< VectorDouble > &values_at_gauss_pts, std::vector< MatrixDouble > &gradient_at_gauss_pts)
data for calculation heat conductivity and heat capacity elements
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > E
virtual MoFEMErrorCode setUserActiveVariables(VectorDouble &activeVariables)
Add additional independent variables More complex physical models depend on gradient of defamation an...
bool fieldDisp
true if displacements instead spatial positions used
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > invF
structure grouping operators and data used for calculation of nonlinear elastic elementIn order to as...
virtual MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
virtual MoFEMErrorCode aSemble(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)