MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
elasticity_operator.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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
16#include "mfem.hpp"
17
18namespace mfem
19{
21{
22public:
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
101 /// Number of degrees of freedom in 1D
102 int d1d_;
103 /// Number of quadrature points in 1D
104 int q1d_;
112 /// Input state L-vector
114 /// Input state E-vector
115 mutable Vector X_el_;
116 /// Output state L-vector
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.
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<real_t> &, const Array<real_t> &,
145 const Array<real_t> &, 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<real_t> &, const Array<real_t> &,
155 const Array<real_t> &, 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<real_t> &, const Array<real_t> &,
165 const Array<real_t> &, 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<real_t> &B_, const Array<real_t> &G_,
188 const Array<real_t> &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<real_t> &B_, const Array<real_t> &G_,
217 const Array<real_t> &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<real_t> &B_, const Array<real_t> &G_,
246 const Array<real_t> &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 */
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 */
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
Structure representing the matrices/tensors needed to evaluate (in reference space) the values,...
Definition fe_base.hpp:137
void SetPrescribedDisplacement(const Array< int > attr)
Set the attributes which mark the degrees of freedom that have a fixed displacement.
const Operator * h1_element_restriction_
int d1d_
Number of degrees of freedom in 1D.
const Array< int > & GetPrescribedDisplacementTDofs()
Return the T-vector degrees of freedom that have been marked as displaced.
Vector Y_local_
Output state L-vector.
Vector X_el_
Input state E-vector.
const int order_
Polynomial order of the FE space.
void SetEssentialAttributes(const Array< int > attr)
Set the essential attributes which mark degrees of freedom for the solving process.
ElasticityOperator(ParMesh &mesh, const int order)
Vector Y_el_
Output state E-Vector.
std::function< void(const int, const Array< real_t > &, const Array< real_t > &, const Array< real_t > &, const Vector &, const Vector &, const Vector &, Vector &)> element_kernel_assemble_diagonal_wrapper
Wrapper for the assembly of the gradient on each diagonal element.
Vector X_local_
Input state L-vector.
const int ne_
Number of elements in the mesh (rank local)
int q1d_
Number of quadrature points in 1D.
std::function< void(const int, const Array< real_t > &, const Array< real_t > &, const Array< real_t > &, const Vector &, const Vector &, const Vector &, Vector &, const Vector &)> element_apply_gradient_kernel_wrapper
Wrapper for the application of the gradient of the residual.
std::function< void(const int, const Array< real_t > &, const Array< real_t > &, const Array< real_t > &, const Vector &, const Vector &, const Vector &, Vector &)> element_apply_kernel_wrapper
Wrapper for the application of the residual R(U).
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.
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...
void SetMaterial(const material_type &material)
Set the material type.
Operator & GetGradient(const Vector &U) const override
Get the Gradient object.
const GeometricFactors * geometric_factors_
ParFiniteElementSpace h1_fes_
H1 finite element space.
H1_FECollection h1_fec_
H1 finite element collection.
void AssembleGradientDiagonal(Vector &Ke_diag, Vector &K_diag_local, Vector &K_diag) const
Assemble the linearization of the residual K = dR(U)/dU.
Structure for storing mesh geometric factors: coordinates, Jacobians, and determinants of the Jacobia...
Definition mesh.hpp:2790
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
Abstract operator.
Definition operator.hpp:25
Abstract parallel finite element space.
Definition pfespace.hpp:29
void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const override
Class for parallel meshes.
Definition pmesh.hpp:34
Vector data type.
Definition vector.hpp:80
int material(Vector &x, Vector &xmin, Vector &xmax)
Definition shaper.cpp:53