PolyFEM
Loading...
Searching...
No Matches
GenericElastic.hpp
Go to the documentation of this file.
1#pragma once
2
5
7
8// non linear NeoHookean material model
9namespace polyfem::assembler
10{
11 enum class AutodiffType
12 {
13 FULL,
14 STRESS,
15 NONE
16 };
17
18 template <typename T>
19 using DefGradMatrix = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3>;
20
21 template <typename Derived>
23 {
24 public:
28
30 virtual ~GenericElastic() = default;
31
32 // energy, gradient, and hessian used in newton method
33 double compute_energy(const NonLinearAssemblerData &data) const override;
34 Eigen::MatrixXd assemble_hessian(const NonLinearAssemblerData &data) const override;
35 Eigen::VectorXd assemble_gradient(const NonLinearAssemblerData &data) const override;
36
37 void assign_stress_tensor(const OutputData &data,
38 const int all_size,
39 const ElasticityTensorType &type,
40 Eigen::MatrixXd &all,
41 const std::function<Eigen::MatrixXd(const Eigen::MatrixXd &)> &fun) const override;
42
44 const Eigen::MatrixXd &mat,
45 Eigen::MatrixXd &stress,
46 Eigen::MatrixXd &result) const override;
47
49 Eigen::MatrixXd &stress,
50 Eigen::MatrixXd &result) const override;
51
53 const Eigen::MatrixXd &vect,
54 Eigen::MatrixXd &stress,
55 Eigen::MatrixXd &result) const override;
56
58 Derived &derived() { return static_cast<Derived &>(*this); }
60 const Derived &derived() const { return static_cast<const Derived &>(*this); }
61
62 // sets material params
63 virtual void add_multimaterial(const int index, const json &params, const Units &units) override = 0;
64
65 bool allow_inversion() const override { return true; }
66
67 virtual bool real_def_grad() const { return true; }
68
69 protected:
71
72 virtual Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> gradient(
73 const RowVectorNd &p,
74 const double t,
75 const int el_id,
76 const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> &F) const
77 {
78 log_and_throw_error("gradient not implemented");
79 }
80
81 virtual Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 9, 9> hessian(
82 const RowVectorNd &p,
83 const double t,
84 const int el_id,
85 const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> &F) const
86 {
87 log_and_throw_error("hessian not implemented");
88 }
89
90 private:
91 // utility function that computes energy, the template is used for double, DScalar1, and DScalar2 in energy, gradient and hessian
92 template <typename T>
94 {
95 typedef Eigen::Matrix<T, Eigen::Dynamic, 1> AutoDiffVect;
96 typedef Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> AutoDiffGradMat;
97 typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> DoubleGradMat;
98
99 AutoDiffVect local_disp;
100 get_local_disp(data, size(), local_disp);
101
102 AutoDiffGradMat def_grad(size(), size());
103
104 T energy = T(0.0);
105
106 const int n_pts = data.da.size();
107 for (long p = 0; p < n_pts; ++p)
108 {
109 compute_disp_grad_at_quad(data, local_disp, p, size(), def_grad);
110
111 // Id + grad d
112 for (int d = 0; d < size(); ++d)
113 def_grad(d, d) += T(1);
114
115 if (!derived().real_def_grad())
116 {
117 DoubleGradMat tmp_jac_it = data.vals.jac_it[p];
118 tmp_jac_it = tmp_jac_it.inverse();
119
120 AutoDiffGradMat jac_it(size(), size());
121 for (long k = 0; k < jac_it.size(); ++k)
122 jac_it(k) = T(tmp_jac_it(k));
123
124 def_grad *= jac_it;
125 }
126
127 const T val = derived().elastic_energy(data.vals.val.row(p), data.t, data.vals.element_id, def_grad);
128
129 energy += val * data.da(p);
130 }
131 return energy;
132 }
133
134 Eigen::MatrixXd assemble_hessian_full_ad(const NonLinearAssemblerData &data) const;
135 Eigen::VectorXd assemble_gradient_full_ad(const NonLinearAssemblerData &data) const;
136
137 Eigen::MatrixXd assemble_hessian_stress_ad(const NonLinearAssemblerData &data) const;
138 Eigen::VectorXd assemble_gradient_stress_ad(const NonLinearAssemblerData &data) const;
139
140 Eigen::MatrixXd assemble_hessian_stress_noad(const NonLinearAssemblerData &data) const;
141 Eigen::VectorXd assemble_gradient_stress_noad(const NonLinearAssemblerData &data) const;
142
143 template <int n_basis, int dim>
144 void compute_gradient_from_stress(const NonLinearAssemblerData &data, Eigen::VectorXd &res) const
145 {
147 typedef Eigen::Matrix<Diff, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> AutoDiffGradMat;
148 typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> DoubleGradMat;
149
150 const int n_pts = data.da.size();
151
152 Eigen::Matrix<double, n_basis, dim> local_disp(data.vals.basis_values.size(), size());
153 local_disp.setZero();
154 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
155 {
156 const auto &bs = data.vals.basis_values[i];
157 for (size_t ii = 0; ii < bs.global.size(); ++ii)
158 {
159 for (int d = 0; d < dim; ++d)
160 {
161 local_disp(i, d) += bs.global[ii].val * data.x(bs.global[ii].index * size() + d);
162 }
163 }
164 }
165
166 Eigen::Matrix<double, dim, dim> def_grad(size(), size());
168 AutoDiffGradMat def_grad_ad(size(), size());
169
170 res.resize(data.vals.basis_values.size() * size());
171 res.setZero();
172
173 for (long p = 0; p < n_pts; ++p)
174 {
175 Eigen::Matrix<double, n_basis, dim> grad(data.vals.basis_values.size(), size());
176
177 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
178 grad.row(i) = data.vals.basis_values[i].grad.row(p);
179
180 Eigen::Matrix<double, dim, dim> jac_it = data.vals.jac_it[p];
181 Eigen::Matrix<double, n_basis, dim> G = grad * jac_it;
182 def_grad = local_disp.transpose() * G + Eigen::Matrix<double, dim, dim>::Identity(size(), size());
183
184 for (int d1 = 0; d1 < size(); ++d1)
185 for (int d2 = 0; d2 < size(); ++d2)
186 def_grad_ad(d1, d2) = Diff(d1 * size() + d2, def_grad(d1, d2));
187
188 if (!derived().real_def_grad())
189 {
190 DoubleGradMat tmp_jac_it = data.vals.jac_it[p];
191 tmp_jac_it = tmp_jac_it.inverse();
192
193 AutoDiffGradMat jac_it(size(), size());
194 for (long k = 0; k < jac_it.size(); ++k)
195 jac_it(k) = Diff(tmp_jac_it(k));
196
197 def_grad_ad *= jac_it;
198 }
199
200 const Diff val = derived().elastic_energy(data.vals.val.row(p), data.t, data.vals.element_id, def_grad_ad);
201
202 Eigen::Matrix<double, dim, dim> P;
203 for (int i = 0; i < size(); ++i)
204 for (int j = 0; j < size(); ++j)
205 P(i, j) = val.getGradient()(i * size() + j);
206
207 const Eigen::Matrix<double, n_basis, dim> Rloc = G * P.transpose() * data.da(p);
208
209 for (int a = 0; a < data.vals.basis_values.size(); ++a)
210 for (int d = 0; d < size(); ++d)
211 res(a * size() + d) += Rloc(a, d);
212 }
213 }
214
215 template <int n_basis, int dim>
216 void compute_gradient_from_stress_noad(const NonLinearAssemblerData &data, Eigen::VectorXd &res) const
217 {
218 typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> DoubleGradMat;
219
220 const int n_pts = data.da.size();
221
222 Eigen::Matrix<double, n_basis, dim> local_disp(data.vals.basis_values.size(), size());
223 local_disp.setZero();
224 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
225 {
226 const auto &bs = data.vals.basis_values[i];
227 for (size_t ii = 0; ii < bs.global.size(); ++ii)
228 {
229 for (int d = 0; d < dim; ++d)
230 {
231 local_disp(i, d) += bs.global[ii].val * data.x(bs.global[ii].index * size() + d);
232 }
233 }
234 }
235
236 Eigen::Matrix<double, dim, dim> def_grad(size(), size());
237 res.resize(data.vals.basis_values.size() * size());
238 res.setZero();
239
240 for (long p = 0; p < n_pts; ++p)
241 {
242 Eigen::Matrix<double, n_basis, dim> grad(data.vals.basis_values.size(), size());
243
244 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
245 grad.row(i) = data.vals.basis_values[i].grad.row(p);
246
247 Eigen::Matrix<double, dim, dim> jac_it = data.vals.jac_it[p];
248 Eigen::Matrix<double, n_basis, dim> G = grad * jac_it;
249 def_grad = local_disp.transpose() * G + Eigen::Matrix<double, dim, dim>::Identity(size(), size());
250
251 if (!derived().real_def_grad())
252 {
253 DoubleGradMat jac_it = data.vals.jac_it[p];
254 jac_it = jac_it.inverse();
255
256 def_grad *= jac_it;
257 }
258
259 const Eigen::Matrix<double, dim, dim> P = derived().gradient(data.vals.val.row(p), data.t, data.vals.element_id, def_grad);
260 const Eigen::Matrix<double, n_basis, dim> Bgrad = derived().real_def_grad() ? G : grad;
261 const Eigen::Matrix<double, n_basis, dim> Rloc = Bgrad * P.transpose() * data.da(p);
262
263 for (int a = 0; a < data.vals.basis_values.size(); ++a)
264 for (int d = 0; d < size(); ++d)
265 res(a * size() + d) += Rloc(a, d);
266 }
267 }
268
269 template <int n_basis, int dim>
270 void compute_hessian_from_stress(const NonLinearAssemblerData &data, Eigen::MatrixXd &H) const
271 {
272 typedef DScalar2<double, Eigen::Matrix<double, Eigen::Dynamic, 1, 0, 9, 1>, Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 9, 9>> Diff2;
273 typedef Eigen::Matrix<Diff2, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> AutoDiffGradMat;
274 typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> DoubleGradMat;
275
276 const int d = size();
277 const int nb = data.vals.basis_values.size();
278 const int n_pts = data.da.size();
279
280 Eigen::Matrix<double, n_basis, dim> local_disp(nb, d);
281 local_disp.setZero();
282 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
283 {
284 const auto &bs = data.vals.basis_values[i];
285 for (size_t ii = 0; ii < bs.global.size(); ++ii)
286 {
287 for (int dd = 0; dd < d; ++dd)
288 local_disp(i, dd) += bs.global[ii].val * data.x(bs.global[ii].index * d + dd);
289 }
290 }
291
293 AutoDiffGradMat def_grad_ad(d, d);
294 Eigen::Matrix<double, dim, dim> def_grad(d, d);
295
296 H.resize(nb * d, nb * d);
297 H.setZero();
298
299 for (long p = 0; p < n_pts; ++p)
300 {
301 Eigen::Matrix<double, n_basis, dim> grad_ref(nb, d);
302 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
303 grad_ref.row(i) = data.vals.basis_values[i].grad.row(p);
304
305 const Eigen::Matrix<double, dim, dim> jac_it = data.vals.jac_it[p];
306 const Eigen::Matrix<double, n_basis, dim> G = grad_ref * jac_it;
307
308 def_grad = local_disp.transpose() * G + Eigen::Matrix<double, dim, dim>::Identity();
309
310 for (int i = 0; i < d; ++i)
311 for (int j = 0; j < d; ++j)
312 def_grad_ad(i, j) = Diff2(i * d + j, def_grad(i, j));
313
314 if (!derived().real_def_grad())
315 {
316 DoubleGradMat tmp_jac_it = data.vals.jac_it[p];
317 tmp_jac_it = tmp_jac_it.inverse();
318
319 AutoDiffGradMat jac_it(size(), size());
320 for (long k = 0; k < jac_it.size(); ++k)
321 jac_it(k) = Diff2(tmp_jac_it(k));
322
323 def_grad_ad *= jac_it;
324 }
325
326 const Diff2 val = derived().elastic_energy(data.vals.val.row(p), data.t, data.vals.element_id, def_grad_ad);
327
328 Eigen::Matrix<double, dim * dim, dim * dim> A(d * d, d * d);
329 A.setZero();
330
331 // Material tangent A = d vec(P) / d vec(F)
332 // Since P = dW/dF, A is just the Hessian of W wrt vec(F).
333 for (int r = 0; r < d * d; ++r)
334 for (int c = 0; c < d * d; ++c)
335 A(r, c) = val.getHessian()(r, c);
336
337 for (int a = 0; a < nb; ++a)
338 {
339 const Eigen::Matrix<double, dim * dim, dim> Ba = compute_B_block<dim>(G.row(a));
340
341 for (int b = 0; b < nb; ++b)
342 {
343 const Eigen::Matrix<double, dim * dim, dim> Bb = compute_B_block<dim>(G.row(b));
344
345 const Eigen::Matrix<double, dim, dim> Kab =
346 Ba.transpose() * A * Bb * data.da(p);
347
348 H.template block<dim, dim>(a * d, b * d) += Kab;
349 }
350 }
351 }
352 }
353
354 template <int n_basis, int dim>
355 void compute_hessian_from_stress_noad(const NonLinearAssemblerData &data, Eigen::MatrixXd &H) const
356 {
357 typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> DoubleGradMat;
358
359 const int d = size();
360 const int nb = data.vals.basis_values.size();
361 const int n_pts = data.da.size();
362
363 Eigen::Matrix<double, n_basis, dim> local_disp(nb, d);
364 local_disp.setZero();
365 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
366 {
367 const auto &bs = data.vals.basis_values[i];
368 for (size_t ii = 0; ii < bs.global.size(); ++ii)
369 {
370 for (int dd = 0; dd < d; ++dd)
371 local_disp(i, dd) += bs.global[ii].val * data.x(bs.global[ii].index * d + dd);
372 }
373 }
374
375 Eigen::Matrix<double, dim, dim> def_grad(d, d);
376
377 H.resize(nb * d, nb * d);
378 H.setZero();
379
380 for (long p = 0; p < n_pts; ++p)
381 {
382 Eigen::Matrix<double, n_basis, dim> grad_ref(nb, d);
383 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
384 grad_ref.row(i) = data.vals.basis_values[i].grad.row(p);
385
386 const Eigen::Matrix<double, dim, dim> jac_it = data.vals.jac_it[p];
387 const Eigen::Matrix<double, n_basis, dim> G = grad_ref * jac_it;
388
389 def_grad = local_disp.transpose() * G + Eigen::Matrix<double, dim, dim>::Identity();
390
391 if (!derived().real_def_grad())
392 {
393 DoubleGradMat jac_it = data.vals.jac_it[p];
394 jac_it = jac_it.inverse();
395
396 def_grad *= jac_it;
397 }
398
399 // Material tangent A = d vec(P) / d vec(F)
400 // Since P = dW/dF, A is just the Hessian of W wrt vec(F).
401 Eigen::Matrix<double, dim * dim, dim * dim> A = derived().hessian(data.vals.val.row(p), data.t, data.vals.element_id, def_grad);
402
403 const Eigen::Matrix<double, n_basis, dim> Bgrad = derived().real_def_grad() ? G : grad_ref;
404
405 for (int a = 0; a < nb; ++a)
406 {
407 const Eigen::Matrix<double, dim * dim, dim> Ba = compute_B_block<dim>(Bgrad.row(a));
408
409 for (int b = 0; b < nb; ++b)
410 {
411 const Eigen::Matrix<double, dim * dim, dim> Bb = compute_B_block<dim>(Bgrad.row(b));
412
413 const Eigen::Matrix<double, dim, dim> Kab =
414 Ba.transpose() * A * Bb * data.da(p);
415
416 H.template block<dim, dim>(a * d, b * d) += Kab;
417 }
418 }
419 }
420 }
421
422 template <int dim>
423 Eigen::Matrix<double, dim * dim, dim> compute_B_block(const Eigen::Matrix<double, 1, dim> &g) const;
424 };
425} // namespace polyfem::assembler
double val
Definition Assembler.cpp:86
virtual void assemble_gradient(const bool is_volume, const int n_basis, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const AssemblyValsCache &cache, const double t, const double dt, const Eigen::MatrixXd &displacement, const Eigen::MatrixXd &displacement_prev, Eigen::MatrixXd &rhs) const
virtual void assemble_hessian(const bool is_volume, const int n_basis, const bool project_to_psd, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const AssemblyValsCache &cache, const double t, const double dt, const Eigen::MatrixXd &displacement, const Eigen::MatrixXd &displacement_prev, utils::MatrixCache &mat_cache, StiffnessMatrix &grad) const
virtual double assemble_energy(const bool is_volume, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const AssemblyValsCache &cache, const double t, const double dt, const Eigen::MatrixXd &displacement, const Eigen::MatrixXd &displacement_prev) const
Definition Assembler.hpp:79
std::vector< Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > > jac_it
void compute_stress_grad_multiply_stress(const OptAssemblerData &data, Eigen::MatrixXd &stress, Eigen::MatrixXd &result) const override
void compute_hessian_from_stress_noad(const NonLinearAssemblerData &data, Eigen::MatrixXd &H) const
void assign_stress_tensor(const OutputData &data, const int all_size, const ElasticityTensorType &type, Eigen::MatrixXd &all, const std::function< Eigen::MatrixXd(const Eigen::MatrixXd &)> &fun) const override
Eigen::Matrix< double, dim *dim, dim > compute_B_block(const Eigen::Matrix< double, 1, dim > &g) const
Eigen::VectorXd assemble_gradient_full_ad(const NonLinearAssemblerData &data) const
virtual Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, 0, 9, 9 > hessian(const RowVectorNd &p, const double t, const int el_id, const Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > &F) const
Derived & derived()
Returns this as a reference to derived class.
Eigen::MatrixXd assemble_hessian_full_ad(const NonLinearAssemblerData &data) const
const Derived & derived() const
Returns this as a const reference to derived class.
Eigen::MatrixXd assemble_hessian_stress_noad(const NonLinearAssemblerData &data) const
Eigen::VectorXd assemble_gradient_stress_noad(const NonLinearAssemblerData &data) const
Eigen::MatrixXd assemble_hessian_stress_ad(const NonLinearAssemblerData &data) const
Eigen::MatrixXd assemble_hessian(const NonLinearAssemblerData &data) const override
virtual void add_multimaterial(const int index, const json &params, const Units &units) override=0
void compute_stress_grad_multiply_vect(const OptAssemblerData &data, const Eigen::MatrixXd &vect, Eigen::MatrixXd &stress, Eigen::MatrixXd &result) const override
Eigen::VectorXd assemble_gradient(const NonLinearAssemblerData &data) const override
bool allow_inversion() const override
Eigen::VectorXd assemble_gradient_stress_ad(const NonLinearAssemblerData &data) const
void compute_hessian_from_stress(const NonLinearAssemblerData &data, Eigen::MatrixXd &H) const
void compute_gradient_from_stress_noad(const NonLinearAssemblerData &data, Eigen::VectorXd &res) const
virtual Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > gradient(const RowVectorNd &p, const double t, const int el_id, const Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > &F) const
void compute_gradient_from_stress(const NonLinearAssemblerData &data, Eigen::VectorXd &res) const
double compute_energy(const NonLinearAssemblerData &data) const override
T compute_energy_aux(const NonLinearAssemblerData &data) const
void compute_stress_grad_multiply_mat(const OptAssemblerData &data, const Eigen::MatrixXd &mat, Eigen::MatrixXd &stress, Eigen::MatrixXd &result) const override
Used for test only.
Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > DefGradMatrix
Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > AutoDiffGradMat
Eigen::Matrix< T, Eigen::Dynamic, 1 > AutoDiffVect
DScalar1< double, Eigen::Matrix< double, Eigen::Dynamic, 1 > > Diff
void compute_disp_grad_at_quad(const assembler::NonLinearAssemblerData &data, const AutoDiffVect &local_disp, const int p, const int size, AutoDiffGradMat &def_grad)
nlohmann::json json
Definition Common.hpp:9
void get_local_disp(const assembler::NonLinearAssemblerData &data, const int size, AutoDiffVect &local_disp)
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd
Definition Types.hpp:13
void log_and_throw_error(const std::string &msg)
Definition Logger.cpp:73
Automatic differentiation scalar with first-order derivatives.
Definition autodiff.h:112
Automatic differentiation scalar with first- and second-order derivatives.
Definition autodiff.h:493
static void setVariableCount(size_t value)
Set the independent variable count used by the automatic differentiation layer.
Definition autodiff.h:54