MFEM  v3.4
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(ir) { }
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 /** The abstract base class BlockNonlinearFormIntegrator is
75  a generalization of the NonlinearFormIntegrator class suitable
76  for block state vectors. */
78 {
79 public:
80  /// Compute the local energy
81  virtual double GetElementEnergy(const Array<const FiniteElement *>&el,
83  const Array<const Vector *>&elfun);
84 
85  /// Perform the local action of the BlockNonlinearFormIntegrator
88  const Array<const Vector *> &elfun,
89  const Array<Vector *> &elvec);
90 
91  virtual void AssembleFaceVector(const Array<const FiniteElement *> &el1,
94  const Array<const Vector *> &elfun,
95  const Array<Vector *> &elvect);
96 
97  /// Assemble the local gradient matrix
98  virtual void AssembleElementGrad(const Array<const FiniteElement*> &el,
100  const Array<const Vector *> &elfun,
101  const Array2D<DenseMatrix *> &elmats);
102 
103  virtual void AssembleFaceGrad(const Array<const FiniteElement *>&el1,
106  const Array<const Vector *> &elfun,
107  const Array2D<DenseMatrix *> &elmats);
108 
110 };
111 
112 
113 /// Abstract class for hyperelastic models
115 {
116 protected:
117  ElementTransformation *Ttr; /**< Reference-element to target-element
118  transformation. */
119 
120 public:
121  HyperelasticModel() : Ttr(NULL) { }
122  virtual ~HyperelasticModel() { }
123 
124  /// A reference-element to target-element transformation that can be used to
125  /// evaluate Coefficient%s.
126  /** @note It is assumed that _Ttr.SetIntPoint() is already called for the
127  point of interest. */
128  void SetTransformation(ElementTransformation &_Ttr) { Ttr = &_Ttr; }
129 
130  /** @brief Evaluate the strain energy density function, W = W(Jpt).
131  @param[in] Jpt Represents the target->physical transformation
132  Jacobian matrix. */
133  virtual double EvalW(const DenseMatrix &Jpt) const = 0;
134 
135  /** @brief Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
136  @param[in] Jpt Represents the target->physical transformation
137  Jacobian matrix.
138  @param[out] P The evaluated 1st Piola-Kirchhoff stress tensor. */
139  virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const = 0;
140 
141  /** @brief Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor
142  and assemble its contribution to the local gradient matrix 'A'.
143  @param[in] Jpt Represents the target->physical transformation
144  Jacobian matrix.
145  @param[in] DS Gradient of the basis matrix (dof x dim).
146  @param[in] weight Quadrature weight coefficient for the point.
147  @param[in,out] A Local gradient matrix where the contribution from this
148  point will be added.
149 
150  Computes weight * d(dW_dxi)_d(xj) at the current point, for all i and j,
151  where x1 ... xn are the FE dofs. This function is usually defined using
152  the matrix invariants and their derivatives.
153  */
154  virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
155  const double weight, DenseMatrix &A) const = 0;
156 };
157 
158 
159 /** Inverse-harmonic hyperelastic model with a strain energy density function
160  given by the formula: W(J) = (1/2) det(J) Tr((J J^t)^{-1}) where J is the
161  deformation gradient. */
163 {
164 protected:
165  mutable DenseMatrix Z, S; // dim x dim
166  mutable DenseMatrix G, C; // dof x dim
167 
168 public:
169  virtual double EvalW(const DenseMatrix &J) const;
170 
171  virtual void EvalP(const DenseMatrix &J, DenseMatrix &P) const;
172 
173  virtual void AssembleH(const DenseMatrix &J, const DenseMatrix &DS,
174  const double weight, DenseMatrix &A) const;
175 };
176 
177 
178 /** Neo-Hookean hyperelastic model with a strain energy density function given
179  by the formula: \f$(\mu/2)(\bar{I}_1 - dim) + (K/2)(det(J)/g - 1)^2\f$ where
180  J is the deformation gradient and \f$\bar{I}_1 = (det(J))^{-2/dim} Tr(J
181  J^t)\f$. The parameters \f$\mu\f$ and K are the shear and bulk moduli,
182  respectively, and g is a reference volumetric scaling. */
184 {
185 protected:
186  mutable double mu, K, g;
189 
190  mutable DenseMatrix Z; // dim x dim
191  mutable DenseMatrix G, C; // dof x dim
192 
193  inline void EvalCoeffs() const;
194 
195 public:
196  NeoHookeanModel(double _mu, double _K, double _g = 1.0)
197  : mu(_mu), K(_K), g(_g), have_coeffs(false) { c_mu = c_K = c_g = NULL; }
198 
200  : mu(0.0), K(0.0), g(1.0), c_mu(&_mu), c_K(&_K), c_g(_g),
201  have_coeffs(true) { }
202 
203  virtual double EvalW(const DenseMatrix &J) const;
204 
205  virtual void EvalP(const DenseMatrix &J, DenseMatrix &P) const;
206 
207  virtual void AssembleH(const DenseMatrix &J, const DenseMatrix &DS,
208  const double weight, DenseMatrix &A) const;
209 };
210 
211 
212 /** Hyperelastic integrator for any given HyperelasticModel.
213 
214  Represents @f$ \int W(Jpt) dx @f$ over a target zone, where W is the
215  @a model's strain energy density function, and Jpt is the Jacobian of the
216  target->physical coordinates transformation. The target configuration is
217  given by the current mesh at the time of the evaluation of the integrator.
218 */
220 {
221 private:
222  HyperelasticModel *model;
223 
224  // Jrt: the Jacobian of the target-to-reference-element transformation.
225  // Jpr: the Jacobian of the reference-to-physical-element transformation.
226  // Jpt: the Jacobian of the target-to-physical-element transformation.
227  // P: represents dW_d(Jtp) (dim x dim).
228  // DSh: gradients of reference shape functions (dof x dim).
229  // DS: gradients of the shape functions in the target (stress-free)
230  // configuration (dof x dim).
231  // PMatI: coordinates of the deformed configuration (dof x dim).
232  // PMatO: reshaped view into the local element contribution to the operator
233  // output - the result of AssembleElementVector() (dof x dim).
234  DenseMatrix DSh, DS, Jrt, Jpr, Jpt, P, PMatI, PMatO;
235 
236 public:
237  /** @param[in] m HyperelasticModel that will be integrated. */
239 
240  /** @brief Computes the integral of W(Jacobian(Trt)) over a target zone
241  @param[in] el Type of FiniteElement.
242  @param[in] Ttr Represents ref->target coordinates transformation.
243  @param[in] elfun Physical coordinates of the zone. */
244  virtual double GetElementEnergy(const FiniteElement &el,
246  const Vector &elfun);
247 
248  virtual void AssembleElementVector(const FiniteElement &el,
250  const Vector &elfun, Vector &elvect);
251 
252  virtual void AssembleElementGrad(const FiniteElement &el,
254  const Vector &elfun, DenseMatrix &elmat);
255 };
256 
257 /** Hyperelastic incompressible Neo-Hookean integrator with the PK1 stress
258  \f$P = \mu F - p F^{-T}\f$ where \f$\mu\f$ is the shear modulus,
259  \f$p\f$ is the pressure, and \f$F\f$ is the deformation gradient */
261 {
262 private:
263  Coefficient *c_mu;
264  DenseMatrix DSh_u, DS_u, J0i, J, J1, Finv, P, F, FinvT;
265  DenseMatrix PMatI_u, PMatO_u, PMatI_p, PMatO_p, Z, G, C;
266  Vector Sh_p;
267 
268 public:
270 
271  virtual double GetElementEnergy(const Array<const FiniteElement *>&el,
273  const Array<const Vector *> &elfun);
274 
275  /// Perform the local action of the NonlinearFormIntegrator
278  const Array<const Vector *> &elfun,
279  const Array<Vector *> &elvec);
280 
281  /// Assemble the local gradient matrix
282  virtual void AssembleElementGrad(const Array<const FiniteElement*> &el,
284  const Array<const Vector *> &elfun,
285  const Array2D<DenseMatrix *> &elmats);
286 };
287 
288 }
289 
290 #endif
Abstract class for Finite Elements.
Definition: fe.hpp:140
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).
virtual void AssembleElementGrad(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun, const Array2D< DenseMatrix * > &elmats)
Assemble the local gradient matrix.
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:83
virtual void AssembleFaceVector(const Array< const FiniteElement * > &el1, const Array< const FiniteElement * > &el2, FaceElementTransformations &Tr, const Array< const Vector * > &elfun, const Array< Vector * > &elvect)
Definition: nonlininteg.cpp:69
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...
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).
virtual void AssembleFaceGrad(const Array< const FiniteElement * > &el1, const Array< const FiniteElement * > &el2, FaceElementTransformations &Tr, const Array< const Vector * > &elfun, const Array2D< DenseMatrix * > &elmats)
Definition: nonlininteg.cpp:90
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const =0
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
virtual void AssembleElementGrad(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun, const Array2D< DenseMatrix * > &elmats)
Assemble the local gradient matrix.
Definition: nonlininteg.cpp:80
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 double GetElementEnergy(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun)
Compute the local energy.
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.
Dynamic 2D array using row-major layout.
Definition: array.hpp:289
void SetTransformation(ElementTransformation &_Ttr)
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
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).
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...
virtual void AssembleElementVector(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun, const Array< Vector * > &elvec)
Perform the local action of the NonlinearFormIntegrator.
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.
IncompressibleNeoHookeanIntegrator(Coefficient &_mu)
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &Ttr, const Vector &elfun, Vector &elvect)
Perform the local action of the NonlinearFormIntegrator.
virtual double GetElementEnergy(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun)
Compute the local energy.
Vector data type.
Definition: vector.hpp:48
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:50
virtual void AssembleElementVector(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun, const Array< Vector * > &elvec)
Perform the local action of the BlockNonlinearFormIntegrator.
Definition: nonlininteg.cpp:59