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