MFEM  v3.3.2
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
nonlininteg.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #ifndef MFEM_NONLININTEG
13 #define MFEM_NONLININTEG
14 
15 #include "../config/config.hpp"
16 #include "fe.hpp"
17 #include "coefficient.hpp"
18 
19 namespace mfem
20 {
21 
22 /** The abstract base class NonlinearFormIntegrator is used to express the
23  local action of a general nonlinear finite element operator. In addition
24  it may provide the capability to assemble the local gradient operator
25  and to compute the local energy. */
27 {
28 protected:
30 
32  : IntRule(NULL) { }
33 
34 public:
35  /** @brief Prescribe a fixed IntegrationRule to use (when @a ir != NULL) or
36  let the integrator choose (when @a ir == NULL). */
37  void SetIntRule(const IntegrationRule *ir) { IntRule = ir; }
38 
39  /// Prescribe a fixed IntegrationRule to use.
40  void SetIntegrationRule(const IntegrationRule &irule) { IntRule = &irule; }
41 
42  /// Perform the local action of the NonlinearFormIntegrator
43  virtual void AssembleElementVector(const FiniteElement &el,
45  const Vector &elfun, Vector &elvect);
46 
47  /// @brief Perform the local action of the NonlinearFormIntegrator resulting
48  /// from a face integral term.
49  virtual void AssembleFaceVector(const FiniteElement &el1,
50  const FiniteElement &el2,
52  const Vector &elfun, Vector &elvect);
53 
54  /// Assemble the local gradient matrix
55  virtual void AssembleElementGrad(const FiniteElement &el,
57  const Vector &elfun, DenseMatrix &elmat);
58 
59  /// @brief Assemble the local action of the gradient of the
60  /// NonlinearFormIntegrator resulting from a face integral term.
61  virtual void AssembleFaceGrad(const FiniteElement &el1,
62  const FiniteElement &el2,
64  const Vector &elfun, DenseMatrix &elmat);
65 
66  /// Compute the local energy
67  virtual double GetElementEnergy(const FiniteElement &el,
69  const Vector &elfun);
70 
72 };
73 
74 
75 /// Abstract class for hyperelastic models
77 {
78 protected:
79  ElementTransformation *Ttr; /**< Reference-element to target-element
80  transformation. */
81 
82 public:
83  HyperelasticModel() : Ttr(NULL) { }
84  virtual ~HyperelasticModel() { }
85 
86  /// A reference-element to target-element transformation that can be used to
87  /// evaluate Coefficient%s.
88  /** @note It is assumed that _Ttr.SetIntPoint() is already called for the
89  point of interest. */
90  void SetTransformation(ElementTransformation &_Ttr) { Ttr = &_Ttr; }
91 
92  /** @brief Evaluate the strain energy density function, W = W(Jpt).
93  @param[in] Jpt Represents the target->physical transformation
94  Jacobian matrix. */
95  virtual double EvalW(const DenseMatrix &Jpt) const = 0;
96 
97  /** @brief Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
98  @param[in] Jpt Represents the target->physical transformation
99  Jacobian matrix.
100  @param[out] P The evaluated 1st Piola-Kirchhoff stress tensor. */
101  virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const = 0;
102 
103  /** @brief Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor
104  and assemble its contribution to the local gradient matrix 'A'.
105  @param[in] Jpt Represents the target->physical transformation
106  Jacobian matrix.
107  @param[in] DS Gradient of the basis matrix (dof x dim).
108  @param[in] weight Quadrature weight coefficient for the point.
109  @param[in,out] A Local gradient matrix where the contribution from this
110  point will be added.
111 
112  Computes weight * d(dW_dxi)_d(xj) at the current point, for all i and j,
113  where x1 ... xn are the FE dofs. This function is usually defined using
114  the matrix invariants and their derivatives.
115  */
116  virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
117  const double weight, DenseMatrix &A) const = 0;
118 };
119 
120 
121 /** Inverse-harmonic hyperelastic model with a strain energy density function
122  given by the formula: W(J) = (1/2) det(J) Tr((J J^t)^{-1}) where J is the
123  deformation gradient. */
125 {
126 protected:
127  mutable DenseMatrix Z, S; // dim x dim
128  mutable DenseMatrix G, C; // dof x dim
129 
130 public:
131  virtual double EvalW(const DenseMatrix &J) const;
132 
133  virtual void EvalP(const DenseMatrix &J, DenseMatrix &P) const;
134 
135  virtual void AssembleH(const DenseMatrix &J, const DenseMatrix &DS,
136  const double weight, DenseMatrix &A) const;
137 };
138 
139 
140 /** Neo-Hookean hyperelastic model with a strain energy density function given
141  by the formula: \f$(\mu/2)(\bar{I}_1 - dim) + (K/2)(det(J)/g - 1)^2\f$ where
142  J is the deformation gradient and \f$\bar{I}_1 = (det(J))^{-2/dim} Tr(J
143  J^t)\f$. The parameters \f$\mu\f$ and K are the shear and bulk moduli,
144  respectively, and g is a reference volumetric scaling. */
146 {
147 protected:
148  mutable double mu, K, g;
151 
152  mutable DenseMatrix Z; // dim x dim
153  mutable DenseMatrix G, C; // dof x dim
154 
155  inline void EvalCoeffs() const;
156 
157 public:
158  NeoHookeanModel(double _mu, double _K, double _g = 1.0)
159  : mu(_mu), K(_K), g(_g), have_coeffs(false) { c_mu = c_K = c_g = NULL; }
160 
162  : mu(0.0), K(0.0), g(1.0), c_mu(&_mu), c_K(&_K), c_g(_g),
163  have_coeffs(true) { }
164 
165  virtual double EvalW(const DenseMatrix &J) const;
166 
167  virtual void EvalP(const DenseMatrix &J, DenseMatrix &P) const;
168 
169  virtual void AssembleH(const DenseMatrix &J, const DenseMatrix &DS,
170  const double weight, DenseMatrix &A) const;
171 };
172 
173 
174 /** Hyperelastic integrator for any given HyperelasticModel.
175 
176  Represents @f$ \int W(Jpt) dx @f$ over a target zone, where W is the
177  @a model's strain energy density function, and Jpt is the Jacobian of the
178  target->physical coordinates transformation. The target configuration is
179  given by the current mesh at the time of the evaluation of the integrator.
180 */
182 {
183 private:
184  HyperelasticModel *model;
185 
186  // Jrt: the Jacobian of the target-to-reference-element transformation.
187  // Jpr: the Jacobian of the reference-to-physical-element transformation.
188  // Jpt: the Jacobian of the target-to-physical-element transformation.
189  // P: represents dW_d(Jtp) (dim x dim).
190  // DSh: gradients of reference shape functions (dof x dim).
191  // DS: gradients of the shape functions in the target (stress-free)
192  // configuration (dof x dim).
193  // PMatI: coordinates of the deformed configuration (dof x dim).
194  // PMatO: reshaped view into the local element contribution to the operator
195  // output - the result of AssembleElementVector() (dof x dim).
196  DenseMatrix DSh, DS, Jrt, Jpr, Jpt, P, PMatI, PMatO;
197 
198 public:
199  /** @param[in] m HyperelasticModel that will be integrated. */
201 
202  /** @brief Computes the integral of W(Jacobian(Trt)) over a target zone
203  @param[in] el Type of FiniteElement.
204  @param[in] Ttr Represents ref->target coordinates transformation.
205  @param[in] elfun Physical coordinates of the zone. */
206  virtual double GetElementEnergy(const FiniteElement &el,
208  const Vector &elfun);
209 
210  virtual void AssembleElementVector(const FiniteElement &el,
212  const Vector &elfun, Vector &elvect);
213 
214  virtual void AssembleElementGrad(const FiniteElement &el,
216  const Vector &elfun, DenseMatrix &elmat);
217 };
218 
219 }
220 
221 #endif
Abstract class for Finite Elements.
Definition: fe.hpp:47
virtual void EvalP(const DenseMatrix &J, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
virtual double EvalW(const DenseMatrix &J) const
Evaluate the strain energy density function, W = W(Jpt).
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:83
NeoHookeanModel(Coefficient &_mu, Coefficient &_K, Coefficient *_g=NULL)
virtual void AssembleH(const DenseMatrix &J, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: nonlininteg.cpp:85
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual double EvalW(const DenseMatrix &Jpt) const =0
Evaluate the strain energy density function, W = W(Jpt).
virtual void EvalP(const DenseMatrix &J, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: nonlininteg.cpp:66
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const =0
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
void EvalCoeffs() const
virtual void AssembleFaceGrad(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Tr, const Vector &elfun, DenseMatrix &elmat)
Assemble the local action of the gradient of the NonlinearFormIntegrator resulting from a face integr...
Definition: nonlininteg.cpp:41
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &Ttr, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun, Vector &elvect)
Perform the local action of the NonlinearFormIntegrator.
Definition: nonlininteg.cpp:17
virtual double GetElementEnergy(const FiniteElement &el, ElementTransformation &Ttr, const Vector &elfun)
Computes the integral of W(Jacobian(Trt)) over a target zone.
void SetTransformation(ElementTransformation &_Ttr)
Definition: nonlininteg.hpp:90
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
Definition: nonlininteg.cpp:33
NeoHookeanModel(double _mu, double _K, double _g=1.0)
Base class Coefficient that may optionally depend on time.
Definition: coefficient.hpp:31
virtual void AssembleFaceVector(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect)
Perform the local action of the NonlinearFormIntegrator resulting from a face integral term...
Definition: nonlininteg.cpp:25
ElementTransformation * Ttr
Definition: nonlininteg.hpp:79
const IntegrationRule * IntRule
Definition: nonlininteg.hpp:29
virtual void AssembleH(const DenseMatrix &J, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
HyperelasticNLFIntegrator(HyperelasticModel *m)
virtual double EvalW(const DenseMatrix &J) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: nonlininteg.cpp:59
void SetIntegrationRule(const IntegrationRule &irule)
Prescribe a fixed IntegrationRule to use.
Definition: nonlininteg.hpp:40
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const =0
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
void SetIntRule(const IntegrationRule *ir)
Prescribe a fixed IntegrationRule to use (when ir != NULL) or let the integrator choose (when ir == N...
Definition: nonlininteg.hpp:37
Abstract class for hyperelastic models.
Definition: nonlininteg.hpp:76
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &Ttr, const Vector &elfun, Vector &elvect)
Perform the local action of the NonlinearFormIntegrator.
Vector data type.
Definition: vector.hpp:41
NonlinearFormIntegrator(const IntegrationRule *ir=NULL)
Definition: nonlininteg.hpp:31
virtual double GetElementEnergy(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun)
Compute the local energy.
Definition: nonlininteg.cpp:51