MFEM  v4.5.2
Finite element discretization library
neohookean.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_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 
20 using mfem::internal::tensor;
21 using 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  */
39 template <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 / 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  * @return stress
71  */
72  MFEM_HOST_DEVICE static void
74  tensor<double, dim, dim> &dudx,
75  tensor<double, dim, dim> &sigma)
76  {
77  sigma = self->stress(dudx);
78  }
79 
80  /**
81  * @brief Compute the gradient.
82  *
83  * This method is used in the ElasticityDiagonalPreconditioner type to
84  * compute the gradient matrix entries of the current quadrature point,
85  * instead of the action.
86  *
87  * @param[in] dudx
88  * @return
89  */
90  MFEM_HOST_DEVICE tensor<double, dim, dim, dim, dim>
91  gradient(tensor<double, dim, dim> dudx) const
92  {
93  static constexpr auto I = mfem::internal::IsotropicIdentity<dim>();
94 
95  tensor<double, dim, dim> F = I + dudx;
96  tensor<double, dim, dim> invF = inv(F);
97  tensor<double, dim, dim> devB =
98  dev(dudx + transpose(dudx) + dot(dudx, transpose(dudx)));
99  double J = det(F);
100  double coef = (C1 / pow(J, 5.0 / 3.0));
101  return make_tensor<dim, dim, dim, dim>([&](int i, int j, int k,
102  int l)
103  {
104  return 2.0 * (D1 * J * (i == j) - (5.0 / 3.0) * coef * devB[i][j]) *
105  invF[l][k] +
106  2.0 * coef *
107  ((i == k) * F[j][l] + F[i][l] * (j == k) -
108  (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<double, dim, dim>
120  action_of_gradient(const tensor<double, dim, dim> &dudx,
121  const tensor<double, 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<double, dim, dim>{};
148  }
149 
150  MFEM_HOST_DEVICE tensor<double, dim, dim>
151  action_of_gradient_dual(const tensor<double, dim, dim> &dudx,
152  const tensor<double, dim, dim> &ddudx) const
153  {
154  auto sigma = stress(make_tensor<dim, dim>([&](int i, int j)
155  {
156  return mfem::internal::dual<double, double> {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<double, dim, dim>
164  action_of_gradient_enzyme_fwd(const tensor<double, dim, dim> &dudx,
165  const tensor<double, dim, dim> &ddudx) const
166  {
167  tensor<double, dim, dim> sigma{};
168  tensor<double, dim, dim> dsigma{};
169 
170  __enzyme_fwddiff<void>(stress_wrapper, enzyme_const, this, enzyme_dup,
171  &dudx, &ddudx, enzyme_dupnoneed, &sigma, &dsigma);
172  return dsigma;
173  }
174 
175  MFEM_HOST_DEVICE tensor<double, dim, dim>
176  action_of_gradient_enzyme_rev(const tensor<double, dim, dim> &dudx,
177  const tensor<double, dim, dim> &ddudx) const
178  {
179  tensor<double, dim, dim, dim, dim> gradient{};
180  tensor<double, dim, dim> sigma{};
181  tensor<double, 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;
188  __enzyme_autodiff<void>(stress_wrapper, enzyme_const, this, enzyme_dup,
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<double, dim, dim>
199  action_of_gradient_finite_diff(const tensor<double, dim, dim> &dudx,
200  const tensor<double, dim, dim> &ddudx) const
201  {
202  return (stress(dudx + 1.0e-8 * ddudx) - stress(dudx - 1.0e-8 * ddudx)) /
203  2.0e-8;
204  }
205 
206  // d(stress)_{ij} := (d(stress)_ij / d(du_dx)_{kl}) * d(du_dx)_{kl}
207  // Only works with 3D stress
208  MFEM_HOST_DEVICE tensor<double, dim, dim>
209  action_of_gradient_symbolic(const tensor<double, dim, dim> &du_dx,
210  const tensor<double, dim, dim> &ddu_dx) const
211  {
212  static constexpr auto I = mfem::internal::IsotropicIdentity<dim>();
213 
214  tensor<double, dim, dim> F = I + du_dx;
215  tensor<double, dim, dim> invFT = inv(transpose(F));
216  tensor<double, dim, dim> devB =
217  dev(du_dx + transpose(du_dx) + dot(du_dx, transpose(du_dx)));
218  double J = det(F);
219  double coef = (C1 / pow(J, 5.0 / 3.0));
220  double a1 = ddot(invFT, ddu_dx);
221  double a2 = ddot(F, ddu_dx);
222 
223  return (2.0 * D1 * J * a1 - (4.0 / 3.0) * coef * a2) * I -
224  ((10.0 / 3.0) * coef * a1) * devB +
225  (2 * coef) * (dot(ddu_dx, transpose(F)) + dot(F, transpose(ddu_dx)));
226  }
227 
228  // Parameters
229  double D1 = 100.0;
230  double C1 = 50.0;
231 };
232 
233 #endif
Neo-Hookean material.
Definition: neohookean.hpp:40
MFEM_HOST_DEVICE tensor< double, dim, dim > action_of_gradient_dual(const tensor< double, dim, dim > &dudx, const tensor< double, dim, dim > &ddudx) const
Definition: neohookean.hpp:151
MFEM_HOST_DEVICE tensor< double, dim, dim > action_of_gradient_enzyme_fwd(const tensor< double, dim, dim > &dudx, const tensor< double, dim, dim > &ddudx) const
Definition: neohookean.hpp:164
MFEM_HOST_DEVICE tensor< T, dim, dim > stress(const tensor< T, dim, dim > &dudx) const
Compute the stress response.
Definition: neohookean.hpp:52
int enzyme_dupnoneed
MFEM_HOST_DEVICE tensor< double, dim, dim, dim, dim > gradient(tensor< double, dim, dim > dudx) const
Compute the gradient.
Definition: neohookean.hpp:91
Implementation of the tensor class.
int enzyme_dup
MFEM_HOST_DEVICE tensor< double, dim, dim > action_of_gradient_enzyme_rev(const tensor< double, dim, dim > &dudx, const tensor< double, dim, dim > &ddudx) const
Definition: neohookean.hpp:176
MFEM_HOST_DEVICE tensor< double, dim, dim > action_of_gradient(const tensor< double, dim, dim > &dudx, const tensor< double, dim, dim > &ddudx) const
Apply the gradient of the stress.
Definition: neohookean.hpp:120
int enzyme_const
static MFEM_HOST_DEVICE void stress_wrapper(NeoHookeanMaterial< dim, gradient_type > *self, tensor< double, dim, dim > &dudx, tensor< double, dim, dim > &sigma)
A method to wrap the stress calculation into a static function.
Definition: neohookean.hpp:73
int dim
Definition: ex24.cpp:53
MFEM_HOST_DEVICE tensor< double, dim, dim > action_of_gradient_symbolic(const tensor< double, dim, dim > &du_dx, const tensor< double, dim, dim > &ddu_dx) const
Definition: neohookean.hpp:209
MFEM_HOST_DEVICE tensor< double, dim, dim > action_of_gradient_finite_diff(const tensor< double, dim, dim > &dudx, const tensor< double, dim, dim > &ddudx) const
Definition: neohookean.hpp:199
double sigma(const Vector &x)
Definition: maxwell.cpp:102