MFEM  v4.5.2
Finite element discretization library
elasticity_operator.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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_ELASTICITY_OP_HPP
13 #define MFEM_ELASTICITY_OP_HPP
14 
15 #include "../kernels/elasticity_kernels.hpp"
16 #include "mfem.hpp"
17 
18 namespace mfem
19 {
21 {
22 public:
23  ElasticityOperator(ParMesh &mesh, const int order);
24 
25  /**
26  * @brief Compute the residual Y = R(U) representing the elasticity equation
27  * with a material model chosen by calling SetMaterial.
28  *
29  * The output vector @a Y has essential degrees of freedom applied by setting
30  * them to zero. This ensures R(U)_i = 0 being satisfied for each essential
31  * dof i.
32  *
33  * @param U U
34  * @param Y Residual R(U)
35  */
36  virtual void Mult(const Vector &U, Vector &Y) const override;
37 
38  /**
39  * @brief Get the Gradient object
40  *
41  * Update and cache the state vector @a U, used to compute the linearization
42  * dR(U)/dU.
43  *
44  * @param U
45  * @return Operator&
46  */
47  Operator &GetGradient(const Vector &U) const override;
48 
49  /**
50  * @brief Multiply the linearization of the residual R(U) wrt to the current
51  * state U by a perturbation @a dX.
52  *
53  * Y = dR(U)/dU * dX = K(U) dX
54  *
55  * @param dX
56  * @param Y
57  */
58  void GradientMult(const Vector &dX, Vector &Y) const;
59 
60  /**
61  * @brief Assemble the linearization of the residual K = dR(U)/dU.
62  *
63  * This method needs three input vectors which also act as output vectors.
64  * They don't have to be the right size on the first call, but it is advised
65  * that memory is kept alive during successive call. The data layout of the
66  * outputs will be
67  *
68  * @a Ke_diag: dofs x dofs x dofs x dim x ne x dim
69  *
70  * @a K_diag_local: width(H1_Restriction) x dim
71  *
72  * @a K_diag: width(H1_Prolongation) x dim
73  *
74  * This data layout is needed due to the Ordering::byNODES. See method
75  * implementation comments for more details. The output @a K_diag has
76  * modified entries when essential boundaries are defined. Each essential dof
77  * row and column are set to zero with it's diagonal entry set to 1.
78  *
79  * @param Ke_diag
80  * @param K_diag_local
81  * @param K_diag
82  */
83  void AssembleGradientDiagonal(Vector &Ke_diag, Vector &K_diag_local,
84  Vector &K_diag) const;
85 
87 
89  /// Polynomial order of the FE space
90  const int order_;
91  const int dim_;
92  const int vdim_;
93  /// Number of elements in the mesh (rank local)
94  const int ne_;
95  /// H1 finite element collection
97  /// H1 finite element space
99  // Integration rule
100  IntegrationRule *ir_ = nullptr;
101  /// Number of degrees of freedom in 1D
102  int d1d_;
103  /// Number of quadrature points in 1D
104  int q1d_;
111  const DofToQuad *maps_;
112  /// Input state L-vector
113  mutable Vector X_local_;
114  /// Input state E-vector
115  mutable Vector X_el_;
116  /// Output state L-vector
117  mutable Vector Y_local_;
118  /// Output state E-Vector
119  mutable Vector Y_el_;
120  /// Cached current state. Used to determine the state on which to compute the
121  /// linearization on during the Newton method.
125  /// Temporary vector for the perturbation of the solution with essential
126  /// boundaries eliminated. Defined as a T-vector.
127  mutable Vector dX_ess_;
128 
129  /// Flag to enable caching of the gradient. If enabled, during linear
130  /// iterations the operator only applies the gradient on each quadrature
131  /// point rather than recompute the action.
132  bool use_cache_ = true;
133  mutable bool recompute_cache_ = false;
135 
136  /**
137  * @brief Wrapper for the application of the residual R(U).
138  *
139  * The wrapper is used in SetMaterial to instantiate the chosen kernel and
140  * erase the material type kernel. This is purely an interface design choice
141  * and could be replaced by an abstract base class for the material including
142  * virtual function calls.
143  */
144  std::function<void(const int, const Array<double> &, const Array<double> &,
145  const Array<double> &, const Vector &, const Vector &,
146  const Vector &, Vector &)>
148 
149  /**
150  * @brief Wrapper for the application of the gradient of the residual
151  *
152  * K(U) dX = dR(U)/dU dX
153  */
154  std::function<void(const int, const Array<double> &, const Array<double> &,
155  const Array<double> &, const Vector &, const Vector &,
156  const Vector &, Vector &, const Vector &)>
158 
159  /**
160  * @brief Wrapper for the assembly of the gradient on each diagonal element
161  *
162  * Ke_ii(U) = dRe_ii(U)/dU
163  */
164  std::function<void(const int, const Array<double> &, const Array<double> &,
165  const Array<double> &, const Vector &, const Vector &,
166  const Vector &, Vector &)>
168 
169  /**
170  * @brief Set the material type.
171  *
172  * This method sets the material type by instantiating the kernels with a
173  * material_type object.
174  *
175  * @tparam material_type
176  * @param[in] material
177  */
178  template <typename material_type>
179  void SetMaterial(const material_type &material)
180  {
181  if (dim_ != 3)
182  {
183  MFEM_ABORT("dim != 3 not implemented");
184  }
185 
187  [=](const int ne, const Array<double> &B_, const Array<double> &G_,
188  const Array<double> &W_, const Vector &Jacobian_,
189  const Vector &detJ_, const Vector &X_, Vector &Y_)
190  {
191  const int id = (d1d_ << 4) | q1d_;
192  switch (id)
193  {
194  case 0x22:
195  {
196  ElasticityKernels::Apply3D<2, 2, material_type>(
197  ne, B_, G_, W_, Jacobian_, detJ_, X_, Y_, material);
198  break;
199  }
200  case 0x33:
201  {
202  ElasticityKernels::Apply3D<3, 3, material_type>(
203  ne, B_, G_, W_, Jacobian_, detJ_, X_, Y_, material);
204  break;
205  }
206  case 0x44:
207  ElasticityKernels::Apply3D<4, 4, material_type>(
208  ne, B_, G_, W_, Jacobian_, detJ_, X_, Y_, material);
209  break;
210  default:
211  MFEM_ABORT("Not implemented: " << std::hex << id << std::dec);
212  }
213  };
214 
216  [=](const int ne, const Array<double> &B_, const Array<double> &G_,
217  const Array<double> &W_, const Vector &Jacobian_,
218  const Vector &detJ_, const Vector &dU_, Vector &dF_,
219  const Vector &U_)
220  {
221  const int id = (d1d_ << 4) | q1d_;
222  switch (id)
223  {
224  case 0x22:
225  ElasticityKernels::ApplyGradient3D<2, 2, material_type>(
226  ne, B_, G_, W_, Jacobian_, detJ_, dU_, dF_, U_, material,
228  break;
229  case 0x33:
230  ElasticityKernels::ApplyGradient3D<3, 3, material_type>(
231  ne, B_, G_, W_, Jacobian_, detJ_, dU_, dF_, U_, material,
233  break;
234  case 0x44:
235  ElasticityKernels::ApplyGradient3D<4, 4, material_type>(
236  ne, B_, G_, W_, Jacobian_, detJ_, dU_, dF_, U_, material,
238  break;
239  default:
240  MFEM_ABORT("Not implemented for D1D=" << d1d_ << " and Q1D=" << q1d_);
241  }
242  };
243 
245  [=](const int ne, const Array<double> &B_, const Array<double> &G_,
246  const Array<double> &W_, const Vector &Jacobian_,
247  const Vector &detJ_, const Vector &X_, Vector &Y_)
248  {
249  const int id = (d1d_ << 4) | q1d_;
250  switch (id)
251  {
252  case 0x22:
253  ElasticityKernels::AssembleGradientDiagonal3D<2, 2, material_type>(
254  ne, B_, G_, W_, Jacobian_, detJ_, X_, Y_, material);
255  break;
256  case 0x33:
257  ElasticityKernels::AssembleGradientDiagonal3D<3, 3, material_type>(
258  ne, B_, G_, W_, Jacobian_, detJ_, X_, Y_, material);
259  break;
260  case 0x44:
261  ElasticityKernels::AssembleGradientDiagonal3D<4, 4, material_type>(
262  ne, B_, G_, W_, Jacobian_, detJ_, X_, Y_, material);
263  break;
264  default:
265  MFEM_ABORT("Not implemented: " << std::hex << id << std::dec);
266  }
267  };
268  }
269 
270  /**
271  * @brief Set the essential attributes which mark degrees of freedom for the
272  * solving process.
273  *
274  * Can be either a fixed boundary or a prescribed displacement.
275  *
276  * @param attr
277  */
279  {
281  }
282 
283  /**
284  * @brief Set the attributes which mark the degrees of freedom that have a
285  * fixed displacement.
286  *
287  * @param[in] attr
288  */
290  {
292  }
293 
294  /**
295  * @brief Return the T-vector degrees of freedom that have been marked as
296  * displaced.
297  *
298  * @return T-vector degrees of freedom that have been marked as displaced
299  */
301 };
302 
303 } // namespace mfem
304 
305 #endif
std::function< void(const int, const Array< double > &, const Array< double > &, const Array< double > &, const Vector &, const Vector &, const Vector &, Vector &)> element_kernel_assemble_diagonal_wrapper
Wrapper for the assembly of the gradient on each diagonal element.
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:1032
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
const GeometricFactors * geometric_factors_
void SetMaterial(const material_type &material)
Set the material type.
Vector X_el_
Input state E-vector.
const int order_
Polynomial order of the FE space.
virtual void Mult(const Vector &U, Vector &Y) const override
Compute the residual Y = R(U) representing the elasticity equation with a material model chosen by ca...
ElasticityOperator(ParMesh &mesh, const int order)
Abstract parallel finite element space.
Definition: pfespace.hpp:28
Structure for storing mesh geometric factors: coordinates, Jacobians, and determinants of the Jacobia...
Definition: mesh.hpp:1915
void AssembleGradientDiagonal(Vector &Ke_diag, Vector &K_diag_local, Vector &K_diag) const
Assemble the linearization of the residual K = dR(U)/dU.
std::function< void(const int, const Array< double > &, const Array< double > &, const Array< double > &, const Vector &, const Vector &, const Vector &, Vector &, const Vector &)> element_apply_gradient_kernel_wrapper
Wrapper for the application of the gradient of the residual.
int material(Vector &x, Vector &xmin, Vector &xmax)
Definition: shaper.cpp:53
Vector Y_el_
Output state E-Vector.
int q1d_
Number of quadrature points in 1D.
const Operator * h1_element_restriction_
Vector X_local_
Input state L-vector.
Vector Y_local_
Output state L-vector.
const int ne_
Number of elements in the mesh (rank local)
Structure representing the matrices/tensors needed to evaluate (in reference space) the values...
Definition: fe_base.hpp:135
const Array< int > & GetPrescribedDisplacementTDofs()
Return the T-vector degrees of freedom that have been marked as displaced.
int d1d_
Number of degrees of freedom in 1D.
H1_FECollection h1_fec_
H1 finite element collection.
std::function< void(const int, const Array< double > &, const Array< double > &, const Array< double > &, const Vector &, const Vector &, const Vector &, Vector &)> element_apply_kernel_wrapper
Wrapper for the application of the residual R(U).
void SetPrescribedDisplacement(const Array< int > attr)
Set the attributes which mark the degrees of freedom that have a fixed displacement.
void SetEssentialAttributes(const Array< int > attr)
Set the essential attributes which mark degrees of freedom for the solving process.
Vector data type.
Definition: vector.hpp:60
void GradientMult(const Vector &dX, Vector &Y) const
Multiply the linearization of the residual R(U) wrt to the current state U by a perturbation dX...
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:252
ParFiniteElementSpace h1_fes_
H1 finite element space.
Abstract operator.
Definition: operator.hpp:24
Class for parallel meshes.
Definition: pmesh.hpp:32
Operator & GetGradient(const Vector &U) const override
Get the Gradient object.