b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2element_continuum_util.H
1//------------------------------------------------------------------------
2// b2element_continuum_util.H --
3//
4// Utility classes and functions for 2D and 3D continuum elements.
5//
6// written by Thomas Ludwig
7//
8// Copyright (c) 2007-2012,2016
9// SMR Engineering & Development SA
10// 2502 Bienne, Switzerland
11//
12// All Rights Reserved. Proprietary source code. The contents of
13// this file may not be disclosed to third parties, copied or
14// duplicated in any form, in whole or in part, without the prior
15// written permission of SMR.
16//
17// CONTENTS:
18//
19// - A few type definitions for commonly used linear_algebra classes.
20//
21// - There are a few utility functions that are often used for tensor
22// manipulation etc.
23//------------------------------------------------------------------------
24
25#ifndef __B2ELEMENT_CONTINUUM_UTIL_H__
26#define __B2ELEMENT_CONTINUUM_UTIL_H__
27
28#include <cassert>
29#include <map>
30#include <vector>
31
32#include "elements/properties/b2solid_mechanics.H"
33#include "model/b2element.H"
34#include "model/b2solution.H"
35#include "utils/b2linear_algebra.H"
36
37// #define TEST
38
39namespace b2000 {
40
41//----------------------------------------------------------------------
42// Shorthands for vector and matrix types.
43//----------------------------------------------------------------------
44
45typedef b2linalg::Vector<double, b2linalg::Vdense> VD;
46typedef b2linalg::Vector<double, b2linalg::Vdense_ref> VR;
47typedef b2linalg::Vector<double, b2linalg::Vdense_constref> VDC;
48typedef b2linalg::Matrix<double, b2linalg::Mrectangle> ME;
49typedef b2linalg::Matrix<double, b2linalg::Mrectangle_constref> MEC;
50typedef b2linalg::Matrix<double, b2linalg::Mpacked> MP;
51typedef b2linalg::Matrix<double, b2linalg::Mupper_packed> MUP;
52typedef b2linalg::Matrix<double, b2linalg::Mupper_packed_constref> MUPC;
53typedef b2linalg::Matrix<double, b2linalg::Mlower_packed> MLP;
54typedef b2linalg::Matrix<double, b2linalg::Mlower_packed_constref> MLPC;
55
56//----------------------------------------------------------------------
57// Utility functions.
58//----------------------------------------------------------------------
59
60namespace ffmv {
61
62template <int SIZE>
63inline void sc_eq_va_mul_vb(double& c, const double A[SIZE], const double B[SIZE]) {
64 c = 0.0;
65 for (int i = 0; i != SIZE; ++i) { c += A[i] * B[i]; }
66}
67
68template <int SIZE_BC, int SIZE_AB>
69inline void vc_eq_va_mul_mb(
70 double C[SIZE_BC], const double A[SIZE_AB], const double B[SIZE_AB][SIZE_BC]) {
71 for (int i = 0; i != SIZE_BC; ++i) {
72 double s = 0.0;
73 for (int j = 0; j != SIZE_AB; ++j) { s += A[j] * B[j][i]; }
74 C[i] = s;
75 }
76}
77
78template <int SIZE_BC, int SIZE_AB>
79inline void vc_eq_va_mul_mbt(
80 double C[SIZE_BC], const double A[SIZE_AB], const double B[SIZE_BC][SIZE_AB]) {
81 for (int i = 0; i != SIZE_BC; ++i) {
82 double s = 0.0;
83 for (int j = 0; j != SIZE_AB; ++j) { s += A[j] * B[i][j]; }
84 C[i] = s;
85 }
86}
87
88template <int SIZE_CR, int SIZE_CC>
89inline void mb_eq_mat(double C[SIZE_CR][SIZE_CC], const double A[SIZE_CC][SIZE_CR]) {
90 for (int i = 0; i != SIZE_CR; ++i) {
91 for (int j = 0; j != SIZE_CC; ++j) { C[i][j] = A[j][i]; }
92 }
93}
94
95template <int SIZE_CR, int SIZE_CC>
96inline void mb_eq_add_mat(double C[SIZE_CR][SIZE_CC], const double A[SIZE_CC][SIZE_CR]) {
97 for (int i = 0; i != SIZE_CR; ++i) {
98 for (int j = 0; j != SIZE_CC; ++j) { C[i][j] += A[j][i]; }
99 }
100}
101
102template <int SIZE_CR, int SIZE_CC>
103inline void mc_eq_ma_mul_sb(
104 double C[SIZE_CR][SIZE_CC], const double A[SIZE_CR][SIZE_CC], const double sb) {
105 for (int j = 0; j != SIZE_CR; ++j) {
106 for (int i = 0; i != SIZE_CC; ++i) { C[j][i] = sb * A[j][i]; }
107 }
108}
109
110template <int SIZE_CR, int SIZE_CC>
111inline void mc_eq_add_ma_mul_sb(
112 double C[SIZE_CR][SIZE_CC], const double A[SIZE_CR][SIZE_CC], const double sb) {
113 double* Cp = C[0];
114 const double* Ap = A[0];
115 for (int i = 0; i != SIZE_CR * SIZE_CC; ++i) { Cp[i] += sb * Ap[i]; }
116}
117
118template <int SIZE_CR, int SIZE_CC>
119inline void mc_eq_ma_plus_mb(
120 double C[SIZE_CR][SIZE_CC], const double A[SIZE_CR][SIZE_CC],
121 const double B[SIZE_CR][SIZE_CC]) {
122 double* Cp = C[0];
123 const double* Ap = A[0];
124 const double* Bp = B[0];
125 for (int i = 0; i != SIZE_CR * SIZE_CC; ++i) { Cp[i] = Ap[i] + Bp[i]; }
126}
127
128template <int SIZE_CR, int SIZE_CC>
129inline void mc_eq_add_ma_plus_mb(
130 double C[SIZE_CR][SIZE_CC], const double A[SIZE_CR][SIZE_CC],
131 const double B[SIZE_CR][SIZE_CC]) {
132 double* Cp = C[0];
133 const double* Ap = A[0];
134 const double* Bp = B[0];
135 for (int i = 0; i != SIZE_CR * SIZE_CC; ++i) { Cp[i] += Ap[i] + Bp[i]; }
136}
137
138template <int SIZE_CR, int SIZE_CC>
139inline void mc_eq_ma_minus_mb(
140 double C[SIZE_CR][SIZE_CC], const double A[SIZE_CR][SIZE_CC],
141 const double B[SIZE_CR][SIZE_CC]) {
142 double* Cp = C[0];
143 const double* Ap = A[0];
144 const double* Bp = B[0];
145 for (int i = 0; i != SIZE_CR * SIZE_CC; ++i) { Cp[i] = Ap[i] - Bp[i]; }
146}
147
148template <int SIZE_CR, int SIZE_CC>
149inline void mc_eq_add_ma_minus_mb(
150 double C[SIZE_CR][SIZE_CC], const double A[SIZE_CR][SIZE_CC],
151 const double B[SIZE_CR][SIZE_CC]) {
152 double* Cp = C[0];
153 const double* Ap = A[0];
154 const double* Bp = B[0];
155 for (int i = 0; i != SIZE_CR * SIZE_CC; ++i) { Cp[i] += Ap[i] - Bp[i]; }
156}
157
158template <int SIZE_CR, int SIZE_CC, int SIZE_AB>
159inline void mc_eq_ma_mul_mb(
160 double C[SIZE_CR][SIZE_CC], const double A[SIZE_CR][SIZE_AB],
161 const double B[SIZE_AB][SIZE_CC]) {
162 for (int i = 0; i != SIZE_CR; ++i) {
163 for (int j = 0; j != SIZE_CC; ++j) {
164 double s = 0.0;
165 for (int k = 0; k != SIZE_AB; ++k) { s += A[i][k] * B[k][j]; }
166 C[i][j] = s;
167 }
168 }
169}
170
171template <int SIZE_CR, int SIZE_CC, int SIZE_AB>
172inline void mc_eq_add_ma_mul_mb(
173 double C[SIZE_CR][SIZE_CC], const double A[SIZE_CR][SIZE_AB],
174 const double B[SIZE_AB][SIZE_CC]) {
175 for (int i = 0; i != SIZE_CR; ++i) {
176 for (int j = 0; j != SIZE_CC; ++j) {
177 double s = 0.0;
178 for (int k = 0; k != SIZE_AB; ++k) { s += A[i][k] * B[k][j]; }
179 C[i][j] += s;
180 }
181 }
182}
183
184template <int SIZE_CR, int SIZE_CC, int SIZE_AB>
185inline void mc_eq_mat_mul_mb(
186 double C[SIZE_CR][SIZE_CC], const double AT[SIZE_AB][SIZE_CR],
187 const double B[SIZE_AB][SIZE_CC]) {
188 for (int i = 0; i != SIZE_CR; ++i) {
189 for (int j = 0; j != SIZE_CC; ++j) {
190 double s = 0.0;
191 for (int k = 0; k != SIZE_AB; ++k) { s += AT[k][i] * B[k][j]; }
192 C[i][j] = s;
193 }
194 }
195}
196
197template <int SIZE_CR, int SIZE_CC, int SIZE_AB>
198inline void mc_eq_add_mat_mul_mb(
199 double C[SIZE_CR][SIZE_CC], const double AT[SIZE_AB][SIZE_CR],
200 const double B[SIZE_AB][SIZE_CC]) {
201 for (int i = 0; i != SIZE_CR; ++i) {
202 for (int j = 0; j != SIZE_CC; ++j) {
203 double s = 0.0;
204 for (int k = 0; k != SIZE_AB; ++k) { s += AT[k][i] * B[k][j]; }
205 C[i][j] += s;
206 }
207 }
208}
209
210template <int SIZE_CR, int SIZE_CC, int SIZE_AB>
211inline void mc_eq_ma_mul_mbt(
212 double C[SIZE_CR][SIZE_CC], const double A[SIZE_CR][SIZE_AB],
213 const double BT[SIZE_CC][SIZE_AB]) {
214 for (int i = 0; i != SIZE_CR; ++i) {
215 for (int j = 0; j != SIZE_CC; ++j) {
216 double s = 0.0;
217 for (int k = 0; k != SIZE_AB; ++k) { s += A[i][k] * BT[j][k]; }
218 C[i][j] = s;
219 }
220 }
221}
222
223template <int SIZE_CR, int SIZE_CC, int SIZE_AB>
224inline void mc_eq_add_ma_mul_mbt(
225 double C[SIZE_CR][SIZE_CC], const double A[SIZE_CR][SIZE_AB],
226 const double BT[SIZE_CC][SIZE_AB]) {
227 for (int i = 0; i != SIZE_CR; ++i) {
228 for (int j = 0; j != SIZE_CC; ++j) {
229 double s = 0.0;
230 for (int k = 0; k != SIZE_AB; ++k) { s += A[i][k] * BT[j][k]; }
231 C[i][j] += s;
232 }
233 }
234}
235
236template <int SIZE_CR, int SIZE_CC, int SIZE_AB>
237inline void mc_eq_mat_mul_mbt(
238 double C[SIZE_CR][SIZE_CC], const double AT[SIZE_AB][SIZE_CR],
239 const double BT[SIZE_CC][SIZE_AB]) {
240 for (int i = 0; i != SIZE_CR; ++i) {
241 for (int j = 0; j != SIZE_CC; ++j) {
242 double s = 0.0;
243 for (int k = 0; k != SIZE_AB; ++k) { s += AT[k][i] * BT[j][k]; }
244 C[i][j] = s;
245 }
246 }
247}
248
249template <int SIZE_CR, int SIZE_CC, int SIZE_AB>
250inline void mc_eq_add_mat_mul_mbt(
251 double C[SIZE_CR][SIZE_CC], const double AT[SIZE_AB][SIZE_CR],
252 const double BT[SIZE_CC][SIZE_AB]) {
253 for (int i = 0; i != SIZE_CR; ++i) {
254 for (int j = 0; j != SIZE_CC; ++j) {
255 double s = 0.0;
256 for (int k = 0; k != SIZE_AB; ++k) { s += AT[k][i] * BT[j][k]; }
257 C[i][j] += s;
258 }
259 }
260}
261
262template <int SIZE_CR, int SIZE_CC, int SIZE_AB>
263inline void mct_eq_ma_mul_mb(
264 double CT[SIZE_CR][SIZE_CC], const double A[SIZE_CC][SIZE_AB],
265 const double B[SIZE_AB][SIZE_CR]) {
266 for (int j = 0; j != SIZE_CR; ++j) {
267 for (int i = 0; i != SIZE_CC; ++i) {
268 double s = 0.0;
269 for (int k = 0; k != SIZE_AB; ++k) { s += A[i][k] * B[k][j]; }
270 CT[j][i] = s;
271 }
272 }
273}
274
275template <int SIZE_CR, int SIZE_CC, int SIZE_AB>
276inline void mct_eq_add_ma_mul_mb(
277 double CT[SIZE_CR][SIZE_CC], const double A[SIZE_CC][SIZE_AB],
278 const double B[SIZE_AB][SIZE_CR]) {
279 for (int j = 0; j != SIZE_CR; ++j) {
280 for (int i = 0; i != SIZE_CC; ++i) {
281 double s = 0.0;
282 for (int k = 0; k != SIZE_AB; ++k) { s += A[i][k] * B[k][j]; }
283 CT[j][i] += s;
284 }
285 }
286}
287
288template <int SIZE_CR, int SIZE_CC, int SIZE_AB>
289inline void mct_eq_mat_mul_mb(
290 double CT[SIZE_CR][SIZE_CC], const double AT[SIZE_AB][SIZE_CC],
291 const double B[SIZE_AB][SIZE_CR]) {
292 for (int j = 0; j != SIZE_CR; ++j) {
293 for (int i = 0; i != SIZE_CC; ++i) {
294 double s = 0.0;
295 for (int k = 0; k != SIZE_AB; ++k) { s += AT[k][i] * B[k][j]; }
296 CT[j][i] = s;
297 }
298 }
299}
300
301template <int SIZE_CR, int SIZE_CC, int SIZE_AB>
302inline void mct_eq_add_mat_mul_mb(
303 double CT[SIZE_CR][SIZE_CC], const double AT[SIZE_AB][SIZE_CC],
304 const double B[SIZE_AB][SIZE_CR]) {
305 for (int j = 0; j != SIZE_CR; ++j) {
306 for (int i = 0; i != SIZE_CC; ++i) {
307 double s = 0.0;
308 for (int k = 0; k != SIZE_AB; ++k) { s += AT[k][i] * B[k][j]; }
309 CT[j][i] += s;
310 }
311 }
312}
313
314template <int SIZE_CR, int SIZE_CC, int SIZE_AB>
315inline void mct_eq_ma_mul_mbt(
316 double CT[SIZE_CR][SIZE_CC], const double A[SIZE_CC][SIZE_AB],
317 const double BT[SIZE_CR][SIZE_AB]) {
318 for (int j = 0; j != SIZE_CR; ++j) {
319 for (int i = 0; i != SIZE_CC; ++i) {
320 double s = 0.0;
321 for (int k = 0; k != SIZE_AB; ++k) { s += A[i][k] * BT[j][k]; }
322 CT[j][i] = s;
323 }
324 }
325}
326
327template <int SIZE_CR, int SIZE_CC, int SIZE_AB>
328inline void mct_eq_add_ma_mul_mbt(
329 double CT[SIZE_CR][SIZE_CC], const double A[SIZE_CC][SIZE_AB],
330 const double BT[SIZE_CR][SIZE_AB]) {
331 for (int j = 0; j != SIZE_CR; ++j) {
332 for (int i = 0; i != SIZE_CC; ++i) {
333 double s = 0.0;
334 for (int k = 0; k != SIZE_AB; ++k) { s += A[i][k] * BT[j][k]; }
335 CT[j][i] += s;
336 }
337 }
338}
339
340template <int SIZE_CR, int SIZE_CC, int SIZE_AB>
341inline void mct_eq_mat_mul_mbt(
342 double CT[SIZE_CR][SIZE_CC], const double AT[SIZE_AB][SIZE_CC],
343 const double BT[SIZE_CR][SIZE_AB]) {
344 for (int j = 0; j != SIZE_CR; ++j) {
345 for (int i = 0; i != SIZE_CC; ++i) {
346 double s = 0.0;
347 for (int k = 0; k != SIZE_AB; ++k) { s += AT[k][i] * BT[j][k]; }
348 CT[j][i] = s;
349 }
350 }
351}
352
353template <int SIZE_CR, int SIZE_CC, int SIZE_AB>
354inline void mct_eq_add_mat_mul_mbt(
355 double CT[SIZE_CR][SIZE_CC], const double AT[SIZE_AB][SIZE_CC],
356 const double BT[SIZE_CR][SIZE_AB]) {
357 for (int j = 0; j != SIZE_CR; ++j) {
358 for (int i = 0; i != SIZE_CC; ++i) {
359 double s = 0.0;
360 for (int k = 0; k != SIZE_AB; ++k) { s += AT[k][i] * BT[j][k]; }
361 CT[j][i] += s;
362 }
363 }
364}
365
366} // namespace ffmv
367
368namespace fmv {
369
370template <int SIZE>
371inline void sc_eq_va_mul_vb(double& c, const VDC& A, const VDC& B) {
372 c = 0.0;
373 for (int i = 0; i < SIZE; ++i) { c += A[i] * B[i]; }
374}
375
376template <int SIZE_CR, int SIZE_CC>
377inline void mc_eq_ma_mul_sb(ME& C, const MEC& A, const double sb) {
378 for (int i = 0; i < SIZE_CR; ++i) {
379 for (int j = 0; j < SIZE_CC; ++j) { C(i, j) = sb * A(i, j); }
380 }
381}
382
383template <int SIZE_CR, int SIZE_CC>
384inline void mc_eq_ma_plus_mb(ME& C, const MEC& A, const MEC& B) {
385 for (int i = 0; i < SIZE_CR; ++i) {
386 for (int j = 0; j < SIZE_CC; ++j) { C(i, j) = A(i, j) + B(i, j); }
387 }
388}
389
390template <int SIZE_CR, int SIZE_CC>
391inline void mc_eq_ma_minus_mb(ME& C, const MEC& A, const MEC& B) {
392 for (int i = 0; i < SIZE_CR; ++i) {
393 for (int j = 0; j < SIZE_CC; ++j) { C(i, j) = A(i, j) - B(i, j); }
394 }
395}
396
397template <int SIZE_CR, int SIZE_CC, int SIZE_AB>
398inline void mc_eq_ma_mul_mb(ME& C, const MEC& A, const MEC& B) {
399 for (int i = 0; i < SIZE_CR; ++i) {
400 for (int j = 0; j < SIZE_CC; ++j) {
401 double s = 0.0;
402 for (int k = 0; k < SIZE_AB; ++k) { s += A(i, k) * B(k, j); }
403 C(i, j) = s;
404 }
405 }
406}
407
408template <int SIZE_CR, int SIZE_CC, int SIZE_AB>
409inline void mc_eq_mat_mul_mb(ME& C, const MEC& A, const MEC& B) {
410 for (int i = 0; i < SIZE_CR; ++i) {
411 for (int j = 0; j < SIZE_CC; ++j) {
412 double s = 0.0;
413 for (int k = 0; k < SIZE_AB; ++k) { s += A(k, i) * B(k, j); }
414 C(i, j) = s;
415 }
416 }
417}
418
419template <int SIZE_CR, int SIZE_CC, int SIZE_AB>
420inline void mc_eq_ma_mul_mbt(ME& C, const MEC& A, const MEC& B) {
421 for (int i = 0; i < SIZE_CR; ++i) {
422 for (int j = 0; j < SIZE_CC; ++j) {
423 double s = 0.0;
424 for (int k = 0; k < SIZE_AB; ++k) { s += A(i, k) * B(j, k); }
425 C(i, j) = s;
426 }
427 }
428}
429
430template <int SIZE_CR, int SIZE_CC, int SIZE_AB>
431inline void mc_eq_mat_mul_mbt(ME& C, const MEC& A, const MEC& B) {
432 for (int i = 0; i < SIZE_CR; ++i) {
433 for (int j = 0; j < SIZE_CC; ++j) {
434 double s = 0.0;
435 for (int k = 0; k < SIZE_AB; ++k) { s += A(k, i) * B(j, k); }
436 C(i, j) = s;
437 }
438 }
439}
440
441} // namespace fmv
442
443//----------------------------------------------------------------------
444// eval_determinant_2x2 --
445//
446// Compute the determinant of a 2x2 matrix.
447//
448// Return value:
449// Returns the determinant of the input matrix.
450//
451// Parameters:
452// a(in) : Input matrix.
453//----------------------------------------------------------------------
454inline double eval_determinant_2x2(const ME& a) {
455 // Get the determinant of the matrix.
456 double det = a[0][0] * a[1][1] - a[0][1] * a[1][0];
457 return det;
458}
459
460//----------------------------------------------------------------------
461// eval_determinant_and_adjoint_2x2 --
462//
463// Compute the determinant and the adjoint of a 2x2 matrix.
464// The inverse can then be computed inv(A) = adj(A) / det(A).
465//
466// Return value:
467// Returns the determinant of the input matrix.
468//
469// Parameters:
470// a(in) : Input matrix.
471// adj(out): On enter, must be a 2x2 matrix. On return,
472// contains the adjoint of the input matrix.
473//----------------------------------------------------------------------
474inline double eval_determinant_and_adjoint_2x2(const ME& a, ME& adj) {
475 // Get the determinant of the matrix.
476 double det = eval_determinant_2x2(a);
477
478 // Get the adjoint of the matrix.
479 adj[0][0] = +a[1][1];
480 adj[0][1] = -a[0][1];
481 adj[1][0] = -a[1][0];
482 adj[1][1] = +a[0][0];
483
484 return det;
485}
486
487//----------------------------------------------------------------------
488// eval_determinant_3x3 --
489//
490// Compute the determinant of a 3x3 matrix.
491//
492// Return value:
493// Returns the determinant of the input matrix.
494//
495// Parameters:
496// a(in) : Input matrix.
497//----------------------------------------------------------------------
498inline double eval_determinant_3x3(const ME& a) {
499 // Get the determinant of the matrix.
500 double det = a[0][0] * a[1][1] * a[2][2] + a[0][1] * a[1][2] * a[2][0]
501 + a[0][2] * a[1][0] * a[2][1] - a[0][0] * a[1][2] * a[2][1]
502 - a[0][1] * a[1][0] * a[2][2] - a[0][2] * a[1][1] * a[2][0];
503 return det;
504}
505
506//----------------------------------------------------------------------
507// eval_determinant_and_adjoint_3x3 --
508//
509// Compute the determinant and the adjoint of a 3x3 matrix.
510// The inverse can then be computed inv(A) = adj(A) / det(A).
511//
512// Return value:
513// Returns the determinant of the input matrix.
514//
515// Parameters:
516// a(in) : Input matrix.
517// adj(out): On enter, must be a 3x3 matrix. On return,
518// contains the adjoint of the input matrix.
519//
520// Maxima commands:
521// (%i1) a: matrix([a00,a01,a02],[a10,a11,a12],[a20,a21,a22]);
522// (%i2) adjoint(a);
523//----------------------------------------------------------------------
524inline double eval_determinant_and_adjoint_3x3(const ME& a, ME& adj) {
525 double a00 = a[0][0];
526 double a01 = a[0][1];
527 double a02 = a[0][2];
528 double a10 = a[1][0];
529 double a11 = a[1][1];
530 double a12 = a[1][2];
531 double a20 = a[2][0];
532 double a21 = a[2][1];
533 double a22 = a[2][2];
534
535 // Get the adjoint of the matrix.
536 adj[0][0] = a11 * a22 - a12 * a21;
537 adj[0][1] = a02 * a21 - a01 * a22;
538 adj[0][2] = a01 * a12 - a02 * a11;
539 adj[1][0] = a12 * a20 - a10 * a22;
540 adj[1][1] = a00 * a22 - a02 * a20;
541 adj[1][2] = a02 * a10 - a00 * a12;
542 adj[2][0] = a10 * a21 - a11 * a20;
543 adj[2][1] = a01 * a20 - a00 * a21;
544 adj[2][2] = a00 * a11 - a01 * a10;
545
546 // Return the determinant of the matrix.
547 return (
548 a00 * a11 * a22 + a01 * a12 * a20 + a02 * a10 * a21 - a00 * a12 * a21 - a01 * a10 * a22
549 - a02 * a11 * a20);
550}
551
552//----------------------------------------------------------------------
553// eval_linear_strain_3 --
554//
555// Compute the linear strain as a 6-vector from the 3x3
556// matrix F containing the deformation gradient tensor.
557//
558// E = 0.5 * (F^T + F) - I
559//
560// The 6-vector follows the Bathe conventions and contains:
561// (00, 11, 22, 01+10, 12+21, 02+20)
562//
563// Parameters:
564// F(in) : Deformation gradient tensor, 3x3 matrix.
565// E(out) : Array of 6 doubles, contains the linear
566// strain tensor.
567//----------------------------------------------------------------------
568inline void eval_linear_strain_3(const ME& F, double* E) {
569 const double f00 = F(0, 0);
570 const double f01 = F(0, 1);
571 const double f02 = F(0, 2);
572 const double f10 = F(1, 0);
573 const double f11 = F(1, 1);
574 const double f12 = F(1, 2);
575 const double f20 = F(2, 0);
576 const double f21 = F(2, 1);
577 const double f22 = F(2, 2);
578
579 E[0] = f00 - 1.0;
580 E[1] = f11 - 1.0;
581 E[2] = f22 - 1.0;
582 E[3] = f01 + f10;
583 E[4] = f12 + f21;
584 E[5] = f02 + f20;
585}
586
587inline void eval_linear_strain_3(const double F[3][3], double E[6]) {
588 const double f00 = F[0][0];
589 const double f01 = F[0][1];
590 const double f02 = F[0][2];
591 const double f10 = F[1][0];
592 const double f11 = F[1][1];
593 const double f12 = F[1][2];
594 const double f20 = F[2][0];
595 const double f21 = F[2][1];
596 const double f22 = F[2][2];
597
598 E[0] = f00 - 1.0;
599 E[1] = f11 - 1.0;
600 E[2] = f22 - 1.0;
601 E[3] = f01 + f10;
602 E[4] = f12 + f21;
603 E[5] = f02 + f20;
604}
605
606//----------------------------------------------------------------------
607// eval_green_lagrange_strain_3 --
608//
609// Compute the Green-Lagrange strain as a 6-vector from the 3x3
610// matrix F containing the deformation gradient tensor.
611//
612// E = 0.5 * ((F^T * F) - I)
613//
614// The 6-vector follows the Bathe conventions and contains:
615// (00, 11, 22, 01+10, 12+21, 02+20)
616//
617// Parameters:
618// F(in) : Deformation gradient tensor, 3x3 matrix.
619// E(out) : Array of 6 doubles, contains the Green-Lagrange
620// strain tensor.
621//----------------------------------------------------------------------
622inline void eval_green_lagrange_strain_3(const ME& F, double* E) {
623 const double f00 = F(0, 0);
624 const double f01 = F(0, 1);
625 const double f02 = F(0, 2);
626 const double f10 = F(1, 0);
627 const double f11 = F(1, 1);
628 const double f12 = F(1, 2);
629 const double f20 = F(2, 0);
630 const double f21 = F(2, 1);
631 const double f22 = F(2, 2);
632
633 E[0] = 0.5 * (f00 * f00 + f10 * f10 + f20 * f20 - 1.0);
634 E[1] = 0.5 * (f01 * f01 + f11 * f11 + f21 * f21 - 1.0);
635 E[2] = 0.5 * (f02 * f02 + f12 * f12 + f22 * f22 - 1.0);
636 E[3] = f00 * f01 + f10 * f11 + f20 * f21;
637 E[4] = f01 * f02 + f11 * f12 + f21 * f22;
638 E[5] = f00 * f02 + f10 * f12 + f20 * f22;
639}
640
641inline void eval_green_lagrange_strain_3t(const double FT[3][3], double E[6]) {
642 const double f00 = FT[0][0];
643 const double f01 = FT[1][0];
644 const double f02 = FT[2][0];
645 const double f10 = FT[0][1];
646 const double f11 = FT[1][1];
647 const double f12 = FT[2][1];
648 const double f20 = FT[0][2];
649 const double f21 = FT[1][2];
650 const double f22 = FT[2][2];
651
652 E[0] = 0.5 * (f00 * f00 + f10 * f10 + f20 * f20 - 1.0);
653 E[1] = 0.5 * (f01 * f01 + f11 * f11 + f21 * f21 - 1.0);
654 E[2] = 0.5 * (f02 * f02 + f12 * f12 + f22 * f22 - 1.0);
655 E[3] = f00 * f01 + f10 * f11 + f20 * f21;
656 E[4] = f01 * f02 + f11 * f12 + f21 * f22;
657 E[5] = f00 * f02 + f10 * f12 + f20 * f22;
658}
659
660//----------------------------------------------------------------------
661// eval_linear_strain_2 --
662//
663// Compute the linear strain as a 3-vector from the 2x2
664// matrix F containing the deformation gradient tensor.
665//
666// E = 0.5 * (F^T + F) - I
667//
668// The 3-vector follows the Bathe conventions and contains:
669// (00, 11, 01+10)
670//
671// Parameters:
672// F(in) : Deformation gradient tensor, 2x2 matrix.
673// E(out) : Array of 3 doubles, contains the linear
674// strain tensor.
675//----------------------------------------------------------------------
676inline void eval_linear_strain_2(const ME& F, double* E) {
677 const double f00 = F(0, 0);
678 const double f01 = F(0, 1);
679 const double f10 = F(1, 0);
680 const double f11 = F(1, 1);
681
682 E[0] = f00 - 1.0;
683 E[1] = f11 - 1.0;
684 E[2] = f01 + f10;
685}
686
687inline void eval_linear_strain_2(const double F[2][2], double E[3]) {
688 const double f00 = F[0][0];
689 const double f01 = F[0][1];
690 const double f10 = F[1][0];
691 const double f11 = F[1][1];
692
693 E[0] = f00 - 1.0;
694 E[1] = f11 - 1.0;
695 E[2] = f01 + f10;
696}
697
698//----------------------------------------------------------------------
699// eval_green_lagrange_strain_2 --
700//
701// Compute the Green-Lagrange strain as a 3-vector from the 2x2
702// matrix F containing the deformation gradient tensor.
703//
704// E = 0.5 * ((F^T * F) - I)
705//
706// The 6-vector follows the Bathe conventions and contains:
707// (00, 11, 01+10)
708//
709// Parameters:
710// F(in) : Deformation gradient tensor, 2x2 matrix.
711// E(out) : Array of 3 doubles, contains the Green-Lagrange
712// strain tensor.
713//----------------------------------------------------------------------
714inline void eval_green_lagrange_strain_2(const ME& F, double* E) {
715 const double f00 = F(0, 0);
716 const double f01 = F(0, 1);
717 const double f10 = F(1, 0);
718 const double f11 = F(1, 1);
719
720 E[0] = 0.5 * (f00 * f00 + f10 * f10 - 1.0);
721 E[1] = 0.5 * (f01 * f01 + f11 * f11 - 1.0);
722 E[2] = f00 * f01 + f10 * f11;
723}
724
725inline void eval_green_lagrange_strain_2t(const double FT[2][2], double E[3]) {
726 const double f00 = FT[0][0];
727 const double f01 = FT[1][0];
728 const double f10 = FT[0][1];
729 const double f11 = FT[1][1];
730
731 E[0] = 0.5 * (f00 * f00 + f10 * f10 - 1.0);
732 E[1] = 0.5 * (f01 * f01 + f11 * f11 - 1.0);
733 E[2] = f00 * f01 + f10 * f11;
734}
735
736// Maps a symmetric 2nd-order tensor to a 6-vector. The order of
737// the tensor components in the 6-vector is (00, 11, 22, 01/10,
738// 12/21, 02/20). This table is used for lookup into the
739// symmetric stress-strain matrix as returned by the property
740// object.
741const int tensor_3x3_to_6_conversion_table[3][3] = {{0, 3, 5}, {3, 1, 4}, {5, 4, 2}};
742
743void tensor_strain_3x3_to_6(const ME& a, double* b);
744
745void tensor_stress_6_to_3x3(const double* a, ME& b);
746
747void tensor_stress_3x3_to_6(const ME& a, double* b);
748
749void tensor_strain_2x2_to_3(const ME& a, double* b);
750
751void tensor_stress_3_to_2x2(const double* a, ME& b);
752
753void make_director_3x3(ME& m);
754
755double get_surface(ME& m);
756
757// From a Jacobian matrix (2x3 or 3x3), get the surface directors
758// (the first two lines), compute the surface spanned by them, and
759// normalize them to unit length.
760void get_surface_directors(const ME& J, ME& d_coor_d_lcoor, double& surface);
761
762//----------------------------------------------------------------------
763// This class describes natural coordinates. It has an ordering
764// property, hence it can be used as key for STL containers.
765//----------------------------------------------------------------------
766class Natural_Coor {
767public:
768 inline Natural_Coor() {
769 rst[0] = 0;
770 rst[1] = 0;
771 rst[2] = 0;
772 }
773
774 inline Natural_Coor(double r_, double s_, double t_) {
775 rst[0] = r_;
776 rst[1] = s_;
777 rst[2] = t_;
778 }
779
780 inline Natural_Coor(const double* rst_) {
781 rst[0] = rst_[0];
782 rst[1] = rst_[1];
783 rst[2] = rst_[2];
784 }
785
786 inline const double& operator[](int index) const {
787 assert(index >= 0 && index < 3);
788 return rst[index];
789 }
790
791 inline double& operator[](int index) {
792 assert(index >= 0 && index < 3);
793 return rst[index];
794 }
795
796 inline const double* data() const { return rst; }
797
798 inline double* data() { return rst; }
799
800 inline bool operator<(const Natural_Coor& n) const {
801 if (rst[2] != n.rst[2]) { return rst[2] < n.rst[2]; }
802 if (rst[1] != n.rst[1]) { return rst[1] < n.rst[1]; }
803 return rst[0] < n.rst[0];
804 }
805
806private:
807 double rst[3];
808};
809
810//----------------------------------------------------------------------
811// Class to hold values for shape functions and their derivatives for
812// given natural coordinates. Objects of this class store the
813// results obtained from the Shape_Evaluator objects.
814//----------------------------------------------------------------------
815class Shape_Value {
816public:
817 Shape_Value(int num_nodes) : h(num_nodes), h_derived(3, num_nodes) { ; }
818
819 VD h;
820 ME h_derived;
821};
822
823//----------------------------------------------------------------------
824// The (pre-)evaluated shape functions and derivatives are stored in
825// an STL map. These maps are static per element-type and per
826// processing node (the latter is not yet implemented, but needed for
827// parallelization).
828//
829// Template parameters:
830// - SE : Shape evaluator.
831//----------------------------------------------------------------------
832template <typename SE>
833class Shape_Value_Map : public std::map<Natural_Coor, Shape_Value> {
834public:
835 Shape_Value_Map(){};
836 virtual ~Shape_Value_Map(){};
837
838 iterator insert(const key_type& natural_coor) {
839 iterator i = find(natural_coor);
840 // if (i != end())
841 // std::cerr << "i != end " << (*i).first[0]
842 // << "," << (*i).first[1]
843 // << "," << (*i).first[2] << std::endl;
844 if (i == end()) {
845 Shape_Value shape_value(SE::num_nodes);
846 SE::eval_h(natural_coor.data(), shape_value.h);
847 SE::eval_h_derived(natural_coor.data(), shape_value.h_derived);
848 std::pair<iterator, bool> result = std::map<Natural_Coor, Shape_Value>::insert(
849 value_type(natural_coor, shape_value));
850 i = result.first;
851 }
852 return i;
853 }
854};
855
856typedef std::map<Natural_Coor, Shape_Value>::const_iterator svm_const_iterator_t;
857
858//----------------------------------------------------------------------
859// This class uniquely identifies an integration point by its
860// natural coordinates and the layer. It has an ordering
861// property, hence it can be used as key for STL containers.
862//
863// The inclusion of the layer-id is due to the fact that certain
864// integration schemes define integration points at the border
865// between two layers, hence for these schemes, there are
866// duplicate integration points. These must be distinguished
867// (different material properties).
868//----------------------------------------------------------------------
869class IP_ID : public Natural_Coor {
870public:
871 inline IP_ID() : f_layer_id(0) {}
872
873 inline IP_ID(double r_, double s_, double t_, int layer_id_)
874 : Natural_Coor(r_, s_, t_), f_layer_id(layer_id_) {
875 ;
876 }
877
878 inline IP_ID(const double* rst_, int layer_id_) : Natural_Coor(rst_), f_layer_id(layer_id_) {
879 ;
880 }
881
882 inline IP_ID(const Natural_Coor& n, int layer_id_) : Natural_Coor(n), f_layer_id(layer_id_) {
883 ;
884 }
885
886 inline int layer_id() const { return f_layer_id; }
887
888 inline int& layer_id() { return f_layer_id; }
889
890 inline bool operator<(const IP_ID& i) const {
891 if (f_layer_id != i.f_layer_id) { return f_layer_id < i.f_layer_id; }
892 if ((*this)[2] != i[2]) { return (*this)[2] < i[2]; }
893 if ((*this)[1] != i[1]) { return (*this)[1] < i[1]; }
894 return (*this)[0] < i[0];
895 }
896
897private:
898 int f_layer_id; // The id of the layer.
899};
900
901//----------------------------------------------------------------------
902// add_w1_to_second_3d
903//
904// Description:
905//
906// This routine adds 'volume * B^T * C * B' to the second
907// variation. The loops are arranged such that the number of
908// floating-point operations is minimized.
909//
910// Remarks:
911//
912// - The matrix B must be of size (6, NDOF).
913//
914// - The stiffness matrix may be in lower- or upper-packed form and
915// must be of size (NDOF, NDOF).
916//
917// - The stiffness matrix may be actually larger than NDOF. This
918// is useful for elements having non-displacement dofs. However
919// it is required that the displacement dofs are at position
920// [0..NDOF-1].
921//
922// - The array C must correspond to a lower-packed 6x6 matrix.
923//----------------------------------------------------------------------
924inline void add_w1_to_second_3d(
925 const ME& B, // Strain-displ. matrix.
926 const double C[21], // Material matrix.
927 const double volume,
928 MP& d_value_d_dof // Stiffness matrix.
929) {
930 assert(B.size1() == 6);
931 assert(d_value_d_dof.size1() == B.size2());
932 const size_t NDOF = B.size2();
933
934#ifdef TEST
935 ME test(NDOF, NDOF);
936 {
937 b2linalg::Matrix<double, b2linalg::Mlower_packed_constref> CC(6, C);
938 ME CCC(CC);
939 ME CB = CCC * B;
940 ME BTBC = transposed(B) * CB;
941 test = BTBC * volume;
942 ME KK(d_value_d_dof);
943 test += KK;
944 }
945#endif /* TEST */
946
947 for (size_t j = 0; j != NDOF; ++j) {
948 VDC Bj = B[j];
949
950 double t[6];
951 t[0] = Bj[0] * C[0] + Bj[1] * C[1] + Bj[2] * C[2] + Bj[3] * C[3] + Bj[4] * C[4]
952 + Bj[5] * C[5];
953
954 t[1] = Bj[0] * C[1] + Bj[1] * C[6] + Bj[2] * C[7] + Bj[3] * C[8] + Bj[4] * C[9]
955 + Bj[5] * C[10];
956
957 t[2] = Bj[0] * C[2] + Bj[1] * C[7] + Bj[2] * C[11] + Bj[3] * C[12] + Bj[4] * C[13]
958 + Bj[5] * C[14];
959
960 t[3] = Bj[0] * C[3] + Bj[1] * C[8] + Bj[2] * C[12] + Bj[3] * C[15] + Bj[4] * C[16]
961 + Bj[5] * C[17];
962
963 t[4] = Bj[0] * C[4] + Bj[1] * C[9] + Bj[2] * C[13] + Bj[3] * C[16] + Bj[4] * C[18]
964 + Bj[5] * C[19];
965
966 t[5] = Bj[0] * C[5] + Bj[1] * C[10] + Bj[2] * C[14] + Bj[3] * C[17] + Bj[4] * C[19]
967 + Bj[5] * C[20];
968
969 t[0] *= volume;
970 t[1] *= volume;
971 t[2] *= volume;
972 t[3] *= volume;
973 t[4] *= volume;
974 t[5] *= volume;
975
976 for (size_t i = 0; i <= j; ++i) {
977 VDC Bi = B[i];
978 const double s = Bi[0] * t[0] + Bi[1] * t[1] + Bi[2] * t[2] + Bi[3] * t[3]
979 + Bi[4] * t[4] + Bi[5] * t[5];
980 d_value_d_dof(i, j) += s;
981 }
982 }
983
984#ifdef TEST
985 {
986 double sd = 0.0;
987 for (size_t j = 0; j != NDOF; ++j) {
988 for (size_t i = 0; i <= j; ++i) { sd += std::abs(test(i, j) - d_value_d_dof(i, j)); }
989 }
990 if (sd > 1e-9) {
991 std::cerr << "Error in add_w1_to_second: sd=" << sd << std::endl;
992 for (size_t j = 0; j != NDOF; ++j) {
993 for (size_t i = 0; i <= j; ++i) {
994 const double diff = std::abs(test(i, j) - d_value_d_dof(i, j));
995 if (diff > 1e-11) {
996 std::cerr << std::setprecision(17) << " K[" << i << "," << j
997 << "]=" << d_value_d_dof(i, j) << " T=" << test(i, j)
998 << " diff=" << diff << std::endl;
999 }
1000 }
1001 }
1002 exit(1);
1003 }
1004 }
1005#endif /* TEST */
1006}
1007
1008//----------------------------------------------------------------------
1009// add_w1_to_second_3d
1010//
1011// Description:
1012//
1013// This routine adds 'volume * B^T * C * B' to the second
1014// variation. The loops are arranged such that the number of
1015// floating-point operations is minimized.
1016//
1017// Remarks:
1018//
1019// - NDOF is the number of degrees-of-freedom that are *used*.
1020// The matrix BT must be of size (NDOF, 6).
1021//
1022// - The stiffness matrix must be in upper-packed form. The
1023// internal storage format and access pattern is assumed.
1024//
1025// - The stiffness matrix may be actually larger than NDOF. This
1026// is useful for elements having non-displacement dofs. However
1027// it is required that the displacement dofs are at position
1028// [0..NDOF-1].
1029//
1030// - The array C must correspond to a lower-packed 6x6 matrix.
1031//----------------------------------------------------------------------
1032template <int NDOF>
1033inline void add_w1_to_second_3d(
1034 const double BT[NDOF][6], const double C[21], const double volume,
1035 MP& d_value_d_dof // Stiffness matrix.
1036) {
1037 double* second_0_0;
1038 double* second_0_j;
1039 const double* b_i;
1040 const double* b_j;
1041 double t[6];
1042 double s;
1043
1044 assert(d_value_d_dof.is_upper());
1045 second_0_0 = &d_value_d_dof(0, 0);
1046
1047 for (int j = 0; j != NDOF; ++j) {
1048 b_j = BT[j];
1049 t[0] = b_j[0] * C[0] + b_j[1] * C[1] + b_j[2] * C[2] + b_j[3] * C[3] + b_j[4] * C[4]
1050 + b_j[5] * C[5];
1051
1052 t[1] = b_j[0] * C[1] + b_j[1] * C[6] + b_j[2] * C[7] + b_j[3] * C[8] + b_j[4] * C[9]
1053 + b_j[5] * C[10];
1054
1055 t[2] = b_j[0] * C[2] + b_j[1] * C[7] + b_j[2] * C[11] + b_j[3] * C[12] + b_j[4] * C[13]
1056 + b_j[5] * C[14];
1057
1058 t[3] = b_j[0] * C[3] + b_j[1] * C[8] + b_j[2] * C[12] + b_j[3] * C[15] + b_j[4] * C[16]
1059 + b_j[5] * C[17];
1060
1061 t[4] = b_j[0] * C[4] + b_j[1] * C[9] + b_j[2] * C[13] + b_j[3] * C[16] + b_j[4] * C[18]
1062 + b_j[5] * C[19];
1063
1064 t[5] = b_j[0] * C[5] + b_j[1] * C[10] + b_j[2] * C[14] + b_j[3] * C[17] + b_j[4] * C[19]
1065 + b_j[5] * C[20];
1066
1067 t[0] *= volume;
1068 t[1] *= volume;
1069 t[2] *= volume;
1070 t[3] *= volume;
1071 t[4] *= volume;
1072 t[5] *= volume;
1073
1074 second_0_j = second_0_0 + (j * (j + 1) / 2);
1075
1076 for (int i = 0; i <= j; ++i) {
1077 b_i = BT[i];
1078 s = b_i[0] * t[0] + b_i[1] * t[1] + b_i[2] * t[2] + b_i[3] * t[3] + b_i[4] * t[4]
1079 + b_i[5] * t[5];
1080 second_0_j[i] += s;
1081 }
1082 }
1083}
1084
1085//----------------------------------------------------------------------
1086// add_w2_to_second_3d
1087//
1088// Description:
1089//
1090// This routine adds the "non-linear" stiffness contribution
1091// to the tangent stiffness (2nd variation) matrix.
1092//
1093// Remarks:
1094//
1095// - NDOF is the number of degrees-of-freedom that are *used*.
1096// The matrix D must be of size (3, NDOF / 3).
1097//
1098// - The stiffness matrix must be in upper-packed form. The
1099// internal storage format and access pattern is assumed.
1100//
1101// - The stiffness matrix may be actually larger than NDOF. This
1102// is useful for elements having non-displacement dofs. However
1103// it is required that the displacement dofs are at position
1104// [0..NDOF-1].
1105//
1106// Assumptions:
1107//
1108// The internal storage formats and access patterns for
1109// upper-packed and rectangular matrices is assumed.
1110//
1111// Complexity:
1112//
1113// Be n the number of nodes of the element. The number of
1114// floating-point additions is
1115//
1116// 3n^2 + 9n
1117//
1118// The number of floating-point multiplications is:
1119//
1120// 2n^2 + 11n
1121//
1122// The number of integer operations (increments, additions,
1123// multiplications, divisions) is
1124//
1125// 4n^2 + 19n
1126//
1127//----------------------------------------------------------------------
1128template <int NDOF>
1129inline void add_w2_to_second_3d(
1130 const ME& D, // Derivatives.
1131 const VD& S, // Stress tensor.
1132 const double volume,
1133 MP& d_value_d_dof // Stiffness matrix.
1134) {
1135 const double* d_i;
1136 const double* d_j;
1137 double* second_0_0;
1138 double* second_0_j0;
1139 double* second_0_j1;
1140 double* second_0_j2;
1141 double b[3];
1142 int i0;
1143 int i1;
1144 int i2;
1145 int j0;
1146 int j1;
1147 int j2;
1148
1149 assert(d_value_d_dof.is_upper());
1150 const int NNEL = NDOF / 3;
1151
1152#ifdef TEST
1153 MP test(NDOF, true);
1154 ME W2(NDOF, NDOF);
1155 W2.set_zero();
1156 {
1157 ME BNL(9, NDOF);
1158 BNL.set_zero();
1159 for (int k = 0; k < NNEL; k++) {
1160 int k0 = 3 * k + 0;
1161 int k1 = 3 * k + 1;
1162 int k2 = 3 * k + 2;
1163
1164 BNL(0, k0) = D(0, k);
1165 BNL(1, k0) = D(1, k);
1166 BNL(2, k0) = D(2, k);
1167 BNL(3, k1) = D(0, k);
1168 BNL(4, k1) = D(1, k);
1169 BNL(5, k1) = D(2, k);
1170 BNL(6, k2) = D(0, k);
1171 BNL(7, k2) = D(1, k);
1172 BNL(8, k2) = D(2, k);
1173 }
1174
1175 ME SM(9, 9);
1176 SM.set_zero();
1177 for (int i = 0; i < 3; i++) {
1178 int j = 3 * i;
1179 SM(j + 0, j + 0) = S[0]; // S(0, 0);
1180 SM(j + 0, j + 1) = S[3]; // S(0, 1);
1181 SM(j + 0, j + 2) = S[5]; // S(0, 2);
1182 SM(j + 1, j + 0) = S[3]; // S(1, 0);
1183 SM(j + 1, j + 1) = S[1]; // S(1, 1);
1184 SM(j + 1, j + 2) = S[4]; // S(1, 2);
1185 SM(j + 2, j + 0) = S[5]; // S(2, 0);
1186 SM(j + 2, j + 1) = S[4]; // S(2, 1);
1187 SM(j + 2, j + 2) = S[2]; // S(2, 2);
1188 }
1189 ME SMBNL(9, NDOF);
1190 SMBNL = SM * BNL;
1191 W2 = volume * transposed(BNL) * SMBNL;
1192
1193 test.set_zero();
1194 for (int j = 0; j < NDOF; j++) {
1195 for (int i = 0; i <= j; i++) { test(i, j) = W2(i, j); }
1196 }
1197 test += d_value_d_dof;
1198 }
1199#endif /* TEST */
1200
1201 second_0_0 = &d_value_d_dof(0, 0);
1202 for (int j = 0; j < NNEL; j++) {
1203 d_j = &D(0, j);
1204
1205 j0 = j * 3 + 0;
1206 j1 = j * 3 + 1;
1207 j2 = j * 3 + 2;
1208
1209 b[0] = S[0] * d_j[0] + S[3] * d_j[1] + S[5] * d_j[2];
1210 b[1] = S[3] * d_j[0] + S[1] * d_j[1] + S[4] * d_j[2];
1211 b[2] = S[5] * d_j[0] + S[4] * d_j[1] + S[2] * d_j[2];
1212
1213 second_0_j0 = second_0_0 + (j0 * (j0 + 1) / 2);
1214 second_0_j1 = second_0_0 + (j1 * (j1 + 1) / 2);
1215 second_0_j2 = second_0_0 + (j2 * (j2 + 1) / 2);
1216
1217 for (int i = 0; i <= j; i++) {
1218 d_i = &D(0, i);
1219
1220 double s = d_i[0] * b[0] + d_i[1] * b[1] + d_i[2] * b[2];
1221 s *= volume;
1222
1223 i0 = i * 3 + 0;
1224 i1 = i * 3 + 1;
1225 i2 = i * 3 + 2;
1226
1227 second_0_j0[i0] += s;
1228 second_0_j1[i1] += s;
1229 second_0_j2[i2] += s;
1230 }
1231 }
1232
1233#ifdef TEST
1234 double sd = 0.0;
1235 for (int j = 0; j < NDOF; j++) {
1236 for (int i = 0; i < NDOF; i++) {
1237 double diff = test(i, j) - d_value_d_dof(i, j);
1238 sd += diff * diff;
1239 }
1240 }
1241 if (sd > 1e-10) {
1242 std::cerr << "Error in add_w2_to_second: sd=" << sd << std::endl;
1243 exit(1);
1244 }
1245#endif /* TEST */
1246}
1247
1248//----------------------------------------------------------------------
1249// add_to_mass_matrix_3d
1250//
1251// Add the contribution of the integration point to
1252// the mass matrix.
1253//
1254// Remarks:
1255//
1256// - NDOF is the number of degrees-of-freedom that are *used*.
1257// The vector H must be of size NDOF / 3.
1258//
1259// - The mass matrix may be in upper- or lower-packed form.
1260//
1261// - The mass matrix may be actually larger than NDOF. This is
1262// useful for elements having non-displacement dofs. However
1263// it is required that the displacement dofs are at position
1264// [0..NDOF-1].
1265//
1266// - The mass matrix is computed for the reference configuration,
1267// thus it is assumed that the deformations are small.
1268//
1269//----------------------------------------------------------------------
1270inline void add_to_mass_matrix_3d(
1271 const VDC& H, // Shape functions.
1272 const double density, // Density at integr. point.
1273 const double volume, // Integration volume.
1274 MP& d_value_d_dofdotdot // Mass matrix.
1275) {
1276 const int NNEL = H.size();
1277 const int NDOF = NNEL * 3;
1278
1279#ifdef TEST
1280 ME HS(3, NDOF);
1281 HS.set_zero();
1282 for (int i = 0; i < 3; i++) {
1283 for (int j = 0; j < NNEL; j++) { HS(i, i + 3 * j) = H[j]; }
1284 }
1285
1286 MP test(NDOF, true);
1287 for (int kl = 0; kl < NDOF; kl++) {
1288 for (int mn = 0; mn <= kl; mn++) { test(kl, mn) = HS[kl] * HS[mn] * volume * density; }
1289 }
1290 test += d_value_d_dofdotdot;
1291#endif /* TEST */
1292
1293 const double dv = density * volume;
1294
1295 for (int kl = 0; kl < NDOF; kl++) {
1296 int k = kl / 3;
1297 int l = kl % 3;
1298 const double a = H[k] * dv;
1299 for (int m = 0; m <= k; ++m) {
1300 int mn = 3 * m + l;
1301 d_value_d_dofdotdot(kl, mn) += a * H[m];
1302 }
1303 }
1304
1305#ifdef TEST
1306 double sd = 0.0;
1307 double avg = 0.0;
1308 for (int j = 0; j < NDOF; j++) {
1309 for (int i = 0; i < NDOF; i++) {
1310 avg += test(i, j) * test(i, j);
1311 double diff = test(i, j) - d_value_d_dofdotdot(i, j);
1312 sd += diff * diff;
1313 }
1314 }
1315 avg /= (NDOF * (NDOF + 1));
1316 // std::cerr << "avg=" << avg << " sd=" << sd << std::endl;
1317 if (sd > 1e-10 * avg) {
1318 std::cerr << "Error in add_to_mass_matrix_3d: sd=" << sd << std::endl;
1319 exit(1);
1320 }
1321#endif /* TEST */
1322}
1323
1324//----------------------------------------------------------------------
1325// add_to_mass_matrix_3d
1326//
1327// Add the contribution of the integration point to
1328// the mass matrix.
1329//
1330// Remarks:
1331//
1332// - NDOF is the number of degrees-of-freedom that are *used*.
1333// The vector H must be of size NDOF / 3.
1334//
1335// - The mass matrix must be in upper-packed form. The internal
1336// storage format and access pattern is assumed.
1337//
1338// - The mass matrix may be actually larger than NDOF. This is
1339// useful for elements having non-displacement dofs. However
1340// it is required that the displacement dofs are at position
1341// [0..NDOF-1].
1342//
1343// - The mass matrix is computed for the reference configuration,
1344// thus it is assumed that the deformations are small.
1345//
1346//----------------------------------------------------------------------
1347template <int NNEL>
1348inline void add_to_mass_matrix_3d(
1349 const double H[NNEL], // Shape functions.
1350 const double density, // Density at integr. point.
1351 const double volume, // Integration volume.
1352 MP& d_value_d_dofdotdot // Mass matrix.
1353) {
1354 assert(d_value_d_dofdotdot.is_upper());
1355 const int NDOF = NNEL * 3;
1356
1357#ifdef TEST
1358 ME HS(3, NDOF);
1359 HS.set_zero();
1360 for (int i = 0; i < 3; i++) {
1361 for (int j = 0; j < NNEL; j++) { HS(i, i + 3 * j) = H[j]; }
1362 }
1363
1364 MP test(NDOF, true);
1365 for (int kl = 0; kl < NDOF; kl++) {
1366 for (int mn = 0; mn <= kl; mn++) { test(kl, mn) = HS[kl] * HS[mn] * volume * density; }
1367 }
1368 test += d_value_d_dofdotdot;
1369#endif /* TEST */
1370
1371 const double dv = density * volume;
1372
1373 double* mass_0_0 = &d_value_d_dofdotdot(0, 0);
1374 for (int kl = 0; kl != NDOF; ++kl) {
1375 const int k = kl / 3;
1376 const int l = kl % 3;
1377 const double a = H[k] * dv;
1378 double* mass_0_kl = mass_0_0 + (kl * (kl + 1) / 2);
1379 for (int m = 0; m <= k; ++m) {
1380 int mn = 3 * m + l;
1381 mass_0_kl[mn] += a * H[m];
1382 }
1383 }
1384
1385#ifdef TEST
1386 double sd = 0.0;
1387 double avg = 0.0;
1388 for (int j = 0; j < NDOF; j++) {
1389 for (int i = 0; i < NDOF; i++) {
1390 avg += test(i, j) * test(i, j);
1391 double diff = test(i, j) - d_value_d_dofdotdot(i, j);
1392 sd += diff * diff;
1393 }
1394 }
1395 avg /= (NDOF * (NDOF + 1));
1396 // std::cerr << "avg=" << avg << " sd=" << sd << std::endl;
1397 if (sd > 1e-10 * avg) {
1398 std::cerr << "Error in add_to_mass_matrix_3d: sd=" << sd << std::endl;
1399 exit(1);
1400 }
1401#endif /* TEST */
1402}
1403
1404//----------------------------------------------------------------------
1405// write_ip_coor_disp
1406//----------------------------------------------------------------------
1407inline void write_ip_coor_disp_3d(
1408 const double* xi, // Element-local coor.
1409 const int layer_id, // ID of material layer.
1410 const VDC& H, // Shape function values.
1411 const MEC& X, // Nodal coordinates.
1412 const MEC& U, // Nodal displacements.
1413 const double area, GradientContainer* gradient_container) {
1414 if (gradient_container == nullptr) { return; }
1415
1416 // Compute and store the integration point coordinate.
1417 const GradientContainer::InternalElementPosition ip = {xi[0], xi[1], xi[2], layer_id};
1418 gradient_container->set_current_position(ip, area);
1419 static const GradientContainer::FieldDescription coor_descr = {
1420 "COOR_IP",
1421 "x.y.z",
1422 "Physical coordinate of the point "
1423 "in the reference configuration",
1424 3,
1425 3,
1426 1,
1427 3,
1428 false,
1429 typeid(double)};
1430 if (gradient_container->compute_field_value(coor_descr.name)) {
1431 VD pcoor(3);
1432 pcoor = H * X;
1433 gradient_container->set_field_value(coor_descr, &pcoor[0]);
1434 }
1435
1436 // Compute and store the displacement at the integration
1437 // point.
1438 static const GradientContainer::FieldDescription disp_descr = {
1439 "DISP_IP", "dx.dy.dz", "Physical displacement at the point", 3, 3, 1, 3,
1440 false, typeid(double)};
1441 if (gradient_container->compute_field_value(disp_descr.name)) {
1442 VD pdisp(3);
1443 pdisp = H * U;
1444 gradient_container->set_field_value(disp_descr, &pdisp[0]);
1445 }
1446}
1447
1448//----------------------------------------------------------------------
1449// write_ip_coor_disp
1450//----------------------------------------------------------------------
1451template <int NNEL>
1452inline void write_ip_coor_disp_3d(
1453 const double* xi, // Element-local coor.
1454 const int layer_id, // ID of material layer.
1455 const double H[NNEL], // Shape function values.
1456 const double X[NNEL][3], // Nodal coordinates.
1457 const double U[NNEL][3], // Nodal displacements.
1458 const double area, GradientContainer* gradient_container) {
1459 if (gradient_container == nullptr) { return; }
1460
1461 // Compute and store the integration point coordinate.
1462 const GradientContainer::InternalElementPosition ip = {xi[0], xi[1], xi[2], layer_id};
1463 gradient_container->set_current_position(ip, area);
1464 static const GradientContainer::FieldDescription coor_descr = {
1465 "COOR_IP",
1466 "x.y.z",
1467 "Physical coordinate of the point "
1468 "in the reference configuration",
1469 3,
1470 3,
1471 1,
1472 3,
1473 false,
1474 typeid(double)};
1475 if (gradient_container->compute_field_value(coor_descr.name)) {
1476 double pcoor[3];
1477 ffmv::vc_eq_va_mul_mb<3, NNEL>(pcoor, H, X);
1478 gradient_container->set_field_value(coor_descr, &pcoor[0]);
1479 }
1480
1481 // Compute and store the displacement at the integration
1482 // point.
1483 static const GradientContainer::FieldDescription disp_descr = {
1484 "DISP_IP", "dx.dy.dz", "Physical displacement at the point", 3, 3, 1, 3,
1485 false, typeid(double)};
1486 if (gradient_container->compute_field_value(disp_descr.name)) {
1487 double pdisp[3];
1488 ffmv::vc_eq_va_mul_mb<3, NNEL>(pdisp, H, U);
1489 gradient_container->set_field_value(disp_descr, &pdisp[0]);
1490 }
1491}
1492
1493//----------------------------------------------------------------------
1494// write_ip_coor_disp
1495//----------------------------------------------------------------------
1496template <int NNEL>
1497inline void write_ip_coor_disp_2d(
1498 const double* xi, // Element-local coor.
1499 const int layer_id, // ID of material layer.
1500 const double H[NNEL], // Shape function values.
1501 const double X[NNEL][2], // Nodal coordinates.
1502 const double U[NNEL][2], // Nodal displacements.
1503 const double area, GradientContainer* gradient_container) {
1504 if (gradient_container == nullptr) { return; }
1505
1506 // Compute and store the integration point coordinate.
1507 const GradientContainer::InternalElementPosition ip = {xi[0], xi[1], xi[2], layer_id};
1508 gradient_container->set_current_position(ip, area);
1509 static const GradientContainer::FieldDescription coor_descr = {
1510 "COOR_IP",
1511 "x.y.z",
1512 "Physical coordinate of the point "
1513 "in the reference configuration",
1514 3,
1515 3,
1516 1,
1517 3,
1518 false,
1519 typeid(double)};
1520 if (gradient_container->compute_field_value(coor_descr.name)) {
1521 double pcoor[3] = {};
1522 ffmv::vc_eq_va_mul_mb<2, NNEL>(pcoor, H, X);
1523 gradient_container->set_field_value(coor_descr, &pcoor[0]);
1524 }
1525
1526 // Compute and store the displacement at the integration
1527 // point.
1528 static const GradientContainer::FieldDescription disp_descr = {
1529 "DISP_IP", "dx.dy.dz", "Physical displacement at the point", 3, 3, 1, 3,
1530 false, typeid(double)};
1531 if (gradient_container->compute_field_value(disp_descr.name)) {
1532 double pdisp[3] = {};
1533 ffmv::vc_eq_va_mul_mb<2, NNEL>(pdisp, H, U);
1534 gradient_container->set_field_value(disp_descr, &pdisp[0]);
1535 }
1536}
1537
1538//----------------------------------------------------------------------
1539// add_w1_to_second_2d
1540//
1541// Description:
1542//
1543// This routine adds 'area * B^T * C * B' to the second
1544// variation. The loops are arranged such that the number of
1545// floating-point operations is minimized.
1546//
1547// Assumptions:
1548//
1549// The internal storage formats and access patterns for
1550// upper- and lower-packed matrices is assumed. The matrix
1551// d_value_d_dof must be upper-packed, while the matrix C
1552// must be lower-packed.
1553//----------------------------------------------------------------------
1554template <int NDOF>
1555inline void add_w1_to_second_2d(
1556 const double BT[NDOF][3], const double C[6], const double volume,
1557 MP& d_value_d_dof // Stiffness matrix.
1558) {
1559 double* second_0_0;
1560 double* second_0_j;
1561 const double* b_i;
1562 const double* b_j;
1563 double t[3];
1564 double s;
1565
1566 assert(d_value_d_dof.is_upper());
1567 second_0_0 = &d_value_d_dof(0, 0);
1568
1569 for (int j = 0; j != NDOF; ++j) {
1570 b_j = BT[j];
1571 t[0] = b_j[0] * C[0] + b_j[1] * C[1] + b_j[2] * C[2];
1572 t[1] = b_j[0] * C[1] + b_j[1] * C[3] + b_j[2] * C[4];
1573 t[2] = b_j[0] * C[2] + b_j[1] * C[4] + b_j[2] * C[5];
1574
1575 t[0] *= volume;
1576 t[1] *= volume;
1577 t[2] *= volume;
1578
1579 second_0_j = second_0_0 + (j * (j + 1) / 2);
1580
1581 for (int i = 0; i <= j; ++i) {
1582 b_i = BT[i];
1583 s = b_i[0] * t[0] + b_i[1] * t[1] + b_i[2] * t[2];
1584 second_0_j[i] += s;
1585 }
1586 }
1587}
1588
1589template <int NDOF>
1590inline void add_w1_to_second_2d(
1591 const ME& B, // Strain-displ. matrix.
1592 const MLP& C, // Material matrix.
1593 const double area,
1594 MP& d_value_d_dof // Stiffness matrix.
1595) {
1596 double* second_0_0;
1597 double* second_0_j;
1598 const double* c;
1599 const double* b_i;
1600 const double* b_j;
1601 double t[3];
1602 double s;
1603
1604#ifdef TEST
1605 MP test(NDOF, true);
1606 test = area * transposed(B) * C * B;
1607 test += d_value_d_dof;
1608#endif /* TEST */
1609
1610 second_0_0 = &d_value_d_dof(0, 0);
1611 c = &C(0, 0);
1612
1613 for (int j = 0; j < NDOF; j++) {
1614 b_j = &B(0, j);
1615 t[0] = b_j[0] * c[0] + b_j[1] * c[1] + b_j[2] * c[2];
1616 t[1] = b_j[0] * c[1] + b_j[1] * c[3] + b_j[2] * c[4];
1617 t[2] = b_j[0] * c[2] + b_j[1] * c[4] + b_j[2] * c[5];
1618
1619 t[0] *= area;
1620 t[1] *= area;
1621 t[2] *= area;
1622
1623 second_0_j = second_0_0 + (j * (j + 1) / 2);
1624
1625 for (int i = 0; i <= j; i++) {
1626 b_i = &B(0, i);
1627 s = b_i[0] * t[0] + b_i[1] * t[1] + b_i[2] * t[2];
1628 second_0_j[i] += s;
1629 }
1630 }
1631
1632#ifdef TEST
1633 double sd = 0.0;
1634 for (int j = 0; j < NDOF; j++) {
1635 for (int i = 0; i < NDOF; i++) {
1636 double diff = test(i, j) - d_value_d_dof(i, j);
1637 sd += diff * diff;
1638 }
1639 }
1640 if (sd > 1e-10) {
1641 std::cerr << "Error in add_w1_to_second: sd=" << sd << std::endl;
1642 exit(1);
1643 }
1644#endif /* TEST */
1645}
1646
1647//----------------------------------------------------------------------
1648// add_w2_to_second_2d
1649//
1650// Description:
1651//
1652// This routine adds the "non-linear" stiffness contribution
1653// to the tangent stiffness (2nd variation) matrix.
1654//
1655// Assumptions:
1656//
1657// The internal storage formats and access patterns for
1658// upper-packed and rectangular matrices is assumed.
1659//
1660// Complexity:
1661//
1662// Be n the number of nodes of the element. The number of
1663// floating-point additions is
1664//
1665// 3n^2 + 9n
1666//
1667// The number of floating-point multiplications is:
1668//
1669// 2n^2 + 11n
1670//
1671// The number of integer operations (increments, additions,
1672// multiplications, divisions) is
1673//
1674// 4n^2 + 19n
1675//
1676//----------------------------------------------------------------------
1677template <int NDOF>
1678inline void add_w2_to_second_2d(
1679 const ME& DC, // Derivatives.
1680 const VD& S, // Stress tensor.
1681 const double area,
1682 MP& d_value_d_dof // Stiffness matrix.
1683) {
1684 const double* dc_i;
1685 const double* dc_j;
1686 double* second_0_0;
1687 double* second_0_j0;
1688 double* second_0_j1;
1689 double b[3];
1690 int i0;
1691 int i1;
1692 int j0;
1693 int j1;
1694
1695 const int NNEL = NDOF / 2;
1696
1697#ifdef TEST
1698 MP test(NDOF, true);
1699 ME W2(NDOF, NDOF);
1700 W2.set_zero();
1701 {
1702 ME BNL(4, NDOF);
1703 BNL.set_zero();
1704 for (int k = 0; k < NNEL; k++) {
1705 int k0 = 2 * k + 0;
1706 int k1 = 2 * k + 1;
1707
1708 BNL(0, k0) = DC(0, k);
1709 BNL(1, k0) = DC(1, k);
1710 BNL(3, k1) = DC(0, k);
1711 BNL(4, k1) = DC(1, k);
1712 }
1713
1714 ME SM(9, 9);
1715 SM.set_zero();
1716 for (int i = 0; i < 2; i++) {
1717 int j = 2 * i;
1718 SM(j + 0, j + 0) = S[0]; // S(0, 0);
1719 SM(j + 0, j + 1) = S[2]; // S(0, 1);
1720 SM(j + 1, j + 0) = S[2]; // S(1, 0);
1721 SM(j + 1, j + 1) = S[1]; // S(1, 1);
1722 }
1723 ME SMBNL(4, NDOF);
1724 SMBNL = SM * BNL;
1725 W2 = area * transposed(BNL) * SMBNL;
1726
1727 test.set_zero();
1728 for (int j = 0; j < NDOF; j++) {
1729 for (int i = 0; i <= j; i++) { test(i, j) = W2(i, j); }
1730 }
1731 test += d_value_d_dof;
1732 }
1733#endif /* TEST */
1734
1735 second_0_0 = &d_value_d_dof(0, 0);
1736 for (int j = 0; j < NNEL; j++) {
1737 dc_j = &DC(0, j);
1738
1739 j0 = j * 2 + 0;
1740 j1 = j * 2 + 1;
1741
1742 b[0] = S[0] * dc_j[0] + S[2] * dc_j[1];
1743 b[1] = S[2] * dc_j[0] + S[1] * dc_j[1];
1744
1745 second_0_j0 = second_0_0 + (j0 * (j0 + 1) / 2);
1746 second_0_j1 = second_0_0 + (j1 * (j1 + 1) / 2);
1747
1748 for (int i = 0; i <= j; i++) {
1749 dc_i = &DC(0, i);
1750
1751 double s = dc_i[0] * b[0] + dc_i[1] * b[1];
1752 s *= area;
1753
1754 i0 = i * 2 + 0;
1755 i1 = i * 2 + 1;
1756
1757 second_0_j0[i0] += s;
1758 second_0_j1[i1] += s;
1759 }
1760 }
1761
1762#ifdef TEST
1763 double sd = 0.0;
1764 for (int j = 0; j < NDOF; j++) {
1765 for (int i = 0; i < NDOF; i++) {
1766 double diff = test(i, j) - d_value_d_dof(i, j);
1767 sd += diff * diff;
1768 }
1769 }
1770 if (sd > 1e-10) {
1771 std::cerr << "Error in add_w2_to_second: sd=" << sd << std::endl;
1772 exit(1);
1773 }
1774#endif /* TEST */
1775}
1776
1777//----------------------------------------------------------------------
1778// add_to_mass_matrix_2d
1779//
1780// Add the contribution of the integration point to
1781// the mass matrix.
1782//
1783// Assumptions:
1784//
1785// - The mass matrix is computed for the reference configuration,
1786// thus it is assumed that the deformations are small.
1787//
1788// - The mass matrix must be in upper-packed form. The internal
1789// storage format and access pattern is assumed.
1790//
1791//----------------------------------------------------------------------
1792template <int NDOF>
1793void add_to_mass_matrix_2d(
1794 const VDC& H, // Shape functions.
1795 const double density, // Density at integr. point.
1796 const double area, // Integration area.
1797 MP& d_value_d_dofdotdot // Mass matrix.
1798) {
1799 const int NNEL = NDOF / 2;
1800
1801 ME HS(2, NDOF);
1802 for (int i = 0; i < 2; i++) {
1803 for (int j = 0; j < NNEL; j++) { HS(i, i + 2 * j) = H[j]; }
1804 }
1805
1806#ifdef TEST
1807 MP test(NDOF, true);
1808 for (int kl = 0; kl < NDOF; kl++) {
1809 for (int mn = 0; mn <= kl; mn++) { test(kl, mn) = HS[kl] * HS[mn] * area * density; }
1810 }
1811 test += d_value_d_dofdotdot;
1812#endif /* TEST */
1813
1814 const double dv = density * area;
1815
1816 double* mass_0_0 = &d_value_d_dofdotdot(0, 0);
1817 for (int kl = 0; kl < NDOF; kl++) {
1818 double* mass_0_kl = mass_0_0 + (kl * (kl + 1) / 2);
1819 for (int mn = 0; mn <= kl; mn++) { mass_0_kl[mn] += HS[kl] * HS[mn] * dv; }
1820 }
1821
1822#ifdef TEST
1823 double sd = 0.0;
1824 double avg = 0.0;
1825 for (int j = 0; j < NDOF; j++) {
1826 for (int i = 0; i < NDOF; i++) {
1827 avg += test(i, j) * test(i, j);
1828 double diff = test(i, j) - d_value_d_dofdotdot(i, j);
1829 sd += diff * diff;
1830 }
1831 }
1832 avg /= (NDOF * (NDOF + 1));
1833 // std::cerr << "avg=" << avg << " sd=" << sd << std::endl;
1834 if (sd > 1e-10 * avg) {
1835 std::cerr << "Error in add_to_mass_matrix: sd=" << sd << std::endl;
1836 exit(1);
1837 }
1838#endif /* TEST */
1839}
1840
1841//----------------------------------------------------------------------
1842// add_to_mass_matrix_2d
1843//
1844// Add the contribution of the integration point to
1845// the mass matrix.
1846//
1847// Assumptions:
1848//
1849// - The mass matrix is computed for the reference configuration,
1850// thus it is assumed that the deformations are small.
1851//
1852// - The mass matrix must be in upper-packed form. The internal
1853// storage format and access pattern is assumed.
1854//
1855//----------------------------------------------------------------------
1856template <int NNEL>
1857void add_to_mass_matrix_2d(
1858 const double H[NNEL], // Shape functions.
1859 const double density, // Density at integr. point.
1860 const double volume, // Integration volume.
1861 MP& d_value_d_dofdotdot // Mass matrix.
1862) {
1863 assert(d_value_d_dofdotdot.is_upper());
1864 const int NDOF = NNEL * 2;
1865
1866 const double dv = density * volume;
1867
1868 double* mass_0_0 = &d_value_d_dofdotdot(0, 0);
1869 for (int kl = 0; kl != NDOF; ++kl) {
1870 const int k = kl / 2;
1871 const int l = kl % 2;
1872 const double a = H[k] * dv;
1873 double* mass_0_kl = mass_0_0 + (kl * (kl + 1) / 2);
1874 for (int m = 0; m <= k; ++m) {
1875 int mn = 2 * m + l;
1876 mass_0_kl[mn] += a * H[m];
1877 }
1878 }
1879}
1880
1881//----------------------------------------------------------------------
1882// write_ip_coor_disp_2d
1883//----------------------------------------------------------------------
1884inline void write_ip_coor_disp_2d(
1885 const double* xi, // Element-local coor.
1886 const int layer_id, // ID of material layer.
1887 const VDC& H, // Shape function values.
1888 const MEC& X, // Nodal coordinates.
1889 const MEC& U, // Nodal displacements.
1890 const double area, GradientContainer* gradient_container) {
1891 if (gradient_container == nullptr) { return; }
1892
1893 // Compute and store the integration point coordinate.
1894 const GradientContainer::InternalElementPosition ip = {xi[0], xi[1], xi[2], layer_id};
1895 gradient_container->set_current_position(ip, area);
1896 static const GradientContainer::FieldDescription coor_descr = {
1897 "COOR_IP",
1898 "x.y.z",
1899 "Physical coordinate of the point "
1900 "in the reference configuration",
1901 3,
1902 3,
1903 1,
1904 3,
1905 false,
1906 typeid(double)};
1907 if (gradient_container->compute_field_value(coor_descr.name)) {
1908 VD pcoor(2);
1909 pcoor = H * X;
1910 VD pcoor_3d(3);
1911 pcoor_3d[0] = pcoor[0];
1912 pcoor_3d[1] = pcoor[1];
1913 pcoor_3d[2] = 0.0;
1914 gradient_container->set_field_value(coor_descr, &pcoor_3d[0]);
1915 }
1916
1917 // Compute and store the displacement at the integration
1918 // point.
1919 static const GradientContainer::FieldDescription disp_descr = {
1920 "DISP_IP", "dx.dy.dz", "Physical displacement at the point", 3, 3, 1, 3,
1921 false, typeid(double)};
1922 if (gradient_container->compute_field_value(disp_descr.name)) {
1923 VD pdisp(2);
1924 pdisp = H * U;
1925 VD pdisp_3d(3);
1926 pdisp_3d[0] = pdisp[0];
1927 pdisp_3d[1] = pdisp[1];
1928 pdisp_3d[2] = 0.0;
1929 gradient_container->set_field_value(disp_descr, &pdisp_3d[0]);
1930 }
1931}
1932
1933} // namespace b2000
1934
1935#endif /* __B2ELEMENT_CONTINUUM_UTIL_H__ */
b2000::details::operator<
bool operator<(const csda< T1 > &a, const csda< T2 > &b)
Comparison of two csda numbers is performed on the real part only.
Definition b2csda.H:262
b2000
Contains the base classes for implementing Finite Elements.
Definition b2boundary_condition.H:32