v0.9.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  * 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  ~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 
91  template <typename TYPE> struct FunctionsToCalculatePiolaKirchhoffI;
92 
93  /** \brief data for calculation het conductivity and heat capacity elements
94  * \ingroup nonlinear_elastic_elem
95  */
96  struct BlockData {
97  int iD;
98  double E;
99  double PoissonRatio;
100  // Eberlein Fibres stiffness properties
101  double k1, k2;
102  Range tEts; ///< constrains elements in block set
103  boost::shared_ptr<FunctionsToCalculatePiolaKirchhoffI<adouble>>
105  boost::shared_ptr<FunctionsToCalculatePiolaKirchhoffI<double>>
109  };
110 
111  std::map<int, BlockData>
112  setOfBlocks; ///< maps block set id with appropriate BlockData
113 
114  /** \brief common data used by volume elements
115  * \ingroup nonlinear_elastic_elem
116  */
117  struct CommonData {
118 
119  std::map<std::string, std::vector<VectorDouble>> dataAtGaussPts;
120  std::map<std::string, std::vector<MatrixDouble>> gradAtGaussPts;
123  std::vector<MatrixDouble3by3> sTress;
124  std::vector<MatrixDouble>
125  jacStress; ///< this is simply material tangent operator
126 
127  // This part can be used to calculate stress directly from potential
128 
129  std::vector<double> eNergy;
130  std::vector<VectorDouble> jacEnergy;
131  std::vector<VectorDouble> hessianEnergy;
132  };
134 
135  /** \brief Implementation of elastic (non-linear) St. Kirchhoff equation
136  * \ingroup nonlinear_elastic_elem
137  */
138  template <typename TYPE> struct FunctionsToCalculatePiolaKirchhoffI {
139 
141 
142  /** \brief Calculate determinant of 3x3 matrix
143  */
145  ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>> &a,
146  TYPE &det) {
148  // a11a22a33
149  //+a21a32a13
150  //+a31a12a23
151  //-a11a32a23
152  //-a31a22a13
153  //-a21a12a33
154  // http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html
155  // http://mathworld.wolfram.com/MatrixInverse.html
156  det = a(0, 0) * a(1, 1) * a(2, 2) + a(1, 0) * a(2, 1) * a(0, 2) +
157  a(2, 0) * a(0, 1) * a(1, 2) - a(0, 0) * a(2, 1) * a(1, 2) -
158  a(2, 0) * a(1, 1) * a(0, 2) - a(1, 0) * a(0, 1) * a(2, 2);
160  }
161 
162  /** \brief Calculate inverse of 3x3 matrix
163  */
165  TYPE det,
166  ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>> &a,
167  ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>>
168  &inv_a) {
170  //
171  inv_a.resize(3, 3);
172  // http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html
173  // http://mathworld.wolfram.com/MatrixInverse.html
174  inv_a(0, 0) = a(1, 1) * a(2, 2) - a(1, 2) * a(2, 1);
175  inv_a(0, 1) = a(0, 2) * a(2, 1) - a(0, 1) * a(2, 2);
176  inv_a(0, 2) = a(0, 1) * a(1, 2) - a(0, 2) * a(1, 1);
177  inv_a(1, 0) = a(1, 2) * a(2, 0) - a(1, 0) * a(2, 2);
178  inv_a(1, 1) = a(0, 0) * a(2, 2) - a(0, 2) * a(2, 0);
179  inv_a(1, 2) = a(0, 2) * a(1, 0) - a(0, 0) * a(1, 2);
180  inv_a(2, 0) = a(1, 0) * a(2, 1) - a(1, 1) * a(2, 0);
181  inv_a(2, 1) = a(0, 1) * a(2, 0) - a(0, 0) * a(2, 1);
182  inv_a(2, 2) = a(0, 0) * a(1, 1) - a(0, 1) * a(1, 0);
183  inv_a /= det;
185  }
186 
187  double lambda, mu;
188  ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>> F, C,
189  E, S, invF, P, SiGma, h, H, invH;
190  TYPE J, eNergy, detH;
191 
192  int gG; ///< Gauss point number
193  CommonData *commonDataPtr; ///< common data shared between entities (f.e.
194  ///< field values at Gauss pts.)
196  *opPtr; ///< pointer to finite element tetrahedral operator
197 
200  C.resize(3, 3);
201  noalias(C) = prod(trans(F), F);
203  }
204 
207  E.resize(3, 3);
208  noalias(E) = C;
209  for (int dd = 0; dd < 3; dd++) {
210  E(dd, dd) -= 1;
211  }
212  E *= 0.5;
214  }
215 
216  // St. Venant–Kirchhoff Material
219  TYPE trE = 0;
220  for (int dd = 0; dd < 3; dd++) {
221  trE += E(dd, dd);
222  }
223  S.resize(3, 3);
224  S.clear();
225  for (int dd = 0; dd < 3; dd++) {
226  S(dd, dd) = trE * lambda;
227  }
228  S += 2 * mu * E;
230  }
231 
232  /** \brief Function overload to implement user material
233  *
234 
235  * Calculation of Piola Kirchhoff I is implemented by user. Tangent matrix
236  * user implemented physical equation is calculated using automatic
237  * differentiation.
238 
239  * \f$\mathbf{S} =
240  \lambda\textrm{tr}[\mathbf{E}]\mathbf{I}+2\mu\mathbf{E}\f$
241 
242  * Notes: <br>
243  * Number of actual Gauss point is accessed from variable gG. <br>
244  * Access to operator data structures is available by variable opPtr. <br>
245  * Access to common data is by commonDataPtr. <br>
246 
247  * \param block_data used to give access to material parameters
248  * \param fe_ptr pointer to element data structures
249 
250  For details look to: <br>
251  NONLINEAR CONTINUUM MECHANICS FOR FINITE ELEMENT ANALYSIS, Javier Bonet,
252  Richard D. Wood
253 
254  */
256  const BlockData block_data,
257  boost::shared_ptr<const NumeredEntFiniteElement> fe_ptr) {
259  lambda = LAMBDA(block_data.E, block_data.PoissonRatio);
260  mu = MU(block_data.E, block_data.PoissonRatio);
264  P.resize(3, 3);
265  noalias(P) = prod(F, S);
267  }
268 
269  /**
270  * \brief add additional active variables
271  *
272  * \note This member function if used should be implement by template member
273  * function Specialization, different implementation needed for TYPE=double
274  * or TYPE=adouble
275  *
276  * More complex physical models depend on gradient of defamation and some
277  * additional variables. For example can depend on temperature. This
278  * function adds additional independent variables to the model.
279  *
280  * @param nb_active_variables number of active variables
281  * @return error code
282  */
283  virtual MoFEMErrorCode setUserActiveVariables(int &nb_active_variables) {
286  }
287 
288  /**
289  * \brief Add additional independent variables
290  * More complex physical models depend on gradient of defamation and some
291  * additional variables. For example can depend on temperature. This
292  * function adds additional independent variables to the model.
293  *
294  * /note First 9 elements are reserved for gradient of deformation.
295  * @param activeVariables vector of deepened variables, values after index
296  * 9 should be add.
297  *
298  * @return error code
299  */
300  virtual MoFEMErrorCode
304  }
305 
306  /** \brief Calculate elastic energy density
307  *
308  * \f[\Psi =
309  * \frac{1}{2}\lambda(\textrm{tr}[\mathbf{E}])^2+\mu\mathbf{E}:\mathbf{E}\f]
310  */
312  const BlockData block_data,
313  boost::shared_ptr<const NumeredEntFiniteElement> fe_ptr) {
315  lambda = LAMBDA(block_data.E, block_data.PoissonRatio);
316  mu = MU(block_data.E, block_data.PoissonRatio);
319  TYPE trace = 0;
320  eNergy = 0;
321  for (int ii = 0; ii < 3; ii++) {
322  trace += E(ii, ii);
323  for (int jj = 0; jj < 3; jj++) {
324  TYPE e = E(ii, jj);
325  eNergy += mu * e * e;
326  }
327  }
328  eNergy += 0.5 * lambda * trace * trace;
330  }
331 
332  /** \brief Calculate Eshelby stress
333  */
335  const BlockData block_data,
336  boost::shared_ptr<const NumeredEntFiniteElement> fe_ptr) {
338  CHKERR calculateP_PiolaKirchhoffI(block_data, fe_ptr);
339  CHKERR calculateElasticEnergy(block_data, fe_ptr);
340  SiGma.resize(3, 3, false);
341  noalias(SiGma) = -prod(trans(F), P);
342  for (int dd = 0; dd < 3; dd++) {
343  SiGma(dd, dd) += eNergy;
344  }
346  }
347 
348  /** \brief Do operations when pre-process
349  */
351  std::map<std::string, std::vector<VectorDouble>> &field_map,
352  std::map<std::string, std::vector<MatrixDouble>> &grad_map) {
355  }
356  };
357 
360 
361  std::vector<VectorDouble> &valuesAtGaussPts;
362  std::vector<MatrixDouble> &gradientAtGaussPts;
363  const EntityType zeroAtType;
364 
365  OpGetDataAtGaussPts(const std::string field_name,
366  std::vector<VectorDouble> &values_at_gauss_pts,
367  std::vector<MatrixDouble> &gradient_at_gauss_pts);
368 
369  /** \brief operator calculating deformation gradient
370  *
371  * temperature gradient is calculated multiplying derivatives of shape
372  * functions by degrees of freedom
373  */
374  MoFEMErrorCode doWork(int side, EntityType type,
376  };
377 
379  OpGetCommonDataAtGaussPts(const std::string field_name,
380  CommonData &common_data);
381  };
382 
383  /**
384  * \brief Operator performs automatic differentiation.
385  */
388 
389  BlockData &dAta; ///< Structure keeping data about problem, like material
390  ///< parameters
391  CommonData &commonData; ///< Structure keeping data abut this particular
392  ///< element, e.g. gradient of deformation at
393  ///< integration points
394  int tAg; //,lastId; ///< ADOL-C tag used for recording operations
395  int adlocReturnValue; ///< return value from ADOL-C, if non-zero that is
396  ///< error.
397  bool jAcobian; ///< if true Jacobian is calculated
398  bool fUnction; ///< if true stress i calculated
399  bool aLe; ///< true if arbitrary Lagrangian-Eulerian formulation
400  bool fieldDisp; ///< true if field of displacements is given, usually
401  ///< spatial positions are given.
402 
403  /**
404  \brief Construct operator to calculate Piola-Kirchhoff stress or its
405  derivatives over gradient deformation
406 
407  \param field_name approximation field name of spatial positions or
408  displacements \param data reference to block data (what is Young modulus,
409  Poisson ratio or what elements are part of the block) \param tag adol-c
410  unique tag of the tape \param jacobian if true derivative of Piola Stress
411  is calculated otherwise just stress is calculated \param field_disp if
412  true approximation field keeps displacements not spatial positions
413 
414  */
415  OpJacobianPiolaKirchhoffStress(const std::string field_name,
416  BlockData &data, CommonData &common_data,
417  int tag, bool jacobian, bool ale,
418  bool field_disp);
419 
422 
423  std::vector<MatrixDouble> *ptrh;
424  std::vector<MatrixDouble> *ptrH;
425 
426  /**
427  * \brief Calculate Paola-Kirchhoff I stress
428  * @return error code
429  */
430  virtual MoFEMErrorCode calculateStress(const int gg);
431 
432  /**
433  * \brief Record ADOL-C tape
434  * @return error code
435  */
436  virtual MoFEMErrorCode recordTag(const int gg);
437 
438  /**
439  * \brief Play ADOL-C tape
440  * @return error code
441  */
442  virtual MoFEMErrorCode playTag(const int gg);
443 
444  /**
445  * \brief Cgeck if tape is recorded for given integration point
446  * @param gg integration point
447  * @return true if tag is recorded
448  */
449  virtual bool recordTagForIntegrationPoint(const int gg) {
450  // return true;
451  if (gg == 0)
452  return true;
453  return false;
454  }
455 
456  /**
457  * \brief Calculate stress or jacobian at gauss points
458  *
459  * @param row_side
460  * @param row_type
461  * @param row_data
462  * @return error code
463  */
464  MoFEMErrorCode doWork(int row_side, EntityType row_type,
466  };
467 
468  /**
469  * \brief Calculate explicit derivative of free energy
470  */
473 
476 
477  int tAg; ///< tape tag
478  bool gRadient; ///< if set true gradient of energy is calculated
479  bool hEssian; ///< if set true hessian of energy is calculated
480  bool aLe; ///< true if arbitrary Lagrangian-Eulerian formulation
481  bool fieldDisp; ///< true if displacements instead spatial positions used
482 
483  OpJacobianEnergy(const std::string field_name, ///< field name for spatial
484  ///< positions or
485  ///< displacements
486  BlockData &data, CommonData &common_data, int tag,
487  bool gradient, bool hessian, bool ale, bool field_disp);
488 
491 
492  std::vector<MatrixDouble> *ptrh;
493  std::vector<MatrixDouble> *ptrH;
494 
495  /**
496  * \brief Check if tape is recorded for given integration point
497  * @param gg integration point
498  * @return true if tag is recorded
499  */
500  virtual bool recordTagForIntegrationPoint(const int gg) {
501  if (gg == 0)
502  return true;
503  return false;
504  }
505 
506  /**
507  * \brief Calculate Paola-Kirchhoff I stress
508  * @return error code
509  */
510  virtual MoFEMErrorCode calculateEnergy(const int gg);
511 
512  /**
513  * \brief Record ADOL-C tape
514  * @return error code
515  */
516  virtual MoFEMErrorCode recordTag(const int gg);
517 
518  /**
519  * \brief Play ADOL-C tape
520  * @return error code
521  */
522  virtual MoFEMErrorCode playTag(const int gg);
523 
524  MoFEMErrorCode doWork(int row_side, EntityType row_type,
526  };
527 
530 
533  bool fieldDisp;
534  bool aLe;
535 
536  ublas::vector<int> iNdices;
537  OpRhsPiolaKirchhoff(const std::string field_name, BlockData &data,
538  CommonData &common_data);
539 
541  MoFEMErrorCode doWork(int row_side, EntityType row_type,
543 
544  virtual MoFEMErrorCode aSemble(int row_side, EntityType row_type,
546  };
547 
548  struct OpEnergy
550 
553  Vec ghostVec;
554  bool fieldDisp;
555 
556  OpEnergy(const std::string field_name, BlockData &data,
557  CommonData &common_data, Vec ghost_vec, bool field_disp);
558  ~OpEnergy();
559 
560  MoFEMErrorCode doWork(int row_side, EntityType row_type,
562  };
563 
566 
569  int tAg;
570  bool aLe;
571 
572  ublas::vector<int> rowIndices;
573  ublas::vector<int> colIndices;
574 
575  OpLhsPiolaKirchhoff_dx(const std::string vel_field,
576  const std::string field_name, BlockData &data,
577  CommonData &common_data);
578 
580 
581  /**
582  \brief Directive of Piola Kirchhoff stress over spatial DOFs
583 
584  This project derivative \f$\frac{\partial P}{\partial F}\f$, that is
585  \f[
586  \frac{\partial P}{\partial x_\textrm{DOF}} = \frac{\partial P}{\partial
587  F}\frac{\partial F}{\partial x_\textrm{DOF}}, \f] where second therm
588  \f$\frac{\partial F}{\partial x_\textrm{DOF}}\f$ is derivative of shape
589  function
590 
591  */
593  int gg);
594 
595  virtual MoFEMErrorCode aSemble(int row_side, int col_side,
596  EntityType row_type, EntityType col_type,
599 
600  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
601  EntityType col_type,
604  };
605 
607 
608  OpLhsPiolaKirchhoff_dX(const std::string vel_field,
609  const std::string field_name, BlockData &data,
610  CommonData &common_data);
611 
612  /// \brief Derivative of Piola Kirchhoff stress over material DOFs
614 
615  MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type,
616  EntityType col_type,
619  };
620 
622 
623  OpJacobianEshelbyStress(const std::string field_name, BlockData &data,
624  CommonData &common_data, int tag, bool jacobian,
625  bool ale);
626 
627  MoFEMErrorCode calculateStress(const int gg);
628  };
629 
631 
632  OpRhsEshelbyStress(const std::string field_name, BlockData &data,
633  CommonData &common_data);
634  };
635 
636  /**
637  * \deprecated name with spelling mistake
638  */
640 
642 
643  OpLhsEshelby_dx(const std::string vel_field, const std::string field_name,
644  BlockData &data, CommonData &common_data);
645 
647  };
648 
650 
651  OpLhsEshelby_dX(const std::string vel_field, const std::string field_name,
652  BlockData &data, CommonData &common_data);
653 
655  };
656 
659  materialDoublePtr,
661  materialAdoublePtr);
662 
664  const std::string element_name,
665  const std::string spatial_position_field_name,
666  const std::string material_position_field_name = "MESH_NODE_POSITIONS",
667  const bool ale = false);
668 
669  /** \brief Set operators to calculate left hand tangent matrix and right hand
670  * residual
671  *
672  * \param fun class needed to calculate Piola Kirchhoff I Stress tensor
673  * \param spatial_position_field_name name of approximation field
674  * \param material_position_field_name name of field to define geometry
675  * \param ale true if arbitrary Lagrangian Eulerian formulation
676  * \param field_disp true if approximation field represents displacements
677  * otherwise it is field of spatial positions
678  */
680  const std::string spatial_position_field_name,
681  const std::string material_position_field_name = "MESH_NODE_POSITIONS",
682  const bool ale = false, const bool field_disp = false);
683 };
684 
685 #endif //__NONLINEAR_ELASTIC_HPP
686 
687 /**
688  * \defgroup nonlinear_elastic_elem NonLinear Elastic Element
689  * \ingroup user_modules
690  * \defgroup user_modules User modules
691  **/
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:506
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:74
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:482
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:513
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)
constexpr int order
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:66
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:601
#define MU(E, NU)
Definition: fem_tools.h:33
ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > h
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:412
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:72
OpGetDataAtGaussPts(const std::string field_name, std::vector< VectorDouble > &values_at_gauss_pts, std::vector< MatrixDouble > &gradient_at_gauss_pts)
data for calculation het 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)