MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
neohookean.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_MAT_NEOHOOKEAN_HPP
13#define MFEM_ELASTICITY_MAT_NEOHOOKEAN_HPP
14
15#include "general/enzyme.hpp"
16#include "gradient_type.hpp"
17#include "linalg/tensor.hpp"
18#include "mfem.hpp"
19
20using mfem::internal::tensor;
21using mfem::internal::make_tensor;
22
23/**
24 * @brief Neo-Hookean material
25 *
26 * Defines a Neo-Hookean material response. It satisfies the material_type
27 * interface for ElasticityOperator::SetMaterial. This material type allows
28 * choosing the method of derivative calculation in `action_of_gradient`.
29 * Choices include methods derived by hand using symbolic calculation and a
30 * variety of automatically computed gradient applications, like
31 * - Enzyme forward mode
32 * - Enzyme reverse mode
33 * - Dual number type forward mode
34 * - Finite difference mode
35 *
36 * @tparam dim
37 * @tparam gradient_type
38 */
39template <int dim = 3, GradientType gradient_type = GradientType::Symbolic>
41{
42 static_assert(dim == 3, "NeoHookean model currently implemented only in 3D");
43
44 /**
45 * @brief Compute the stress response.
46 *
47 * @param[in] dudx derivative of the displacement
48 * @return
49 */
50 template <typename T>
51 MFEM_HOST_DEVICE tensor<T, dim, dim>
52 stress(const tensor<T, dim, dim> &dudx) const
53 {
54 static constexpr auto I = mfem::internal::IsotropicIdentity<dim>();
55 T J = det(I + dudx);
56 T p = -2.0 * D1 * J * (J - 1);
57 auto devB = dev(dudx + transpose(dudx) + dot(dudx, transpose(dudx)));
58 auto sigma = -(p / J) * I + 2 * (C1 / (T) pow(J, 5.0 / 3.0)) * devB;
59 return sigma;
60 }
61
62 /**
63 * @brief A method to wrap the stress calculation into a static function.
64 *
65 * This is necessary for Enzyme to access the class pointer (self).
66 *
67 * @param[in] self
68 * @param[in] dudx
69 * @param[in] sigma
70 */
71 MFEM_HOST_DEVICE static void
73 tensor<mfem::real_t, dim, dim> &dudx,
74 tensor<mfem::real_t, dim, dim> &sigma)
75 {
76 sigma = self->stress(dudx);
77 }
78
79 /**
80 * @brief Compute the gradient.
81 *
82 * This method is used in the ElasticityDiagonalPreconditioner type to
83 * compute the gradient matrix entries of the current quadrature point,
84 * instead of the action.
85 *
86 * @param[in] dudx
87 * @return
88 */
89 MFEM_HOST_DEVICE tensor<mfem::real_t, dim, dim, dim, dim>
90 gradient(tensor<mfem::real_t, dim, dim> dudx) const
91 {
92 static constexpr auto I = mfem::internal::IsotropicIdentity<dim>();
93
94 tensor<mfem::real_t, dim, dim> F = I + dudx;
95 tensor<mfem::real_t, dim, dim> invF = inv(F);
96 tensor<mfem::real_t, dim, dim> devB =
97 dev(dudx + transpose(dudx) + dot(dudx, transpose(dudx)));
98 mfem::real_t J = det(F);
99 mfem::real_t coef = (C1 / pow(J, 5.0 / 3.0));
100 return make_tensor<dim, dim, dim, dim>([&](int i, int j, int k,
101 int l)
102 {
103 return 2 * (D1 * J * (i == j) -
104 mfem::real_t(5.0 / 3.0) * coef * devB[i][j]) *
105 invF[l][k] +
106 2 * coef *
107 ((i == k) * F[j][l] + F[i][l] * (j == k) -
108 mfem::real_t(2.0 / 3.0) * ((i == j) * F[k][l]));
109 });
110 }
111
112 /**
113 * @brief Apply the gradient of the stress.
114 *
115 * @param[in] dudx
116 * @param[in] ddudx
117 * @return
118 */
119 MFEM_HOST_DEVICE tensor<mfem::real_t, dim, dim>
120 action_of_gradient(const tensor<mfem::real_t, dim, dim> &dudx,
121 const tensor<mfem::real_t, dim, dim> &ddudx) const
122 {
123 if (gradient_type == GradientType::Symbolic)
124 {
125 return action_of_gradient_symbolic(dudx, ddudx);
126 }
127#ifdef MFEM_USE_ENZYME
128 else if (gradient_type == GradientType::EnzymeFwd)
129 {
130 return action_of_gradient_enzyme_fwd(dudx, ddudx);
131 }
132 else if (gradient_type == GradientType::EnzymeRev)
133 {
134 return action_of_gradient_enzyme_rev(dudx, ddudx);
135 }
136#endif
137 else if (gradient_type == GradientType::FiniteDiff)
138 {
139 return action_of_gradient_finite_diff(dudx, ddudx);
140 }
141 else if (gradient_type == GradientType::InternalFwd)
142 {
143 return action_of_gradient_dual(dudx, ddudx);
144 }
145 // Getting to this point is an error.
146 // For now we just return a zero tensor to suppress a warning:
147 return tensor<mfem::real_t, dim, dim> {};
148 }
149
150 MFEM_HOST_DEVICE tensor<mfem::real_t, dim, dim>
151 action_of_gradient_dual(const tensor<mfem::real_t, dim, dim> &dudx,
152 const tensor<mfem::real_t, dim, dim> &ddudx) const
153 {
154 auto sigma = stress(make_tensor<dim, dim>([&](int i, int j)
155 {
156 return mfem::internal::dual<mfem::real_t, mfem::real_t> {dudx[i][j], ddudx[i][j]};
157 }));
158 return make_tensor<dim, dim>(
159 [&](int i, int j) { return sigma[i][j].gradient; });
160 }
161
162#ifdef MFEM_USE_ENZYME
163 MFEM_HOST_DEVICE tensor<mfem::real_t, dim, dim>
164 action_of_gradient_enzyme_fwd(const tensor<mfem::real_t, dim, dim> &dudx,
165 const tensor<mfem::real_t, dim, dim> &ddudx) const
166 {
167 tensor<mfem::real_t, dim, dim> sigma{};
168 tensor<mfem::real_t, dim, dim> dsigma{};
169
171 &dudx, &ddudx, enzyme_dupnoneed, &sigma, &dsigma);
172 return dsigma;
173 }
174
175 MFEM_HOST_DEVICE tensor<mfem::real_t, dim, dim>
176 action_of_gradient_enzyme_rev(const tensor<mfem::real_t, dim, dim> &dudx,
177 const tensor<mfem::real_t, dim, dim> &ddudx) const
178 {
179 tensor<mfem::real_t, dim, dim, dim, dim> gradient{};
180 tensor<mfem::real_t, dim, dim> sigma{};
181 tensor<mfem::real_t, dim, dim> dir{};
182
183 for (int i = 0; i < dim; i++)
184 {
185 for (int j = 0; j < dim; j++)
186 {
187 dir[i][j] = 1;
189 &dudx, &gradient[i][j], enzyme_dupnoneed,
190 &sigma, &dir);
191 dir[i][j] = 0;
192 }
193 }
194 return ddot(gradient, ddudx);
195 }
196#endif
197
198 MFEM_HOST_DEVICE tensor<mfem::real_t, dim, dim>
199 action_of_gradient_finite_diff(const tensor<mfem::real_t, dim, dim> &dudx,
200 const tensor<mfem::real_t, dim, dim> &ddudx) const
201 {
202 return (stress(dudx + mfem::real_t(1.0e-8) * ddudx) -
203 stress(dudx - mfem::real_t(1.0e-8) * ddudx)) /
204 mfem::real_t(2.0e-8);
205 }
206
207 // d(stress)_{ij} := (d(stress)_ij / d(du_dx)_{kl}) * d(du_dx)_{kl}
208 // Only works with 3D stress
209 MFEM_HOST_DEVICE tensor<mfem::real_t, dim, dim>
210 action_of_gradient_symbolic(const tensor<mfem::real_t, dim, dim> &du_dx,
211 const tensor<mfem::real_t, dim, dim> &ddu_dx) const
212 {
213 static constexpr auto I = mfem::internal::IsotropicIdentity<dim>();
214
215 tensor<mfem::real_t, dim, dim> F = I + du_dx;
216 tensor<mfem::real_t, dim, dim> invFT = inv(transpose(F));
217 tensor<mfem::real_t, dim, dim> devB =
218 dev(du_dx + transpose(du_dx) + dot(du_dx, transpose(du_dx)));
219 mfem::real_t J = det(F);
220 mfem::real_t coef = (C1 / pow(J, 5.0 / 3.0));
221 mfem::real_t a1 = ddot(invFT, ddu_dx);
222 mfem::real_t a2 = ddot(F, ddu_dx);
223
224 return ((2 * D1 * J * a1 - mfem::real_t(4.0 / 3.0) * coef * a2) * I -
225 (mfem::real_t(10.0 / 3.0) * coef * a1) * devB +
226 (2 * coef) * (dot(ddu_dx, transpose(F)) + dot(F, transpose(ddu_dx))));
227 }
228
229 // Parameters
232};
233
234#endif
real_t sigma(const Vector &x)
Definition maxwell.cpp:102
int enzyme_dup
int enzyme_dupnoneed
return_type __enzyme_autodiff(Args...)
int enzyme_const
return_type __enzyme_fwddiff(Args...)
int dim
Definition ex24.cpp:53
float real_t
Definition config.hpp:43
real_t p(const Vector &x, real_t t)
Neo-Hookean material.
MFEM_HOST_DEVICE tensor< mfem::real_t, dim, dim > action_of_gradient(const tensor< mfem::real_t, dim, dim > &dudx, const tensor< mfem::real_t, dim, dim > &ddudx) const
Apply the gradient of the stress.
MFEM_HOST_DEVICE tensor< mfem::real_t, dim, dim > action_of_gradient_symbolic(const tensor< mfem::real_t, dim, dim > &du_dx, const tensor< mfem::real_t, dim, dim > &ddu_dx) const
MFEM_HOST_DEVICE tensor< mfem::real_t, dim, dim > action_of_gradient_enzyme_rev(const tensor< mfem::real_t, dim, dim > &dudx, const tensor< mfem::real_t, dim, dim > &ddudx) const
MFEM_HOST_DEVICE tensor< mfem::real_t, dim, dim > action_of_gradient_finite_diff(const tensor< mfem::real_t, dim, dim > &dudx, const tensor< mfem::real_t, dim, dim > &ddudx) const
static MFEM_HOST_DEVICE void stress_wrapper(NeoHookeanMaterial< dim, gradient_type > *self, tensor< mfem::real_t, dim, dim > &dudx, tensor< mfem::real_t, dim, dim > &sigma)
A method to wrap the stress calculation into a static function.
MFEM_HOST_DEVICE tensor< mfem::real_t, dim, dim > action_of_gradient_dual(const tensor< mfem::real_t, dim, dim > &dudx, const tensor< mfem::real_t, dim, dim > &ddudx) const
MFEM_HOST_DEVICE tensor< mfem::real_t, dim, dim > action_of_gradient_enzyme_fwd(const tensor< mfem::real_t, dim, dim > &dudx, const tensor< mfem::real_t, dim, dim > &ddudx) const
MFEM_HOST_DEVICE tensor< T, dim, dim > stress(const tensor< T, dim, dim > &dudx) const
Compute the stress response.
MFEM_HOST_DEVICE tensor< mfem::real_t, dim, dim, dim, dim > gradient(tensor< mfem::real_t, dim, dim > dudx) const
Compute the gradient.
Implementation of the tensor class.