MFEM v4.9.0
Finite element discretization library
Loading...
Searching...
No Matches
invariants.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, 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_INVARIANTS_HPP
13#define MFEM_INVARIANTS_HPP
14
15#include "../general/error.hpp"
16#include <cmath>
17
18namespace mfem
19{
20
21// Matrix invariants and their derivatives for 2x2 and 3x3 matrices.
22
23/** @brief Auxiliary class used as the default for the second template parameter
24 in the classes InvariantsEvaluator2D and InvariantsEvaluator3D. */
25template <typename scalar_t>
26struct ScalarOps
27{
28 static scalar_t sign(const scalar_t &a)
29 { return (a >= scalar_t(0)) ? scalar_t(1) : scalar_t(-1); }
30
31 static scalar_t pow(const scalar_t &x, int m, int n)
32 { return std::pow(x, scalar_t(m)/n); }
33};
34
35
36/** @brief Auxiliary class for evaluating the 2x2 matrix invariants and their
37 first and second derivatives. */
38/**
39 The type `scalar_t` must support the standard operations:
40
41 =, +=, -=, +, -, *, /, unary -, int*scalar_t, int/scalar_t, scalar_t/int
42
43 The type `scalar_ops` must define the static method:
44
45 scalar_t sign(const scalar_t &);
46*/
47template <typename scalar_t, typename scalar_ops = ScalarOps<scalar_t> >
48class InvariantsEvaluator2D
49{
50protected:
51 // Transformation Jacobian
52 const scalar_t *J;
53
54 // Invariants: I_1 = ||J||_F^2, \bar{I}_1 = I_1/det(J), \bar{I}_2 = det(J).
55 scalar_t I1, I1b, I2b;
56
57 // Derivatives of I1, I1b, I2, and I2b using column-major storage.
58 scalar_t dI1[4], dI1b[4], dI2[4], dI2b[4];
59
60 int D_height, alloc_height;
61 const scalar_t *D; // Always points to external data or is empty
62 scalar_t *DaJ, *DJt, *DXt, *DYt;
63
64 enum EvalMasks
65 {
66 HAVE_I1 = 1,
67 HAVE_I1b = 2,
68 HAVE_I2b = 4,
69 HAVE_dI1 = 8,
70 HAVE_dI1b = 16,
71 HAVE_dI2 = 32,
72 HAVE_dI2b = 64,
73 HAVE_DaJ = 128, // D adj(J) = D dI2b^t
74 HAVE_DJt = 256 // D J^t
75 };
76
77 // Bitwise OR of EvalMasks
78 int eval_state;
79
80 bool dont(int have_mask) const { return !(eval_state & have_mask); }
81
82 void Eval_I1()
83 {
84 eval_state |= HAVE_I1;
85 I1 = J[0]*J[0] + J[1]*J[1] + J[2]*J[2] + J[3]*J[3];
86 }
87 void Eval_I1b()
88 {
89 eval_state |= HAVE_I1b;
90 I1b = Get_I1()/Get_I2b();
91 }
92 void Eval_I2b()
93 {
94 eval_state |= HAVE_I2b;
95 const scalar_t det = J[0]*J[3] - J[1]*J[2];
96 I2b = det;
97 }
98 void Eval_dI1()
99 {
100 eval_state |= HAVE_dI1;
101 dI1[0] = 2*J[0]; dI1[2] = 2*J[2];
102 dI1[1] = 2*J[1]; dI1[3] = 2*J[3];
103 }
104 void Eval_dI1b()
105 {
106 eval_state |= HAVE_dI1b;
107 // I1b = I1/I2b
108 // dI1b = (1/I2b)*dI1 - (I1/I2b^2)*dI2b = (2/I2b)*[J - (I1b/2)*dI2b]
109 const scalar_t c1 = 2/Get_I2b();
110 const scalar_t c2 = Get_I1b()/2;
111 Get_dI2b();
112 dI1b[0] = c1*(J[0] - c2*dI2b[0]);
113 dI1b[1] = c1*(J[1] - c2*dI2b[1]);
114 dI1b[2] = c1*(J[2] - c2*dI2b[2]);
115 dI1b[3] = c1*(J[3] - c2*dI2b[3]);
116 }
117 void Eval_dI2()
118 {
119 eval_state |= HAVE_dI2;
120 // I2 = I2b^2
121 // dI2 = 2*I2b*dI2b = 2*det(J)*adj(J)^T
122 const scalar_t c1 = 2*Get_I2b();
123 Get_dI2b();
124 dI2[0] = c1*dI2b[0];
125 dI2[1] = c1*dI2b[1];
126 dI2[2] = c1*dI2b[2];
127 dI2[3] = c1*dI2b[3];
128 }
129 void Eval_dI2b()
130 {
131 eval_state |= HAVE_dI2b;
132 // I2b = det(J)
133 // dI2b = adj(J)^T
134 Get_I2b();
135 dI2b[0] = J[3];
136 dI2b[1] = -J[2];
137 dI2b[2] = -J[1];
138 dI2b[3] = J[0];
139 }
140 void Eval_DaJ() // D adj(J) = D dI2b^t
141 {
142 eval_state |= HAVE_DaJ;
143 Get_dI2b();
144 Eval_DZt(dI2b, &DaJ);
145 }
146 void Eval_DJt() // D J^t
147 {
148 eval_state |= HAVE_DJt;
149 Eval_DZt(J, &DJt);
150 }
151 void Eval_DZt(const scalar_t *Z, scalar_t **DZt_ptr)
152 {
153 MFEM_ASSERT(D != NULL, "");
154 const int nd = D_height;
155 scalar_t *DZt = *DZt_ptr;
156 if (DZt == NULL) { *DZt_ptr = DZt = new scalar_t[2*alloc_height]; }
157 for (int i = 0; i < nd; i++)
158 {
159 const int i0 = i+nd*0, i1 = i+nd*1;
160 DZt[i0] = D[i0]*Z[0] + D[i1]*Z[2];
161 DZt[i1] = D[i0]*Z[1] + D[i1]*Z[3];
162 }
163 }
164
165public:
166 /// The Jacobian should use column-major storage.
167 InvariantsEvaluator2D(const scalar_t *Jac = NULL)
168 : J(Jac), D_height(), alloc_height(), D(), DaJ(), DJt(), DXt(), DYt(),
169 eval_state(0) { }
170
171 ~InvariantsEvaluator2D()
172 {
173 delete [] DYt;
174 delete [] DXt;
175 delete [] DJt;
176 delete [] DaJ;
177 }
178
179 /// The Jacobian should use column-major storage.
180 void SetJacobian(const scalar_t *Jac) { J = Jac; eval_state = 0; }
181
182 /// The @a Deriv matrix is `dof x 2`, using column-major storage.
183 void SetDerivativeMatrix(int height, const scalar_t *Deriv)
184 {
185 eval_state &= ~(HAVE_DaJ | HAVE_DJt);
186 if (alloc_height < height)
187 {
188 delete [] DYt; DYt = NULL;
189 delete [] DXt; DXt = NULL;
190 delete [] DJt; DJt = NULL;
191 delete [] DaJ; DaJ = NULL;
192 alloc_height = height;
193 }
194 D_height = height;
195 D = Deriv;
196 }
197
198 scalar_t Get_I1() { if (dont(HAVE_I1 )) { Eval_I1(); } return I1; }
199 scalar_t Get_I1b() { if (dont(HAVE_I1b)) { Eval_I1b(); } return I1b; }
200 scalar_t Get_I2() { if (dont(HAVE_I2b)) { Eval_I2b(); } return I2b*I2b; }
201 scalar_t Get_I2b() { if (dont(HAVE_I2b)) { Eval_I2b(); } return I2b; }
202
203 const scalar_t *Get_dI1()
204 {
205 if (dont(HAVE_dI1 )) { Eval_dI1(); } return dI1;
206 }
207 const scalar_t *Get_dI1b()
208 {
209 if (dont(HAVE_dI1b)) { Eval_dI1b(); } return dI1b;
210 }
211 const scalar_t *Get_dI2()
212 {
213 if (dont(HAVE_dI2)) { Eval_dI2(); } return dI2;
214 }
215 const scalar_t *Get_dI2b()
216 {
217 if (dont(HAVE_dI2b)) { Eval_dI2b(); } return dI2b;
218 }
219
220 // Assemble operation for tensor X with components X_jslt:
221 // A(i+nd*j,k+nd*l) += (\sum_st w D_is X_jslt D_kt)
222 // 0 <= i,k < nd, 0 <= j,l,s,t < 2
223 // where nd is the height of D, i.e. the number of DOFs in one component.
224
225 void Assemble_ddI1(scalar_t w, scalar_t *A)
226 {
227 // ddI1_jslt = 2 I_jslt = 2 δ_jl δ_st
228 // A(i+nd*j,k+nd*l) += (\sum_st 2 w D_is δ_jl δ_st D_kt)
229 // or
230 // A(i+nd*j,k+nd*l) += (2 w) (\sum_s D_is D_ks) δ_jl
231 // A(i+nd*j,k+nd*l) += (2 w) (D D^t)_ik δ_jl
232
233 const int nd = D_height;
234 const int ah = 2*nd;
235 const scalar_t a = 2*w;
236 for (int i = 0; i < nd; i++)
237 {
238 const int i0 = i+nd*0, i1 = i+nd*1;
239 const scalar_t aDi[2] = { a*D[i0], a*D[i1] };
240 // k == i
241 const scalar_t aDDt_ii = aDi[0]*D[i0] + aDi[1]*D[i1];
242 A[i0+ah*i0] += aDDt_ii;
243 A[i1+ah*i1] += aDDt_ii;
244 // 0 <= k < i
245 for (int k = 0; k < i; k++)
246 {
247 const int k0 = k+nd*0, k1 = k+nd*1;
248 const scalar_t aDDt_ik = aDi[0]*D[k0] + aDi[1]*D[k1];
249 A[i0+ah*k0] += aDDt_ik;
250 A[k0+ah*i0] += aDDt_ik;
251 A[i1+ah*k1] += aDDt_ik;
252 A[k1+ah*i1] += aDDt_ik;
253 }
254 }
255 }
256 void Assemble_ddI1b(scalar_t w, scalar_t *A)
257 {
258 // ddI1b = X1 + X2 + X3, where
259 // X1_ijkl = (I1b/I2) [ (δ_ks δ_it + δ_kt δ_si) dI2b_tj dI2b_sl ]
260 // = (I1b/I2) [ dI2b_ij dI2b_kl + dI2b_kj dI2b_il ]
261 // X2_ijkl = (2/I2b) δ_ik δ_jl = (1/I2b) ddI1_ijkl
262 // X3_ijkl = -(2/I2) (δ_ks δ_it) (J_tj dI2b_sl + dI2b_tj J_sl)
263 // = -(2/I2) (J_ij dI2b_kl + dI2b_ij J_kl)
264 //
265 // A(i+nd*j,k+nd*l) += (\sum_st w D_is ddI1b_jslt D_kt)
266 // or
267 // A(i+nd*j,k+nd*l) +=
268 // w (I1b/I2) [(D dI2b^t)_ij (D dI2b^t)_kl +
269 // (D dI2b^t)_il (D dI2b^t)_kj]
270 // + w (2/I2b) δ_jl (D D^t)_ik
271 // - w (2/I2) [(D J^t)_ij (D dI2b^t)_kl + (D dI2b^t)_ij (D J^t)_kl]
272
273 if (dont(HAVE_DaJ)) { Eval_DaJ(); }
274 if (dont(HAVE_DJt)) { Eval_DJt(); }
275 const int nd = D_height;
276 const int ah = 2*nd;
277 const scalar_t a = w*Get_I1b()/Get_I2();
278 const scalar_t b = 2*w/Get_I2b();
279 const scalar_t c = -2*w/Get_I2();
280 for (int i = 0; i < nd; i++)
281 {
282 const int i0 = i+nd*0, i1 = i+nd*1;
283 const scalar_t aDaJ_i[2] = { a*DaJ[i0], a*DaJ[i1] };
284 const scalar_t bD_i[2] = { b*D[i0], b*D[i1] };
285 const scalar_t cDJt_i[2] = { c*DJt[i0], c*DJt[i1] };
286 const scalar_t cDaJ_i[2] = { c*DaJ[i0], c*DaJ[i1] };
287 // k == i
288 {
289 // Symmetries: A2_ii_00 = A2_ii_11
290 const scalar_t A2_ii = bD_i[0]*D[i0] + bD_i[1]*D[i1];
291
292 A[i0+ah*i0] += 2*(aDaJ_i[0] + cDJt_i[0])*DaJ[i0] + A2_ii;
293
294 // Symmetries: A_ii_01 = A_ii_10
295 const scalar_t A_ii_01 =
296 (2*aDaJ_i[0] + cDJt_i[0])*DaJ[i1] + cDaJ_i[0]*DJt[i1];
297 A[i0+ah*i1] += A_ii_01;
298 A[i1+ah*i0] += A_ii_01;
299
300 A[i1+ah*i1] += 2*(aDaJ_i[1] + cDJt_i[1])*DaJ[i1] + A2_ii;
301 }
302 // 0 <= k < i
303 for (int k = 0; k < i; k++)
304 {
305 const int k0 = k+nd*0, k1 = k+nd*1;
306 // Symmetries: A1_ik_01 = A1_ik_10 = A1_ki_01 = A1_ki_10
307 const scalar_t A1_ik_01 = aDaJ_i[0]*DaJ[k1] + aDaJ_i[1]*DaJ[k0];
308
309 // Symmetries: A2_ik_00 = A2_ik_11 = A2_ki_00 = A2_ki_11
310 const scalar_t A2_ik = bD_i[0]*D[k0] + bD_i[1]*D[k1];
311
312 const scalar_t A_ik_00 =
313 (2*aDaJ_i[0] + cDJt_i[0])*DaJ[k0] + A2_ik + cDaJ_i[0]*DJt[k0];
314 A[i0+ah*k0] += A_ik_00;
315 A[k0+ah*i0] += A_ik_00;
316
317 const scalar_t A_ik_01 =
318 A1_ik_01 + cDJt_i[0]*DaJ[k1] + cDaJ_i[0]*DJt[k1];
319 A[i0+ah*k1] += A_ik_01;
320 A[k1+ah*i0] += A_ik_01;
321
322 const scalar_t A_ik_10 =
323 A1_ik_01 + cDJt_i[1]*DaJ[k0] + cDaJ_i[1]*DJt[k0];
324 A[i1+ah*k0] += A_ik_10;
325 A[k0+ah*i1] += A_ik_10;
326
327 const scalar_t A_ik_11 =
328 (2*aDaJ_i[1] + cDJt_i[1])*DaJ[k1] + A2_ik + cDaJ_i[1]*DJt[k1];
329 A[i1+ah*k1] += A_ik_11;
330 A[k1+ah*i1] += A_ik_11;
331 }
332 }
333 }
334 void Assemble_ddI2(scalar_t w, scalar_t *A)
335 {
336 // ddI2_ijkl = 2 (2 δ_ks δ_it - δ_kt δ_si) dI2b_tj dI2b_sl
337 // = 4 dI2b_ij dI2b_kl - 2 dI2b_kj dI2b_il
338 // = 2 dI2b_ij dI2b_kl + 2 (dI2b_ij dI2b_kl - dI2b_kj dI2b_il)
339 //
340 // A(i+nd*j,k+nd*l) += (\sum_st w D_is ddI2_jslt D_kt)
341 // or
342 // A(i+nd*j,k+nd*l) +=
343 // (\sum_st w D_is (4 dI2b_js dI2b_lt - 2 dI2b_ls dI2b_jt) D_kt)
344 // A(i+nd*j,k+nd*l) +=
345 // 2 w [2 (D dI2b^t)_ij (D dI2b^t)_kl - (D dI2b^t)_il (D dI2b^t)_kj]
346 //
347 // Note: the expression
348 // (D dI2b^t)_ij (D dI2b^t)_kl - (D dI2b^t)_il (D dI2b^t)_kj
349 // is the determinant of the 2x2 matrix formed by rows {i,k} and columns
350 // {j,l} from the matrix (D dI2b^t).
351
352 if (dont(HAVE_DaJ)) { Eval_DaJ(); }
353 const int nd = D_height;
354 const int ah = 2*nd;
355 const scalar_t a = 2*w;
356 for (int i = 0; i < ah; i++)
357 {
358 const scalar_t avi = a*DaJ[i];
359 A[i+ah*i] += avi*DaJ[i];
360 for (int j = 0; j < i; j++)
361 {
362 const scalar_t aVVt_ij = avi*DaJ[j];
363 A[i+ah*j] += aVVt_ij;
364 A[j+ah*i] += aVVt_ij;
365 }
366 }
367 const int j = 1, l = 0;
368 for (int i = 0; i < nd; i++)
369 {
370 const int ij = i+nd*j, il = i+nd*l;
371 const scalar_t aDaJ_ij = a*DaJ[ij], aDaJ_il = a*DaJ[il];
372 for (int k = 0; k < i; k++)
373 {
374 const int kj = k+nd*j, kl = k+nd*l;
375 const scalar_t A_ijkl = aDaJ_ij*DaJ[kl] - aDaJ_il*DaJ[kj];
376 A[ij+ah*kl] += A_ijkl;
377 A[kl+ah*ij] += A_ijkl;
378 A[kj+ah*il] -= A_ijkl;
379 A[il+ah*kj] -= A_ijkl;
380 }
381 }
382 }
383 void Assemble_ddI2b(scalar_t w, scalar_t *A)
384 {
385 // ddI2b_ijkl = (1/I2b) (δ_ks δ_it - δ_kt δ_si) dI2b_tj dI2b_sl
386 // [j -> u], [l -> v], [i -> j], [k -> l]
387 // ddI2b_julv = (1/I2b) (δ_ls δ_jt - δ_lt δ_sj) dI2b_tu dI2b_sv
388 //
389 // A(i+nd*j,k+nd*l) += (\sum_st w D_is ddI2b_jslt D_kt)
390 // or
391 // A(i+nd*j,k+nd*l) += (\sum_uv w D_iu ddI2b_julv D_kv)
392 // A(i+nd*j,k+nd*l) +=
393 // (\sum_uvst (w/I2b)
394 // D_iu (δ_ls δ_jt - δ_lt δ_sj) dI2b_tu dI2b_sv D_kv)
395 // A(i+nd*j,k+nd*l) +=
396 // (\sum_st (w/I2b)
397 // (D dI2b^t)_it (δ_ls δ_jt - δ_lt δ_sj) (D dI2b^t)_ks)
398 // A(i+nd*j,k+nd*l) += (w/I2b)
399 // [ (D dI2b^t)_ij (D dI2b^t)_kl - (D dI2b^t)_il (D dI2b^t)_kj ]
400
401 if (dont(HAVE_DaJ)) { Eval_DaJ(); }
402 const int nd = D_height;
403 const int ah = 2*nd;
404 const int j = 1, l = 0;
405 const scalar_t a = w/Get_I2b();
406 for (int i = 0; i < nd; i++)
407 {
408 const int ij = i+nd*j, il = i+nd*l;
409 const scalar_t aDaJ_ij = a*DaJ[ij], aDaJ_il = a*DaJ[il];
410 for (int k = 0; k < i; k++)
411 {
412 const int kj = k+nd*j, kl = k+nd*l;
413 const scalar_t A_ijkl = aDaJ_ij*DaJ[kl] - aDaJ_il*DaJ[kj];
414 A[ij+ah*kl] += A_ijkl;
415 A[kl+ah*ij] += A_ijkl;
416 A[kj+ah*il] -= A_ijkl;
417 A[il+ah*kj] -= A_ijkl;
418 }
419 }
420 }
421 // Assemble the contribution from the term: T_ijkl = X_ij Y_kl + Y_ij X_kl,
422 // where X and Y are pointers to 2x2 matrices stored in column-major layout.
423 //
424 // The contribution to the matrix A is given by:
425 // A(i+nd*j,k+nd*l) += \sum_st w D_is T_jslt D_kt
426 // or
427 // A(i+nd*j,k+nd*l) += \sum_st w D_is (X_js Y_lt + Y_js X_lt) D_kt
428 // or
429 // A(i+nd*j,k+nd*l) +=
430 // \sum_st w [ (D X^t)_ij (D Y^t)_kl + (D Y^t)_ij (D X^t)_kl ]
431 void Assemble_TProd(scalar_t w, const scalar_t *X, const scalar_t *Y,
432 scalar_t *A)
433 {
434 Eval_DZt(X, &DXt);
435 Eval_DZt(Y, &DYt);
436 const int nd = D_height;
437 const int ah = 2*nd;
438
439 for (int i = 0; i < ah; i++)
440 {
441 const scalar_t axi = w*DXt[i], ayi = w*DYt[i];
442 A[i+ah*i] += 2*axi*DYt[i];
443 for (int j = 0; j < i; j++)
444 {
445 const scalar_t A_ij = axi*DYt[j] + ayi*DXt[j];
446 A[i+ah*j] += A_ij;
447 A[j+ah*i] += A_ij;
448 }
449 }
450 }
451
452 // Assemble the contribution from the term: T_ijkl = X_ij X_kl, where X is a
453 // pointer to a 2x2 matrix stored in column-major layout.
454 //
455 // The contribution to the matrix A is given by:
456 // A(i+nd*j,k+nd*l) += \sum_st w D_is X_js X_lt D_kt
457 // or
458 // A(i+nd*j,k+nd*l) += \sum_st w [ (D X^t)_ij (D X^t)_kl ]
459 void Assemble_TProd(scalar_t w, const scalar_t *X, scalar_t *A)
460 {
461 Eval_DZt(X, &DXt);
462 const int nd = D_height;
463 const int ah = 2*nd;
464
465 for (int i = 0; i < ah; i++)
466 {
467 const scalar_t axi = w*DXt[i];
468 A[i+ah*i] += axi*DXt[i];
469 for (int j = 0; j < i; j++)
470 {
471 const scalar_t A_ij = axi*DXt[j];
472 A[i+ah*j] += A_ij;
473 A[j+ah*i] += A_ij;
474 }
475 }
476 }
477};
478
479
480/** @brief Auxiliary class for evaluating the 3x3 matrix invariants and their
481 first and second derivatives. */
482/**
483 The type `scalar_t` must support the standard operations:
484
485 =, +=, -=, +, -, *, /, unary -, int*scalar_t, int/scalar_t, scalar_t/int
486
487 The type `scalar_ops` must define the static methods:
488
489 scalar_t sign(const scalar_t &);
490 scalar_t pow(const scalar_t &x, int a, int b); // x^(a/b)
491*/
492template <typename scalar_t, typename scalar_ops = ScalarOps<scalar_t> >
493class InvariantsEvaluator3D
494{
495protected:
496 // Transformation Jacobian
497 const scalar_t *J;
498
499 // Invatiants: I1b = det(J)^{-2/3} * ||J||_F^2
500 // I2b = det(J)^{ 2/3} * ||J^{-1}||_F^2
501 // I3b = det(J)
502 // Computed as:
503 // I_1 = ||J||_F^2, \bar{I}_1 = det(J)^{-2/3}*I_1,
504 // I_2 = (1/2)*(||J||_F^4-||J J^t||_F^2) = (1/2)*(I_1^2-||J J^t||_F^2),
505 // \bar{I}_2 = det(J)^{-4/3}*I_2,
506 // I_3 = det(J)^2, \bar{I}_3 = det(J).
507 scalar_t I1, I1b, I2, I2b, I3b;
508 scalar_t I3b_p; // I3b^{-2/3}
509
510 // Derivatives of I1, I1b, I2, I2b, I3, and I3b using column-major storage.
511 scalar_t dI1[9], dI1b[9], dI2[9], dI2b[9], dI3[9], dI3b[9];
512 scalar_t B[6]; // B = J J^t (diagonal entries first, then off-diagonal)
513
514 int D_height, alloc_height;
515 const scalar_t *D; // Always points to external data or is empty
516 scalar_t *DaJ, *DJt, *DdI2t, *DXt, *DYt;
517
518 enum EvalMasks
519 {
520 HAVE_I1 = 1,
521 HAVE_I1b = 2,
522 HAVE_B_offd = 4,
523 HAVE_I2 = 8,
524 HAVE_I2b = 16,
525 HAVE_I3b = 1<<5,
526 HAVE_I3b_p = 1<<6,
527 HAVE_dI1 = 1<<7,
528 HAVE_dI1b = 1<<8,
529 HAVE_dI2 = 1<<9,
530 HAVE_dI2b = 1<<10,
531 HAVE_dI3 = 1<<11,
532 HAVE_dI3b = 1<<12,
533 HAVE_DaJ = 1<<13, // D adj(J) = D dI3b^t
534 HAVE_DJt = 1<<14, // D J^t
535 HAVE_DdI2t = 1<<15 // D dI2^t
536 };
537
538 // Bitwise OR of EvalMasks
539 int eval_state;
540
541 bool dont(int have_mask) const { return !(eval_state & have_mask); }
542
543 void Eval_I1()
544 {
545 eval_state |= HAVE_I1;
546 B[0] = J[0]*J[0] + J[3]*J[3] + J[6]*J[6];
547 B[1] = J[1]*J[1] + J[4]*J[4] + J[7]*J[7];
548 B[2] = J[2]*J[2] + J[5]*J[5] + J[8]*J[8];
549 I1 = B[0] + B[1] + B[2];
550 }
551 void Eval_I1b() // det(J)^{-2/3}*I_1 = I_1/I_3^{1/3}
552 {
553 eval_state |= HAVE_I1b;
554 I1b = Get_I1()*Get_I3b_p();
555 }
556 void Eval_B_offd()
557 {
558 eval_state |= HAVE_B_offd;
559 // B = J J^t
560 // B[3]=B(0,1), B[4]=B(0,2), B[5]=B(1,2)
561 B[3] = J[0]*J[1] + J[3]*J[4] + J[6]*J[7]; // B(0,1)
562 B[4] = J[0]*J[2] + J[3]*J[5] + J[6]*J[8]; // B(0,2)
563 B[5] = J[1]*J[2] + J[4]*J[5] + J[7]*J[8]; // B(1,2)
564 }
565 void Eval_I2()
566 {
567 eval_state |= HAVE_I2;
568 Get_I1();
569 if (dont(HAVE_B_offd)) { Eval_B_offd(); }
570 const scalar_t BF2 = B[0]*B[0] + B[1]*B[1] + B[2]*B[2] +
571 2*(B[3]*B[3] + B[4]*B[4] + B[5]*B[5]);
572 I2 = (I1*I1 - BF2)/2;
573 }
574 void Eval_I2b() // I2b = I2*I3b^{-4/3}
575 {
576 eval_state |= HAVE_I2b;
577 Get_I3b_p();
578 I2b = Get_I2()*I3b_p*I3b_p;
579 }
580 void Eval_I3b() // det(J)
581 {
582 eval_state |= HAVE_I3b;
583 I3b = J[0]*(J[4]*J[8] - J[7]*J[5]) - J[1]*(J[3]*J[8] - J[5]*J[6]) +
584 J[2]*(J[3]*J[7] - J[4]*J[6]);
585 }
586 scalar_t Get_I3b_p() // I3b^{-2/3}
587 {
588 if (dont(HAVE_I3b_p))
589 {
590 eval_state |= HAVE_I3b_p;
591 const scalar_t i3b = Get_I3b();
592 I3b_p = scalar_ops::pow(i3b, -2, 3);
593 }
594 return I3b_p;
595 }
596 void Eval_dI1()
597 {
598 eval_state |= HAVE_dI1;
599 for (int i = 0; i < 9; i++)
600 {
601 dI1[i] = 2*J[i];
602 }
603 }
604 void Eval_dI1b()
605 {
606 eval_state |= HAVE_dI1b;
607 // I1b = I3b^{-2/3}*I1
608 // dI1b = 2*I3b^{-2/3}*(J - (1/3)*I1/I3b*dI3b)
609 const scalar_t c1 = 2*Get_I3b_p();
610 const scalar_t c2 = Get_I1()/(3*I3b);
611 Get_dI3b();
612 for (int i = 0; i < 9; i++)
613 {
614 dI1b[i] = c1*(J[i] - c2*dI3b[i]);
615 }
616 }
617 void Eval_dI2()
618 {
619 eval_state |= HAVE_dI2;
620 // dI2 = 2 I_1 J - 2 J J^t J = 2 (I_1 I - B) J
621 Get_I1();
622 if (dont(HAVE_B_offd)) { Eval_B_offd(); }
623 // B[0]=B(0,0), B[1]=B(1,1), B[2]=B(2,2)
624 // B[3]=B(0,1), B[4]=B(0,2), B[5]=B(1,2)
625 const scalar_t C[6] =
626 {
627 2*(I1 - B[0]), 2*(I1 - B[1]), 2*(I1 - B[2]),
628 -2*B[3], -2*B[4], -2*B[5]
629 };
630 // | C[0] C[3] C[4] | | J[0] J[3] J[6] |
631 // dI2 = | C[3] C[1] C[5] | | J[1] J[4] J[7] |
632 // | C[4] C[5] C[2] | | J[2] J[5] J[8] |
633 dI2[0] = C[0]*J[0] + C[3]*J[1] + C[4]*J[2];
634 dI2[1] = C[3]*J[0] + C[1]*J[1] + C[5]*J[2];
635 dI2[2] = C[4]*J[0] + C[5]*J[1] + C[2]*J[2];
636
637 dI2[3] = C[0]*J[3] + C[3]*J[4] + C[4]*J[5];
638 dI2[4] = C[3]*J[3] + C[1]*J[4] + C[5]*J[5];
639 dI2[5] = C[4]*J[3] + C[5]*J[4] + C[2]*J[5];
640
641 dI2[6] = C[0]*J[6] + C[3]*J[7] + C[4]*J[8];
642 dI2[7] = C[3]*J[6] + C[1]*J[7] + C[5]*J[8];
643 dI2[8] = C[4]*J[6] + C[5]*J[7] + C[2]*J[8];
644 }
645 void Eval_dI2b()
646 {
647 eval_state |= HAVE_dI2b;
648 // I2b = det(J)^{-4/3}*I2 = I3b^{-4/3}*I2
649 // dI2b = (-4/3)*I3b^{-7/3}*I2*dI3b + I3b^{-4/3}*dI2
650 // = I3b^{-4/3} * [ dI2 - (4/3)*I2/I3b*dI3b ]
651 Get_I3b_p();
652 const scalar_t c1 = I3b_p*I3b_p;
653 const scalar_t c2 = (4*Get_I2()/I3b)/3;
654 Get_dI2();
655 Get_dI3b();
656 for (int i = 0; i < 9; i++)
657 {
658 dI2b[i] = c1*(dI2[i] - c2*dI3b[i]);
659 }
660 }
661 void Eval_dI3()
662 {
663 eval_state |= HAVE_dI3;
664 // I3 = I3b^2
665 // dI3 = 2*I3b*dI3b = 2*det(J)*adj(J)^T
666 const scalar_t c1 = 2*Get_I3b();
667 Get_dI3b();
668 for (int i = 0; i < 9; i++)
669 {
670 dI3[i] = c1*dI3b[i];
671 }
672 }
673 void Eval_dI3b()
674 {
675 eval_state |= HAVE_dI3b;
676 // I3b = det(J)
677 // dI3b = adj(J)^T
678 dI3b[0] = J[4]*J[8] - J[5]*J[7]; // 0 3 6
679 dI3b[1] = J[5]*J[6] - J[3]*J[8]; // 1 4 7
680 dI3b[2] = J[3]*J[7] - J[4]*J[6]; // 2 5 8
681 dI3b[3] = J[2]*J[7] - J[1]*J[8];
682 dI3b[4] = J[0]*J[8] - J[2]*J[6];
683 dI3b[5] = J[1]*J[6] - J[0]*J[7];
684 dI3b[6] = J[1]*J[5] - J[2]*J[4];
685 dI3b[7] = J[2]*J[3] - J[0]*J[5];
686 dI3b[8] = J[0]*J[4] - J[1]*J[3];
687 }
688 void Eval_DZt(const scalar_t *Z, scalar_t **DZt_ptr)
689 {
690 MFEM_ASSERT(D != NULL, "");
691 const int nd = D_height;
692 scalar_t *DZt = *DZt_ptr;
693 if (DZt == NULL) { *DZt_ptr = DZt = new scalar_t[3*alloc_height]; }
694 for (int i = 0; i < nd; i++)
695 {
696 const int i0 = i+nd*0, i1 = i+nd*1, i2 = i+nd*2;
697 DZt[i0] = D[i0]*Z[0] + D[i1]*Z[3] + D[i2]*Z[6];
698 DZt[i1] = D[i0]*Z[1] + D[i1]*Z[4] + D[i2]*Z[7];
699 DZt[i2] = D[i0]*Z[2] + D[i1]*Z[5] + D[i2]*Z[8];
700 }
701 }
702 void Eval_DaJ() // DaJ = D adj(J) = D dI3b^t
703 {
704 eval_state |= HAVE_DaJ;
705 Get_dI3b();
706 Eval_DZt(dI3b, &DaJ);
707 }
708 void Eval_DJt() // DJt = D J^t
709 {
710 eval_state |= HAVE_DJt;
711 Eval_DZt(J, &DJt);
712 }
713 void Eval_DdI2t() // DdI2t = D dI2^t
714 {
715 eval_state |= HAVE_DdI2t;
716 Get_dI2();
717 Eval_DZt(dI2, &DdI2t);
718 }
719
720public:
721 /// The Jacobian should use column-major storage.
722 InvariantsEvaluator3D(const scalar_t *Jac = NULL)
723 : J(Jac), D_height(), alloc_height(),
724 D(), DaJ(), DJt(), DdI2t(), DXt(), DYt(), eval_state(0) { }
725
726 ~InvariantsEvaluator3D()
727 {
728 delete [] DYt;
729 delete [] DXt;
730 delete [] DdI2t;
731 delete [] DJt;
732 delete [] DaJ;
733 }
734
735 /// The Jacobian should use column-major storage.
736 void SetJacobian(const scalar_t *Jac) { J = Jac; eval_state = 0; }
737
738 /// The @a Deriv matrix is `dof x 3`, using column-major storage.
739 void SetDerivativeMatrix(int height, const scalar_t *Deriv)
740 {
741 eval_state &= ~(HAVE_DaJ | HAVE_DJt | HAVE_DdI2t);
742 if (alloc_height < height)
743 {
744 delete [] DYt; DYt = NULL;
745 delete [] DXt; DXt = NULL;
746 delete [] DdI2t; DdI2t = NULL;
747 delete [] DJt; DJt = NULL;
748 delete [] DaJ; DaJ = NULL;
749 alloc_height = height;
750 }
751 D_height = height;
752 D = Deriv;
753 }
754
755 scalar_t Get_I1() { if (dont(HAVE_I1 )) { Eval_I1(); } return I1; }
756 scalar_t Get_I1b() { if (dont(HAVE_I1b)) { Eval_I1b(); } return I1b; }
757 scalar_t Get_I2() { if (dont(HAVE_I2 )) { Eval_I2(); } return I2; }
758 scalar_t Get_I2b() { if (dont(HAVE_I2b)) { Eval_I2b(); } return I2b; }
759 scalar_t Get_I3() { if (dont(HAVE_I3b)) { Eval_I3b(); } return I3b*I3b; }
760 scalar_t Get_I3b() { if (dont(HAVE_I3b)) { Eval_I3b(); } return I3b; }
761
762 const scalar_t *Get_dI1()
763 {
764 if (dont(HAVE_dI1 )) { Eval_dI1(); } return dI1;
765 }
766 const scalar_t *Get_dI1b()
767 {
768 if (dont(HAVE_dI1b)) { Eval_dI1b(); } return dI1b;
769 }
770 const scalar_t *Get_dI2()
771 {
772 if (dont(HAVE_dI2)) { Eval_dI2(); } return dI2;
773 }
774 const scalar_t *Get_dI2b()
775 {
776 if (dont(HAVE_dI2b)) { Eval_dI2b(); } return dI2b;
777 }
778 const scalar_t *Get_dI3()
779 {
780 if (dont(HAVE_dI3)) { Eval_dI3(); } return dI3;
781 }
782 const scalar_t *Get_dI3b()
783 {
784 if (dont(HAVE_dI3b)) { Eval_dI3b(); } return dI3b;
785 }
786
787 // Assemble operation for tensor X with components X_jslt:
788 // A(i+nd*j,k+nd*l) += (\sum_st w D_is X_jslt D_kt)
789 // 0 <= i,k < nd, 0 <= j,l,s,t < 3
790 // where nd is the height of D, i.e. the number of DOFs in one component.
791
792 void Assemble_ddI1(scalar_t w, scalar_t *A)
793 {
794 // ddI1_jslt = 2 I_jslt = 2 δ_jl δ_st
795 // A(i+nd*j,k+nd*l) += (\sum_st 2 w D_is δ_jl δ_st D_kt)
796 // or
797 // A(i+nd*j,k+nd*l) += (2 w) (\sum_s D_is D_ks) δ_jl
798 // A(i+nd*j,k+nd*l) += (2 w) (D D^t)_ik δ_jl
799
800 const int nd = D_height;
801 const int ah = 3*nd;
802 const scalar_t a = 2*w;
803 for (int i = 0; i < nd; i++)
804 {
805 const int i0 = i+nd*0, i1 = i+nd*1, i2 = i+nd*2;
806 const scalar_t aDi[3] = { a*D[i0], a*D[i1], a*D[i2] };
807 // k == i
808 const scalar_t aDDt_ii = aDi[0]*D[i0] + aDi[1]*D[i1] + aDi[2]*D[i2];
809 A[i0+ah*i0] += aDDt_ii;
810 A[i1+ah*i1] += aDDt_ii;
811 A[i2+ah*i2] += aDDt_ii;
812 // 0 <= k < i
813 for (int k = 0; k < i; k++)
814 {
815 const int k0 = k+nd*0, k1 = k+nd*1, k2 = k+nd*2;
816 const scalar_t aDDt_ik = aDi[0]*D[k0] + aDi[1]*D[k1] + aDi[2]*D[k2];
817 A[i0+ah*k0] += aDDt_ik;
818 A[k0+ah*i0] += aDDt_ik;
819 A[i1+ah*k1] += aDDt_ik;
820 A[k1+ah*i1] += aDDt_ik;
821 A[i2+ah*k2] += aDDt_ik;
822 A[k2+ah*i2] += aDDt_ik;
823 }
824 }
825 }
826 void Assemble_ddI1b(scalar_t w, scalar_t *A)
827 {
828 // Similar to InvariantsEvaluator2D::Assemble_ddI1b():
829 //
830 // ddI1b = X1 + X2 + X3, where
831 // X1_ijkl = (2/3*I1b/I3) [ (2/3 δ_ks δ_it + δ_kt δ_si) dI3b_tj dI3b_sl ]
832 // = (2/3*I1b/I3) [ 2/3 dI3b_ij dI3b_kl + dI3b_kj dI3b_il ]
833 // X2_ijkl = (2*I3b^{-2/3}) δ_ik δ_jl = (I3b^{-2/3}) ddI1_ijkl
834 // X3_ijkl = -(4/3*I3b^{-5/3}) (δ_ks δ_it) (J_tj dI3b_sl + dI3b_tj J_sl)
835 // = -(4/3*I3b^{-5/3}) (J_ij dI3b_kl + dI3b_ij J_kl)
836 //
837 // A(i+nd*j,k+nd*l) += (\sum_st w D_is ddI1b_jslt D_kt)
838 // or
839 // A(i+nd*j,k+nd*l) +=
840 // w (2/3*I1b/I3) [ 2/3 DaJ_ij DaJ_kl + DaJ_il DaJ_kj ]
841 // + w (2*I3b^{-2/3}) (D D^t)_ik δ_jl
842 // - w (4/3*I3b^{-5/3}) [ DJt_ij DaJ_kl + DaJ_ij DJt_kl ]
843
844 if (dont(HAVE_DaJ)) { Eval_DaJ(); }
845 if (dont(HAVE_DJt)) { Eval_DJt(); }
846 const int nd = D_height;
847 const int ah = 3*nd;
848 const scalar_t r23 = scalar_t(2)/3;
849 const scalar_t r53 = scalar_t(5)/3;
850 const scalar_t a = r23*w*Get_I1b()/Get_I3();
851 const scalar_t b = 2*w*Get_I3b_p();
852 const scalar_t c = -r23*b/I3b;
853 for (int i = 0; i < nd; i++)
854 {
855 // A1a_ik_jl = 2/3 a DaJ_ij DaJ_kl, A1b_ik_jl = a DaJ_il DaJ_kj
856 // Symmetries: A1a_ik_jl = A1a_ki_lj = 2/3 A1b_ik_lj = 2/3 A1b_ki_jl
857 // A1_ik_jl = A1_ki_lj = A1b_ik_jl + 2/3 A1b_ik_lj
858 // A1_ik_lj = A1_ki_jl = 2/3 A1b_ik_jl + A1b_ik_lj
859 // k == i:
860 // A1_ii_jl = A1_ii_lj = (5/3) a DaJ_ij DaJ_il
861 // l == j:
862 // A1_ik_jj = A1_ki_jj = (5/3) a DaJ_ij DaJ_kj
863 // k == i && l == j:
864 // A1_ii_jj = (5/3) a DaJ_ij^2
865
866 // A2_ik_jl = b (D D^t)_ik δ_jl
867 // Symmetries:
868
869 // A3_ik_jl = c [ DJt_ij DaJ_kl + DaJ_ij DJt_kl ]
870 // Symmetries:
871 // A3_ik_jl = A3_ki_lj = c [ DJt_ij DaJ_kl + DaJ_ij DJt_kl ]
872 // A3_ik_lj = A3_ki_jl = c [ DJt_il DaJ_kj + DaJ_il DJt_kj ]
873 // k == i:
874 // A3_ii_jl = A3_ii_lj = c [ DJt_ij DaJ_il + DaJ_ij DJt_il ]
875 // l == j:
876 // A3_ik_jj = A3_ki_jj = c [ DJt_ij DaJ_kj + DaJ_ij DJt_kj ]
877 // k == i && l == j:
878 // A3_ii_jj = 2 c DJt_ij DaJ_ij
879
880 const int i0 = i+nd*0, i1 = i+nd*1, i2 = i+nd*2;
881 const scalar_t aDaJ_i[3] = { a*DaJ[i0], a*DaJ[i1], a*DaJ[i2] };
882 const scalar_t bD_i[3] = { b*D[i0], b*D[i1], b*D[i2] };
883 const scalar_t cDJt_i[3] = { c*DJt[i0], c*DJt[i1], c*DJt[i2] };
884 const scalar_t cDaJ_i[3] = { c*DaJ[i0], c*DaJ[i1], c*DaJ[i2] };
885 // k == i
886 {
887 // Symmetries: A2_ii_00 = A2_ii_11 = A2_ii_22
888 const scalar_t A2_ii = bD_i[0]*D[i0]+bD_i[1]*D[i1]+bD_i[2]*D[i2];
889 A[i0+ah*i0] += (r53*aDaJ_i[0] + 2*cDJt_i[0])*DaJ[i0] + A2_ii;
890 A[i1+ah*i1] += (r53*aDaJ_i[1] + 2*cDJt_i[1])*DaJ[i1] + A2_ii;
891 A[i2+ah*i2] += (r53*aDaJ_i[2] + 2*cDJt_i[2])*DaJ[i2] + A2_ii;
892
893 // Symmetries: A_ii_jl = A_ii_lj
894 for (int j = 1; j < 3; j++)
895 {
896 const int ij = i+nd*j;
897 for (int l = 0; l < j; l++)
898 {
899 const int il = i+nd*l;
900 const scalar_t A_ii_jl =
901 (r53*aDaJ_i[j] + cDJt_i[j])*DaJ[il] + cDaJ_i[j]*DJt[il];
902 A[ij+ah*il] += A_ii_jl;
903 A[il+ah*ij] += A_ii_jl;
904 }
905 }
906 }
907 // 0 <= k < i
908 for (int k = 0; k < i; k++)
909 {
910 const int k0 = k+nd*0, k1 = k+nd*1, k2 = k+nd*2;
911 // Symmetries: A2_ik_jj = A2_ki_ll
912 const scalar_t A2_ik = bD_i[0]*D[k0]+bD_i[1]*D[k1]+bD_i[2]*D[k2];
913
914 // l == j
915 for (int j = 0; j < 3; j++)
916 {
917 const int ij = i+nd*j, kj = k+nd*j;
918 const scalar_t A_ik_jj = (r53*aDaJ_i[j] + cDJt_i[j])*DaJ[kj] +
919 cDaJ_i[j]*DJt[kj] + A2_ik;
920 A[ij+ah*kj] += A_ik_jj;
921 A[kj+ah*ij] += A_ik_jj;
922 }
923
924 // 0 <= l < j
925 for (int j = 1; j < 3; j++)
926 {
927 const int ij = i+nd*j, kj = k+nd*j;
928 for (int l = 0; l < j; l++)
929 {
930 const int il = i+nd*l, kl = k+nd*l;
931 // A1b_ik_jl = a DaJ_il DaJ_kj
932 const scalar_t A1b_ik_jl = aDaJ_i[l]*DaJ[kj];
933 const scalar_t A1b_ik_lj = aDaJ_i[j]*DaJ[kl];
934 // A1_ik_jl = A1_ki_lj = A1b_ik_jl + 2/3 A1b_ik_lj
935 // A1_ik_lj = A1_ki_jl = 2/3 A1b_ik_jl + A1b_ik_lj
936 // A3_ik_jl = c [ DJt_ij DaJ_kl + DaJ_ij DJt_kl ]
937 const scalar_t A_ik_jl = A1b_ik_jl + r23*A1b_ik_lj +
938 cDJt_i[j]*DaJ[kl]+cDaJ_i[j]*DJt[kl];
939 A[ij+ah*kl] += A_ik_jl;
940 A[kl+ah*ij] += A_ik_jl;
941 const scalar_t A_ik_lj = r23*A1b_ik_jl + A1b_ik_lj +
942 cDJt_i[l]*DaJ[kj]+cDaJ_i[l]*DJt[kj];
943 A[il+ah*kj] += A_ik_lj;
944 A[kj+ah*il] += A_ik_lj;
945 }
946 }
947 }
948 }
949 }
950 void Assemble_ddI2(scalar_t w, scalar_t *A)
951 {
952 // dI2 = 2 I_1 J - 2 J J^t J = 2 (I_1 I - B) J
953 //
954 // ddI2 = X1 + X2 + X3
955 // X1_ijkl = (2 I_1) δ_ik δ_jl
956 // X2_ijkl = 2 ( 2 δ_ku δ_iv - δ_ik δ_uv - δ_kv δ_iu ) J_vj J_ul
957 // X3_ijkl = -2 (J J^t)_ik δ_jl = -2 B_ik δ_jl
958 //
959 // Apply: j->s, i->j, l->t, k->l
960 // X2_jslt = 2 ( δ_lu δ_jv - δ_jl δ_uv +
961 // δ_lu δ_jv - δ_lv δ_ju ) J_vs J_ut
962 //
963 // A(i+nd*j,k+nd*l) += (\sum_st w D_is ddI2_jslt D_kt)
964 //
965 // \sum_st w D_is X1_jslt D_kt =
966 // \sum_st w D_is [ (2 I_1) δ_jl δ_st ] D_kt =
967 // (2 w I_1) D_is δ_jl D_ks = (2 w I_1) (D D^t)_ik δ_jl
968 //
969 // \sum_st w D_is X2_jslt D_kt =
970 // \sum_stuv w D_is [ 2 ( δ_lu δ_jv - δ_jl δ_uv +
971 // δ_lu δ_jv - δ_lv δ_ju ) J_vs J_ut ] D_kt =
972 // \sum_uv 2 w [ δ_lu δ_jv - δ_jl δ_uv +
973 // δ_lu δ_jv - δ_lv δ_ju ] (D J^t)_iv (D J^t)_ku =
974 // 2 w ( DJt_ij DJt_kl - δ_jl (DJt DJt^t)_ik ) +
975 // 2 w ( DJt_ij DJt_kl - DJt_il DJt_kj )
976 //
977 // \sum_st w D_is X3_jslt D_kt = \sum_st w D_is [ -2 B_jl δ_st ] D_kt =
978 // -2 w (D D^t)_ik B_jl
979 //
980 // A(i+nd*j,k+nd*l) +=
981 // (2 w I_1) (D D^t)_ik δ_jl - 2 w (D D^t)_ik B_jl +
982 // 2 w DJt_ij DJt_kl - 2 w (DJt DJt^t)_ik δ_jl +
983 // 2 w ( DJt_ij DJt_kl - DJt_il DJt_kj )
984 //
985 // The last term is a determinant: rows {i,k} and columns {j,l} of DJt:
986 // | DJt_ij DJt_il |
987 // | DJt_kj DJt_kl | = DJt_ij DJt_kl - DJt_il DJt_kj
988
989 if (dont(HAVE_DJt)) { Eval_DJt(); }
990 Get_I1(); // evaluates I1 and the diagonal of B
991 if (dont(HAVE_B_offd)) { Eval_B_offd(); }
992 const int nd = D_height;
993 const int ah = 3*nd;
994 const scalar_t a = 2*w;
995 for (int i = 0; i < ah; i++)
996 {
997 const scalar_t avi = a*DJt[i];
998 A[i+ah*i] += avi*DJt[i];
999 for (int j = 0; j < i; j++)
1000 {
1001 const scalar_t aVVt_ij = avi*DJt[j];
1002 A[i+ah*j] += aVVt_ij;
1003 A[j+ah*i] += aVVt_ij;
1004 }
1005 }
1006
1007 for (int i = 0; i < nd; i++)
1008 {
1009 const int i0 = i+nd*0, i1 = i+nd*1, i2 = i+nd*2;
1010 const scalar_t aD_i[3] = { a*D[i0], a*D[i1], a*D[i2] };
1011 const scalar_t aDJt_i[3] = { a*DJt[i0], a*DJt[i1], a*DJt[i2] };
1012 // k == i
1013 {
1014 const scalar_t aDDt_ii =
1015 aD_i[0]*D[i0] + aD_i[1]*D[i1] + aD_i[2]*D[i2];
1016 const scalar_t Z1_ii =
1017 I1*aDDt_ii - (aDJt_i[0]*DJt[i0] + aDJt_i[1]*DJt[i1] +
1018 aDJt_i[2]*DJt[i2]);
1019 // l == j
1020 for (int j = 0; j < 3; j++)
1021 {
1022 const int ij = i+nd*j;
1023 A[ij+ah*ij] += Z1_ii - aDDt_ii*B[j];
1024 }
1025 // l != j
1026 const scalar_t Z2_ii_01 = aDDt_ii*B[3];
1027 const scalar_t Z2_ii_02 = aDDt_ii*B[4];
1028 const scalar_t Z2_ii_12 = aDDt_ii*B[5];
1029 A[i0+ah*i1] -= Z2_ii_01;
1030 A[i1+ah*i0] -= Z2_ii_01;
1031 A[i0+ah*i2] -= Z2_ii_02;
1032 A[i2+ah*i0] -= Z2_ii_02;
1033 A[i1+ah*i2] -= Z2_ii_12;
1034 A[i2+ah*i1] -= Z2_ii_12;
1035 }
1036 // 0 <= k < i
1037 for (int k = 0; k < i; k++)
1038 {
1039 const int k0 = k+nd*0, k1 = k+nd*1, k2 = k+nd*2;
1040 const scalar_t aDDt_ik =
1041 aD_i[0]*D[k0] + aD_i[1]*D[k1] + aD_i[2]*D[k2];
1042 const scalar_t Z1_ik =
1043 I1*aDDt_ik - (aDJt_i[0]*DJt[k0] + aDJt_i[1]*DJt[k1] +
1044 aDJt_i[2]*DJt[k2]);
1045 // l == j
1046 for (int j = 0; j < 3; j++)
1047 {
1048 const int ij = i+nd*j, kj = k+nd*j;
1049 const scalar_t Z2_ik_jj = Z1_ik - aDDt_ik*B[j];
1050 A[ij+ah*kj] += Z2_ik_jj;
1051 A[kj+ah*ij] += Z2_ik_jj;
1052 }
1053 // l != j
1054 {
1055 const scalar_t Z2_ik_01 = aDDt_ik*B[3];
1056 A[i0+ah*k1] -= Z2_ik_01;
1057 A[i1+ah*k0] -= Z2_ik_01;
1058 A[k0+ah*i1] -= Z2_ik_01;
1059 A[k1+ah*i0] -= Z2_ik_01;
1060 const scalar_t Z2_ik_02 = aDDt_ik*B[4];
1061 A[i0+ah*k2] -= Z2_ik_02;
1062 A[i2+ah*k0] -= Z2_ik_02;
1063 A[k0+ah*i2] -= Z2_ik_02;
1064 A[k2+ah*i0] -= Z2_ik_02;
1065 const scalar_t Z2_ik_12 = aDDt_ik*B[5];
1066 A[i1+ah*k2] -= Z2_ik_12;
1067 A[i2+ah*k1] -= Z2_ik_12;
1068 A[k1+ah*i2] -= Z2_ik_12;
1069 A[k2+ah*i1] -= Z2_ik_12;
1070 }
1071 // 0 <= l < j
1072 for (int j = 1; j < 3; j++)
1073 {
1074 const int ij = i+nd*j, kj = k+nd*j;
1075 for (int l = 0; l < j; l++)
1076 {
1077 const int il = i+nd*l, kl = k+nd*l;
1078 const scalar_t Z3_ik_jl =
1079 aDJt_i[j]*DJt[kl] - aDJt_i[l]*DJt[kj];
1080 A[ij+ah*kl] += Z3_ik_jl;
1081 A[kl+ah*ij] += Z3_ik_jl;
1082 A[il+ah*kj] -= Z3_ik_jl;
1083 A[kj+ah*il] -= Z3_ik_jl;
1084 }
1085 }
1086 }
1087 }
1088 }
1089 void Assemble_ddI2b(scalar_t w, scalar_t *A)
1090 {
1091 // dI2b = (-4/3)*I3b^{-7/3}*I2*dI3b + I3b^{-4/3}*dI2
1092 // = I3b^{-4/3} * [ dI2 - (4/3)*I2/I3b*dI3b ]
1093 //
1094 // ddI2b = X1 + X2 + X3
1095 // X1_ijkl = 16/9 det(J)^{-10/3} I2 dI3b_ij dI3b_kl +
1096 // 4/3 det(J)^{-10/3} I2 dI3b_il dI3b_kj
1097 // X2_ijkl = -4/3 det(J)^{-7/3} (dI2_ij dI3b_kl + dI2_kl dI3b_ij)
1098 // X3_ijkl = det(J)^{-4/3} ddI2_ijkl
1099 //
1100 // Apply: j->s, i->j, l->t, k->l
1101 // X1_jslt = 16/9 det(J)^{-10/3} I2 dI3b_js dI3b_lt +
1102 // 4/3 det(J)^{-10/3} I2 dI3b_jt dI3b_ls
1103 // X2_jslt = -4/3 det(J)^{-7/3} (dI2_js dI3b_lt + dI2_lt dI3b_js)
1104 //
1105 // A(i+nd*j,k+nd*l) += (\sum_st w D_is ddI2b_jslt D_kt)
1106 //
1107 // (\sum_st w D_is X1_jslt D_kt) =
1108 // 16/9 w det(J)^{-10/3} I2 DaJ_ij DaJ_kl +
1109 // 4/3 w det(J)^{-10/3} I2 DaJ_il DaJ_kj
1110 //
1111 // (\sum_st w D_is X1_jslt D_kt) =
1112 // -4/3 w det(J)^{-7/3} D_is (dI2_js dI3b_lt + dI2_lt dI3b_js) D_kt =
1113 // -4/3 w det(J)^{-7/3} [ (D dI2^t)_ij DaJ_kl + DaJ_ij (D dI2^t)_kl ]
1114 //
1115 // A(i+nd*j,k+nd*l) +=
1116 // 16/9 w det(J)^{-10/3} I2 DaJ_ij DaJ_kl +
1117 // 4/3 w det(J)^{-10/3} I2 DaJ_il DaJ_kj -
1118 // 4/3 w det(J)^{-7/3} [ DdI2t_ij DaJ_kl + DaJ_ij DdI2t_kl ] +
1119 // w det(J)^{-4/3} D_is D_kt ddI2_jslt
1120
1121 Get_I3b_p(); // = det(J)^{-2/3}, evaluates I3b
1122 if (dont(HAVE_DaJ)) { Eval_DaJ(); }
1123 if (dont(HAVE_DdI2t)) { Eval_DdI2t(); }
1124 const int nd = D_height;
1125 const int ah = 3*nd;
1126 const scalar_t a = w*I3b_p*I3b_p;
1127 const scalar_t b = (-4*a)/(3*I3b);
1128 const scalar_t c = -b*Get_I2()/I3b;
1129 const scalar_t d = (4*c)/3;
1130
1131 for (int i = 0; i < ah; i++)
1132 {
1133 const scalar_t dvi = d*DaJ[i];
1134 A[i+ah*i] += dvi*DaJ[i];
1135 for (int j = 0; j < i; j++)
1136 {
1137 const scalar_t dVVt_ij = dvi*DaJ[j];
1138 A[i+ah*j] += dVVt_ij;
1139 A[j+ah*i] += dVVt_ij;
1140 }
1141 }
1142 Assemble_ddI2(a, A);
1143 for (int i = 0; i < nd; i++)
1144 {
1145 const int i0 = i+nd*0, i1 = i+nd*1, i2 = i+nd*2;
1146 const scalar_t cDaJ_i[3] = { c*DaJ[i0], c*DaJ[i1], c*DaJ[i2] };
1147 const scalar_t bDaJ_i[3] = { b*DaJ[i0], b*DaJ[i1], b*DaJ[i2] };
1148 const scalar_t bDdI2t_i[3] = { b*DdI2t[i0], b*DdI2t[i1], b*DdI2t[i2] };
1149 // k == i
1150 {
1151 // l == j
1152 for (int j = 0; j < 3; j++)
1153 {
1154 const int ij = i+nd*j;
1155 A[ij+ah*ij] += (cDaJ_i[j] + 2*bDdI2t_i[j])*DaJ[ij];
1156 }
1157 // 0 <= l < j
1158 for (int j = 1; j < 3; j++)
1159 {
1160 const int ij = i+nd*j;
1161 for (int l = 0; l < j; l++)
1162 {
1163 const int il = i+nd*l;
1164 const scalar_t Z_ii_jl =
1165 (cDaJ_i[l] + bDdI2t_i[l])*DaJ[ij] + bDdI2t_i[j]*DaJ[il];
1166 A[ij+ah*il] += Z_ii_jl;
1167 A[il+ah*ij] += Z_ii_jl;
1168 }
1169 }
1170 }
1171 // 0 <= k < i
1172 for (int k = 0; k < i; k++)
1173 {
1174 // l == j
1175 for (int j = 0; j < 3; j++)
1176 {
1177 const int ij = i+nd*j, kj = k+nd*j;
1178 const scalar_t Z_ik_jj =
1179 (cDaJ_i[j] + bDdI2t_i[j])*DaJ[kj] + bDaJ_i[j]*DdI2t[kj];
1180 A[ij+ah*kj] += Z_ik_jj;
1181 A[kj+ah*ij] += Z_ik_jj;
1182 }
1183 // 0 <= l < j
1184 for (int j = 1; j < 3; j++)
1185 {
1186 const int ij = i+nd*j, kj = k+nd*j;
1187 for (int l = 0; l < j; l++)
1188 {
1189 const int il = i+nd*l, kl = k+nd*l;
1190 const scalar_t Z_ik_jl = cDaJ_i[l]*DaJ[kj] +
1191 bDdI2t_i[j]*DaJ[kl] +
1192 bDaJ_i[j]*DdI2t[kl];
1193 A[ij+ah*kl] += Z_ik_jl;
1194 A[kl+ah*ij] += Z_ik_jl;
1195 const scalar_t Z_ik_lj = cDaJ_i[j]*DaJ[kl] +
1196 bDdI2t_i[l]*DaJ[kj] +
1197 bDaJ_i[l]*DdI2t[kj];
1198 A[il+ah*kj] += Z_ik_lj;
1199 A[kj+ah*il] += Z_ik_lj;
1200 }
1201 }
1202 }
1203 }
1204 }
1205 void Assemble_ddI3(scalar_t w, scalar_t *A)
1206 {
1207 // Similar to InvariantsEvaluator2D::Assemble_ddI2():
1208 //
1209 // A(i+nd*j,k+nd*l) += 2 w [ 2 DaJ_ij DaJ_kl - DaJ_il DaJ_kj ]
1210 //
1211 // Note: the expression ( DaJ_ij DaJ_kl - DaJ_il DaJ_kj ) is the
1212 // determinant of the 2x2 matrix formed by rows {i,k} and columns {j,l}
1213 // from the matrix DaJ = D dI3b^t.
1214
1215 if (dont(HAVE_DaJ)) { Eval_DaJ(); }
1216 const int nd = D_height;
1217 const int ah = 3*nd;
1218 const scalar_t a = 2*w;
1219
1220 for (int i = 0; i < ah; i++)
1221 {
1222 const scalar_t avi = a*DaJ[i];
1223 A[i+ah*i] += avi*DaJ[i];
1224 for (int j = 0; j < i; j++)
1225 {
1226 const scalar_t aVVt_ij = avi*DaJ[j];
1227 A[i+ah*j] += aVVt_ij;
1228 A[j+ah*i] += aVVt_ij;
1229 }
1230 }
1231 for (int j = 1; j < 3; j++)
1232 {
1233 for (int l = 0; l < j; l++)
1234 {
1235 for (int i = 0; i < nd; i++)
1236 {
1237 const int ij = i+nd*j, il = i+nd*l;
1238 const scalar_t aDaJ_ij = a*DaJ[ij], aDaJ_il = a*DaJ[il];
1239 for (int k = 0; k < i; k++)
1240 {
1241 const int kj = k+nd*j, kl = k+nd*l;
1242 const scalar_t A_ijkl = aDaJ_ij*DaJ[kl] - aDaJ_il*DaJ[kj];
1243 A[ij+ah*kl] += A_ijkl;
1244 A[kl+ah*ij] += A_ijkl;
1245 A[kj+ah*il] -= A_ijkl;
1246 A[il+ah*kj] -= A_ijkl;
1247 }
1248 }
1249 }
1250 }
1251 }
1252 void Assemble_ddI3b(scalar_t w, scalar_t *A)
1253 {
1254 // Similar to InvariantsEvaluator2D::Assemble_ddI2b():
1255 //
1256 // A(i+nd*j,k+nd*l) += (w/I3b) [ DaJ_ij DaJ_kl - DaJ_il DaJ_kj ]
1257 //
1258 // | DaJ_ij DaJ_il | = determinant of rows {i,k}, columns {j,l} from DaJ
1259 // | DaJ_kj DaJ_kl |
1260
1261 if (dont(HAVE_DaJ)) { Eval_DaJ(); }
1262 const int nd = D_height;
1263 const int ah = 3*nd;
1264 const scalar_t a = w/Get_I3b();
1265 for (int j = 1; j < 3; j++)
1266 {
1267 for (int l = 0; l < j; l++)
1268 {
1269 for (int i = 0; i < nd; i++)
1270 {
1271 const int ij = i+nd*j, il = i+nd*l;
1272 const scalar_t aDaJ_ij = a*DaJ[ij], aDaJ_il = a*DaJ[il];
1273 for (int k = 0; k < i; k++)
1274 {
1275 const int kj = k+nd*j, kl = k+nd*l;
1276 const scalar_t A_ijkl = aDaJ_ij*DaJ[kl] - aDaJ_il*DaJ[kj];
1277 A[ij+ah*kl] += A_ijkl;
1278 A[kl+ah*ij] += A_ijkl;
1279 A[kj+ah*il] -= A_ijkl;
1280 A[il+ah*kj] -= A_ijkl;
1281 }
1282 }
1283 }
1284 }
1285 }
1286 // Assemble the contribution from the term: T_ijkl = X_ij Y_kl + Y_ij X_kl,
1287 // where X and Y are pointers to 3x3 matrices stored in column-major layout.
1288 //
1289 // The contribution to the matrix A is given by:
1290 // A(i+nd*j,k+nd*l) += \sum_st w D_is T_jslt D_kt
1291 // or
1292 // A(i+nd*j,k+nd*l) += \sum_st w D_is (X_js Y_lt + Y_js X_lt) D_kt
1293 // or
1294 // A(i+nd*j,k+nd*l) +=
1295 // \sum_st w [ (D X^t)_ij (D Y^t)_kl + (D Y^t)_ij (D X^t)_kl ]
1296 void Assemble_TProd(scalar_t w, const scalar_t *X, const scalar_t *Y,
1297 scalar_t *A)
1298 {
1299 Eval_DZt(X, &DXt);
1300 Eval_DZt(Y, &DYt);
1301 const int nd = D_height;
1302 const int ah = 3*nd;
1303
1304 for (int i = 0; i < ah; i++)
1305 {
1306 const scalar_t axi = w*DXt[i], ayi = w*DYt[i];
1307 A[i+ah*i] += 2*axi*DYt[i];
1308 for (int j = 0; j < i; j++)
1309 {
1310 const scalar_t A_ij = axi*DYt[j] + ayi*DXt[j];
1311 A[i+ah*j] += A_ij;
1312 A[j+ah*i] += A_ij;
1313 }
1314 }
1315 }
1316 // Assemble the contribution from the term: T_ijkl = X_ij X_kl, where X is a
1317 // pointer to a 3x3 matrix stored in column-major layout.
1318 //
1319 // The contribution to the matrix A is given by:
1320 // A(i+nd*j,k+nd*l) += \sum_st w D_is X_js X_lt D_kt
1321 // or
1322 // A(i+nd*j,k+nd*l) += \sum_st w [ (D X^t)_ij (D X^t)_kl ]
1323 void Assemble_TProd(scalar_t w, const scalar_t *X, scalar_t *A)
1324 {
1325 Eval_DZt(X, &DXt);
1326 const int nd = D_height;
1327 const int ah = 3*nd;
1328
1329 for (int i = 0; i < ah; i++)
1330 {
1331 const scalar_t axi = w*DXt[i];
1332 A[i+ah*i] += axi*DXt[i];
1333 for (int j = 0; j < i; j++)
1334 {
1335 const scalar_t A_ij = axi*DXt[j];
1336 A[i+ah*j] += A_ij;
1337 A[j+ah*i] += A_ij;
1338 }
1339 }
1340 }
1341};
1342
1343}
1344
1345#endif
mfem::InvariantsEvaluator2D
Auxiliary class for evaluating the 2x2 matrix invariants and their first and second derivatives.
Definition invariants.hpp:49
mfem::InvariantsEvaluator2D::Get_dI1b
const scalar_t * Get_dI1b()
Definition invariants.hpp:207
mfem::InvariantsEvaluator2D::Assemble_ddI1
void Assemble_ddI1(scalar_t w, scalar_t *A)
Definition invariants.hpp:225
mfem::InvariantsEvaluator2D::dont
bool dont(int have_mask) const
Definition invariants.hpp:80
mfem::InvariantsEvaluator2D::SetJacobian
void SetJacobian(const scalar_t *Jac)
The Jacobian should use column-major storage.
Definition invariants.hpp:180
mfem::InvariantsEvaluator2D::Assemble_ddI1b
void Assemble_ddI1b(scalar_t w, scalar_t *A)
Definition invariants.hpp:256
mfem::InvariantsEvaluator2D::Get_dI2
const scalar_t * Get_dI2()
Definition invariants.hpp:211
mfem::InvariantsEvaluator2D::SetDerivativeMatrix
void SetDerivativeMatrix(int height, const scalar_t *Deriv)
The Deriv matrix is dof x 2, using column-major storage.
Definition invariants.hpp:183
mfem::InvariantsEvaluator2D::Assemble_ddI2b
void Assemble_ddI2b(scalar_t w, scalar_t *A)
Definition invariants.hpp:383
mfem::InvariantsEvaluator2D::Eval_DZt
void Eval_DZt(const scalar_t *Z, scalar_t **DZt_ptr)
Definition invariants.hpp:151
mfem::InvariantsEvaluator2D::Get_dI2b
const scalar_t * Get_dI2b()
Definition invariants.hpp:215
mfem::InvariantsEvaluator2D::Assemble_TProd
void Assemble_TProd(scalar_t w, const scalar_t *X, const scalar_t *Y, scalar_t *A)
Definition invariants.hpp:431
mfem::InvariantsEvaluator2D::Assemble_ddI2
void Assemble_ddI2(scalar_t w, scalar_t *A)
Definition invariants.hpp:334
mfem::InvariantsEvaluator2D::InvariantsEvaluator2D
InvariantsEvaluator2D(const scalar_t *Jac=NULL)
The Jacobian should use column-major storage.
Definition invariants.hpp:167
mfem::InvariantsEvaluator2D::Assemble_TProd
void Assemble_TProd(scalar_t w, const scalar_t *X, scalar_t *A)
Definition invariants.hpp:459
mfem::InvariantsEvaluator2D::Get_dI1
const scalar_t * Get_dI1()
Definition invariants.hpp:203
mfem::InvariantsEvaluator3D
Auxiliary class for evaluating the 3x3 matrix invariants and their first and second derivatives.
Definition invariants.hpp:494
mfem::InvariantsEvaluator3D::SetJacobian
void SetJacobian(const scalar_t *Jac)
The Jacobian should use column-major storage.
Definition invariants.hpp:736
mfem::InvariantsEvaluator3D::Get_dI2
const scalar_t * Get_dI2()
Definition invariants.hpp:770
mfem::InvariantsEvaluator3D::Get_dI2b
const scalar_t * Get_dI2b()
Definition invariants.hpp:774
mfem::InvariantsEvaluator3D::Get_dI3
const scalar_t * Get_dI3()
Definition invariants.hpp:778
mfem::InvariantsEvaluator3D::dont
bool dont(int have_mask) const
Definition invariants.hpp:541
mfem::InvariantsEvaluator3D::Assemble_TProd
void Assemble_TProd(scalar_t w, const scalar_t *X, const scalar_t *Y, scalar_t *A)
Definition invariants.hpp:1296
mfem::InvariantsEvaluator3D::Get_dI1
const scalar_t * Get_dI1()
Definition invariants.hpp:762
mfem::InvariantsEvaluator3D::SetDerivativeMatrix
void SetDerivativeMatrix(int height, const scalar_t *Deriv)
The Deriv matrix is dof x 3, using column-major storage.
Definition invariants.hpp:739
mfem::InvariantsEvaluator3D::Get_dI3b
const scalar_t * Get_dI3b()
Definition invariants.hpp:782
mfem::InvariantsEvaluator3D::Eval_DZt
void Eval_DZt(const scalar_t *Z, scalar_t **DZt_ptr)
Definition invariants.hpp:688
mfem::InvariantsEvaluator3D::InvariantsEvaluator3D
InvariantsEvaluator3D(const scalar_t *Jac=NULL)
The Jacobian should use column-major storage.
Definition invariants.hpp:722
mfem::InvariantsEvaluator3D::Assemble_ddI1
void Assemble_ddI1(scalar_t w, scalar_t *A)
Definition invariants.hpp:792
mfem::InvariantsEvaluator3D::Assemble_TProd
void Assemble_TProd(scalar_t w, const scalar_t *X, scalar_t *A)
Definition invariants.hpp:1323
mfem::InvariantsEvaluator3D::Assemble_ddI2
void Assemble_ddI2(scalar_t w, scalar_t *A)
Definition invariants.hpp:950
mfem::InvariantsEvaluator3D::Get_dI1b
const scalar_t * Get_dI1b()
Definition invariants.hpp:766
mfem::InvariantsEvaluator3D::Assemble_ddI3
void Assemble_ddI3(scalar_t w, scalar_t *A)
Definition invariants.hpp:1205
mfem::InvariantsEvaluator3D::Assemble_ddI2b
void Assemble_ddI2b(scalar_t w, scalar_t *A)
Definition invariants.hpp:1089
mfem::InvariantsEvaluator3D::Assemble_ddI3b
void Assemble_ddI3b(scalar_t w, scalar_t *A)
Definition invariants.hpp:1252
mfem::InvariantsEvaluator3D::Assemble_ddI1b
void Assemble_ddI1b(scalar_t w, scalar_t *A)
Definition invariants.hpp:826
b
real_t b
Definition lissajous.cpp:42
a
real_t a
Definition lissajous.cpp:41
mfem::ScalarOps
Auxiliary class used as the default for the second template parameter in the classes InvariantsEvalua...
Definition invariants.hpp:27
mfem::ScalarOps::pow
static scalar_t pow(const scalar_t &x, int m, int n)
Definition invariants.hpp:31
mfem::ScalarOps::sign
static scalar_t sign(const scalar_t &a)
Definition invariants.hpp:28
Generated on Thu Dec 11 2025 11:52:55 for MFEM by doxygen 1.11.0