MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
nonlininteg.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12#ifndef MFEM_NONLININTEG
13#define MFEM_NONLININTEG
14
15#include "../config/config.hpp"
16#include "fe.hpp"
17#include "coefficient.hpp"
18#include "fespace.hpp"
20#include "integrator.hpp"
21
22namespace mfem
23{
24
25/** @brief This class is used to express the local action of a general nonlinear
26 finite element operator. In addition it may provide the capability to
27 assemble the local gradient operator and to compute the local energy. */
29{
30public:
31 enum Mode
32 {
33 ELEMENTWISE = 0, /**< Element-wise integration (default) */
34 PATCHWISE = 1, /**< Patch-wise integration (NURBS meshes) */
35 PATCHWISE_REDUCED = 2, /**< Patch-wise integration (NURBS meshes) with
36 reduced integration rules. */
37 };
38
39protected:
41
42 // CEED extension
44
46
48 : Integrator(ir), ceedOp(NULL) { }
49
50public:
51
53
54
55 bool Patchwise() const { return integrationMode != Mode::ELEMENTWISE; }
56
57 /// Set the memory type used for GeometricFactors and other large allocations
58 /// in PA extensions.
59 void SetPAMemoryType(MemoryType mt) { pa_mt = mt; }
60
61
62 /// Perform the local action of the NonlinearFormIntegrator
63 virtual void AssembleElementVector(const FiniteElement &el,
65 const Vector &elfun, Vector &elvect);
66
67 /// @brief Perform the local action of the NonlinearFormIntegrator resulting
68 /// from a face integral term.
69 virtual void AssembleFaceVector(const FiniteElement &el1,
70 const FiniteElement &el2,
72 const Vector &elfun, Vector &elvect);
73
74 /// Assemble the local gradient matrix
75 virtual void AssembleElementGrad(const FiniteElement &el,
77 const Vector &elfun, DenseMatrix &elmat);
78
79 /// @brief Assemble the local action of the gradient of the
80 /// NonlinearFormIntegrator resulting from a face integral term.
81 virtual void AssembleFaceGrad(const FiniteElement &el1,
82 const FiniteElement &el2,
84 const Vector &elfun, DenseMatrix &elmat);
85
86 /// Compute the local energy
87 virtual real_t GetElementEnergy(const FiniteElement &el,
89 const Vector &elfun);
90
91 /// Method defining partial assembly.
92 /** The result of the partial assembly is stored internally so that it can be
93 used later in the methods AddMultPA(). */
94 virtual void AssemblePA(const FiniteElementSpace &fes);
95
96 /** The result of the partial assembly is stored internally so that it can be
97 used later in the methods AddMultPA().
98 Used with BilinearFormIntegrators that have different spaces. */
99 virtual void AssemblePA(const FiniteElementSpace &trial_fes,
100 const FiniteElementSpace &test_fes);
101
102 /** @brief Prepare the integrator for partial assembly (PA) gradient
103 evaluations on the given FE space @a fes at the state @a x. */
104 /** The result of the partial assembly is stored internally so that it can be
105 used later in the methods AddMultGradPA() and AssembleGradDiagonalPA().
106 The state Vector @a x is an E-vector. */
107 virtual void AssembleGradPA(const Vector &x, const FiniteElementSpace &fes);
108
109 /// Compute the local (to the MPI rank) energy with partial assembly.
110 /** Here the state @a x is an E-vector. This method can be called only after
111 the method AssemblePA() has been called. */
112 virtual real_t GetLocalStateEnergyPA(const Vector &x) const;
113
114 /// Method for partially assembled action.
115 /** Perform the action of integrator on the input @a x and add the result to
116 the output @a y. Both @a x and @a y are E-vectors, i.e. they represent
117 the element-wise discontinuous version of the FE space.
118
119 This method can be called only after the method AssemblePA() has been
120 called. */
121 virtual void AddMultPA(const Vector &x, Vector &y) const;
122
123 /// Method for partially assembled gradient action.
124 /** All arguments are E-vectors. This method can be called only after the
125 method AssembleGradPA() has been called.
126
127 @param[in] x The gradient Operator is applied to the Vector @a x.
128 @param[in,out] y The result Vector: $ y += G x $. */
129 virtual void AddMultGradPA(const Vector &x, Vector &y) const;
130
131 /// Method for computing the diagonal of the gradient with partial assembly.
132 /** The result Vector @a diag is an E-Vector. This method can be called only
133 after the method AssembleGradPA() has been called.
134
135 @param[in,out] diag The result Vector: $ diag += diag(G) $. */
136 virtual void AssembleGradDiagonalPA(Vector &diag) const;
137
138 /// Indicates whether this integrator can use a Ceed backend.
139 virtual bool SupportsCeed() const { return false; }
140
141 /// Method defining fully unassembled operator.
142 virtual void AssembleMF(const FiniteElementSpace &fes);
143
144 /** Perform the action of integrator on the input @a x and add the result to
145 the output @a y. Both @a x and @a y are E-vectors, i.e. they represent
146 the element-wise discontinuous version of the FE space.
147
148 This method can be called only after the method AssembleMF() has been
149 called. */
150 virtual void AddMultMF(const Vector &x, Vector &y) const;
151
153
155 {
156 delete ceedOp;
157 }
158};
159
160/** The abstract base class BlockNonlinearFormIntegrator is
161 a generalization of the NonlinearFormIntegrator class suitable
162 for block state vectors. */
164{
165public:
166 /// Compute the local energy
169 const Array<const Vector *>&elfun);
170
171 /// Perform the local action of the BlockNonlinearFormIntegrator
174 const Array<const Vector *> &elfun,
175 const Array<Vector *> &elvec);
176
177 virtual void AssembleFaceVector(const Array<const FiniteElement *> &el1,
180 const Array<const Vector *> &elfun,
181 const Array<Vector *> &elvect);
182
183 /// Assemble the local gradient matrix
186 const Array<const Vector *> &elfun,
187 const Array2D<DenseMatrix *> &elmats);
188
189 virtual void AssembleFaceGrad(const Array<const FiniteElement *>&el1,
192 const Array<const Vector *> &elfun,
193 const Array2D<DenseMatrix *> &elmats);
194
196};
197
198
199/// Abstract class for hyperelastic models
201{
202protected:
203 ElementTransformation *Ttr; /**< Reference-element to target-element
204 transformation. */
205
206public:
207 HyperelasticModel() : Ttr(NULL) { }
208 virtual ~HyperelasticModel() { }
209
210 /// A reference-element to target-element transformation that can be used to
211 /// evaluate Coefficient%s.
212 /** @note It is assumed that Ttr_.SetIntPoint() is already called for the
213 point of interest. */
215
216 /** @brief Evaluate the strain energy density function, W = W(Jpt).
217 @param[in] Jpt Represents the target->physical transformation
218 Jacobian matrix. */
219 virtual real_t EvalW(const DenseMatrix &Jpt) const = 0;
220
221 /** @brief Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
222 @param[in] Jpt Represents the target->physical transformation
223 Jacobian matrix.
224 @param[out] P The evaluated 1st Piola-Kirchhoff stress tensor. */
225 virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const = 0;
226
227 /** @brief Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor
228 and assemble its contribution to the local gradient matrix 'A'.
229 @param[in] Jpt Represents the target->physical transformation
230 Jacobian matrix.
231 @param[in] DS Gradient of the basis matrix (dof x dim).
232 @param[in] weight Quadrature weight coefficient for the point.
233 @param[in,out] A Local gradient matrix where the contribution from this
234 point will be added.
235
236 Computes weight * d(dW_dxi)_d(xj) at the current point, for all i and j,
237 where x1 ... xn are the FE dofs. This function is usually defined using
238 the matrix invariants and their derivatives.
239 */
240 virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
241 const real_t weight, DenseMatrix &A) const = 0;
242};
243
244
245/** Inverse-harmonic hyperelastic model with a strain energy density function
246 given by the formula: W(J) = (1/2) det(J) Tr((J J^t)^{-1}) where J is the
247 deformation gradient. */
249{
250protected:
251 mutable DenseMatrix Z, S; // dim x dim
252 mutable DenseMatrix G, C; // dof x dim
253
254public:
255 real_t EvalW(const DenseMatrix &J) const override;
256
257 void EvalP(const DenseMatrix &J, DenseMatrix &P) const override;
258
259 void AssembleH(const DenseMatrix &J, const DenseMatrix &DS,
260 const real_t weight, DenseMatrix &A) const override;
261};
262
263
264/** Neo-Hookean hyperelastic model with a strain energy density function given
265 by the formula: $(\mu/2)(\bar{I}_1 - dim) + (K/2)(det(J)/g - 1)^2$ where
266 J is the deformation gradient and $$\bar{I}_1 = (det(J))^{-2/dim} Tr(J
267 J^t)$$. The parameters $\mu$ and K are the shear and bulk moduli,
268 respectively, and g is a reference volumetric scaling. */
270{
271protected:
272 mutable real_t mu, K, g;
275
276 mutable DenseMatrix Z; // dim x dim
277 mutable DenseMatrix G, C; // dof x dim
278
279 inline void EvalCoeffs() const;
280
281public:
283 : mu(mu_), K(K_), g(g_), have_coeffs(false) { c_mu = c_K = c_g = NULL; }
284
286 : mu(0.0), K(0.0), g(1.0), c_mu(&mu_), c_K(&K_), c_g(g_),
287 have_coeffs(true) { }
288
289 real_t EvalW(const DenseMatrix &J) const override;
290
291 void EvalP(const DenseMatrix &J, DenseMatrix &P) const override;
292
293 void AssembleH(const DenseMatrix &J, const DenseMatrix &DS,
294 const real_t weight, DenseMatrix &A) const override;
295};
296
297
298/** Hyperelastic integrator for any given HyperelasticModel.
299
300 Represents $ \int W(Jpt) dx $ over a target zone, where W is the
301 @a model's strain energy density function, and Jpt is the Jacobian of the
302 target->physical coordinates transformation. The target configuration is
303 given by the current mesh at the time of the evaluation of the integrator.
304*/
306{
307private:
308 HyperelasticModel *model;
309
310 // Jrt: the Jacobian of the target-to-reference-element transformation.
311 // Jpr: the Jacobian of the reference-to-physical-element transformation.
312 // Jpt: the Jacobian of the target-to-physical-element transformation.
313 // P: represents dW_d(Jtp) (dim x dim).
314 // DSh: gradients of reference shape functions (dof x dim).
315 // DS: gradients of the shape functions in the target (stress-free)
316 // configuration (dof x dim).
317 // PMatI: coordinates of the deformed configuration (dof x dim).
318 // PMatO: reshaped view into the local element contribution to the operator
319 // output - the result of AssembleElementVector() (dof x dim).
320 DenseMatrix DSh, DS, Jrt, Jpr, Jpt, P, PMatI, PMatO;
321
322public:
323 /** @param[in] m HyperelasticModel that will be integrated. */
325
326 /** @brief Computes the integral of W(Jacobian(Trt)) over a target zone
327 @param[in] el Type of FiniteElement.
328 @param[in] Ttr Represents ref->target coordinates transformation.
329 @param[in] elfun Physical coordinates of the zone. */
332 const Vector &elfun) override;
333
336 const Vector &elfun, Vector &elvect) override;
337
338 void AssembleElementGrad(const FiniteElement &el,
340 const Vector &elfun, DenseMatrix &elmat) override;
341protected:
343 const FiniteElement& trial_fe,
344 const FiniteElement& test_fe,
345 const ElementTransformation& trans) const override;
346};
347
348/** Hyperelastic incompressible Neo-Hookean integrator with the PK1 stress
349 $P = \mu F - p F^{-T}$ where $\mu$ is the shear modulus,
350 $p$ is the pressure, and $F$ is the deformation gradient */
352{
353private:
354 Coefficient *c_mu;
355 DenseMatrix DSh_u, DS_u, J0i, J, J1, Finv, P, F, FinvT;
356 DenseMatrix PMatI_u, PMatO_u, PMatI_p, PMatO_p, Z, G, C;
357 Vector Sh_p;
358
359public:
361
364 const Array<const Vector *> &elfun) override;
365
366 /// Perform the local action of the NonlinearFormIntegrator
369 const Array<const Vector *> &elfun,
370 const Array<Vector *> &elvec) override;
371
372 /// Assemble the local gradient matrix
375 const Array<const Vector *> &elfun,
376 const Array2D<DenseMatrix *> &elmats) override;
377};
378
379
381{
382private:
383 Coefficient *Q{};
384 DenseMatrix dshape, dshapex, EF, gradEF, ELV, elmat_comp;
385 Vector shape;
386 // PA extension
387 Vector pa_data;
388 const DofToQuad *maps; ///< Not owned
389 const GeometricFactors *geom; ///< Not owned
390 int dim, ne, nq;
391
392public:
394
396
397 static const IntegrationRule &GetRule(const FiniteElement &fe,
398 const ElementTransformation &T);
399
402 const Vector &elfun,
403 Vector &elvect) override;
404
405 void AssembleElementGrad(const FiniteElement &el,
407 const Vector &elfun,
408 DenseMatrix &elmat) override;
409
411
412 void AssemblePA(const FiniteElementSpace &fes) override;
413
414 void AssembleMF(const FiniteElementSpace &fes) override;
415
416 void AddMultPA(const Vector &x, Vector &y) const override;
417
418 void AddMultMF(const Vector &x, Vector &y) const override;
419
420
421protected:
423 const FiniteElement& trial_fe,
424 const FiniteElement& test_fe,
425 const ElementTransformation& trans) const override
426 {
427 return &GetRule(test_fe, trans);
428 }
429};
430
431
432/** This class is used to assemble the convective form of the nonlinear term
433 arising in the Navier-Stokes equations $(u \cdot \nabla v, w )$ */
436{
437private:
438 Coefficient *Q{};
439 DenseMatrix dshape, dshapex, EF, gradEF, ELV, elmat_comp;
440 Vector shape;
441
442public:
444
446
447 void AssembleElementGrad(const FiniteElement &el,
449 const Vector &elfun,
450 DenseMatrix &elmat) override;
451};
452
453
454/** This class is used to assemble the skew-symmetric form of the nonlinear term
455 arising in the Navier-Stokes equations
456 $.5*(u \cdot \nabla v, w ) - .5*(u \cdot \nabla w, v )$ */
459{
460private:
461 Coefficient *Q{};
462 DenseMatrix dshape, dshapex, EF, gradEF, ELV, elmat_comp;
463 Vector shape;
464
465public:
467
469
470 void AssembleElementGrad(const FiniteElement &el,
472 const Vector &elfun,
473 DenseMatrix &elmat) override;
474};
475
476}
477
478#endif
Dynamic 2D array using row-major layout.
Definition array.hpp:392
virtual real_t GetElementEnergy(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun)
Compute the local energy.
virtual void AssembleFaceVector(const Array< const FiniteElement * > &el1, const Array< const FiniteElement * > &el2, FaceElementTransformations &Tr, const Array< const Vector * > &elfun, const Array< Vector * > &elvect)
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.
virtual void AssembleFaceGrad(const Array< const FiniteElement * > &el1, const Array< const FiniteElement * > &el2, FaceElementTransformations &Tr, const Array< const Vector * > &elfun, const Array2D< DenseMatrix * > &elmats)
virtual void AssembleElementGrad(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun, const Array2D< DenseMatrix * > &elmats)
Assemble the local gradient matrix.
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
void AssembleElementGrad(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, DenseMatrix &elmat) override
Assemble the local gradient matrix.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
Structure representing the matrices/tensors needed to evaluate (in reference space) the values,...
Definition fe_base.hpp:141
A specialized ElementTransformation class representing a face and its two neighboring elements.
Definition eltrans.hpp:750
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:244
Abstract class for all finite elements.
Definition fe_base.hpp:244
Structure for storing mesh geometric factors: coordinates, Jacobians, and determinants of the Jacobia...
Definition mesh.hpp:2937
Abstract class for hyperelastic models.
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const real_t weight, DenseMatrix &A) const =0
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
virtual real_t EvalW(const DenseMatrix &Jpt) const =0
Evaluate the strain energy density function, W = W(Jpt).
ElementTransformation * Ttr
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const =0
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
void SetTransformation(ElementTransformation &Ttr_)
void AssembleElementVector(const FiniteElement &el, ElementTransformation &Ttr, const Vector &elfun, Vector &elvect) override
Perform the local action of the NonlinearFormIntegrator.
real_t GetElementEnergy(const FiniteElement &el, ElementTransformation &Ttr, const Vector &elfun) override
Computes the integral of W(Jacobian(Trt)) over a target zone.
const IntegrationRule * GetDefaultIntegrationRule(const FiniteElement &trial_fe, const FiniteElement &test_fe, const ElementTransformation &trans) const override
Subclasses should override to choose a default integration rule.
HyperelasticNLFIntegrator(HyperelasticModel *m)
void AssembleElementGrad(const FiniteElement &el, ElementTransformation &Ttr, const Vector &elfun, DenseMatrix &elmat) override
Assemble the local gradient matrix.
real_t GetElementEnergy(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun) override
Compute the local energy.
void AssembleElementVector(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun, const Array< Vector * > &elvec) override
Perform the local action of the NonlinearFormIntegrator.
void AssembleElementGrad(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun, const Array2D< DenseMatrix * > &elmats) override
Assemble the local gradient matrix.
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
This base class implements some shared functionality between linear and nonlinear form integrators.
void EvalP(const DenseMatrix &J, DenseMatrix &P) const override
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
void AssembleH(const DenseMatrix &J, const DenseMatrix &DS, const real_t weight, DenseMatrix &A) const override
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
real_t EvalW(const DenseMatrix &J) const override
Evaluate the strain energy density function, W = W(Jpt).
void AssembleH(const DenseMatrix &J, const DenseMatrix &DS, const real_t weight, DenseMatrix &A) const override
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
NeoHookeanModel(real_t mu_, real_t K_, real_t g_=1.0)
real_t EvalW(const DenseMatrix &J) const override
Evaluate the strain energy density function, W = W(Jpt).
void EvalP(const DenseMatrix &J, DenseMatrix &P) const override
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
NeoHookeanModel(Coefficient &mu_, Coefficient &K_, Coefficient *g_=NULL)
This class is used to express the local action of a general nonlinear finite element operator....
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun, Vector &elvect)
Perform the local action of the NonlinearFormIntegrator.
virtual real_t GetLocalStateEnergyPA(const Vector &x) const
Compute the local (to the MPI rank) energy with partial assembly.
virtual void AddMultMF(const Vector &x, Vector &y) const
virtual void AssembleMF(const FiniteElementSpace &fes)
Method defining fully unassembled operator.
virtual void AssembleGradPA(const Vector &x, const FiniteElementSpace &fes)
Prepare the integrator for partial assembly (PA) gradient evaluations on the given FE space fes at th...
NonlinearFormIntegrator(const IntegrationRule *ir=NULL)
virtual void AssemblePA(const FiniteElementSpace &fes)
Method defining partial assembly.
void SetPAMemoryType(MemoryType mt)
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...
virtual void AssembleGradDiagonalPA(Vector &diag) const
Method for computing the diagonal of the gradient with partial assembly.
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
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.
virtual void AddMultPA(const Vector &x, Vector &y) const
Method for partially assembled action.
virtual void AddMultGradPA(const Vector &x, Vector &y) const
Method for partially assembled gradient action.
virtual real_t GetElementEnergy(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun)
Compute the local energy.
virtual bool SupportsCeed() const
Indicates whether this integrator can use a Ceed backend.
void AssembleElementGrad(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, DenseMatrix &elmat) override
Assemble the local gradient matrix.
void AddMultPA(const Vector &x, Vector &y) const override
Method for partially assembled action.
void AssembleElementVector(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, Vector &elvect) override
Perform the local action of the NonlinearFormIntegrator.
void AssemblePA(const FiniteElementSpace &fes) override
Method defining partial assembly.
void AssembleMF(const FiniteElementSpace &fes) override
Method defining fully unassembled operator.
void AddMultMF(const Vector &x, Vector &y) const override
void AssembleElementGrad(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, DenseMatrix &elmat) override
Assemble the local gradient matrix.
static const IntegrationRule & GetRule(const FiniteElement &fe, const ElementTransformation &T)
const IntegrationRule * GetDefaultIntegrationRule(const FiniteElement &trial_fe, const FiniteElement &test_fe, const ElementTransformation &trans) const override
Subclasses should override to choose a default integration rule.
Vector data type.
Definition vector.hpp:82
void trans(const Vector &u, Vector &x)
Definition ex27.cpp:412
float real_t
Definition config.hpp:43
MemoryType
Memory types supported by MFEM.