b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2element_klt_compute_a3.H
1//------------------------------------------------------------------------
2// b2element_klt_compute_a3.H --
3//
4// Compute the surface normal and its derivatives.
5//
6// written by Thomas Ludwig
7// Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
8//
9// (c) 2015,2016,2017 SMR Engineering & Development SA
10// 2502 Bienne, Switzerland
11//
12// (c) 2023 Deutsches Zentrum für Luft- und Raumfahrt (DLR) e.V.
13// Linder Höhe, 51147 Köln
14//
15// All Rights Reserved. Proprietary source code. The contents of
16// this file may not be disclosed to third parties, copied or
17// duplicated in any form, in whole or in part, without the prior
18// written permission of SMR.
19//------------------------------------------------------------------------
20
21#ifndef B2ELEMENT_KLT_COMPUTE_A3_H_
22#define B2ELEMENT_KLT_COMPUTE_A3_H_
23
24#include "b2element_klt_util.H"
25#include "utils/b2linear_algebra.H"
27
28namespace b2000::klt {
29
30void compute_A3(const double A1[3][3], const double A2[3][3], double A3[3][3]);
31
32void compute_A3(const Vec3d A1[3], const Vec3d A2[3], Vec3d A3[3]);
33
34class ComputeA3d {
35public:
36 ComputeA3d() {}
37
38 ComputeA3d(const size_t ndof_, const bool first, const bool second) {
39 resize(ndof_, first, second);
40 }
41
42 size_t get_ndof() const { return ndof; }
43
44 void resize(const size_t ndof_, const bool first, const bool second) {
45 ndof = ndof_;
46 c.resize(ndof, first, second);
47 n.resize(ndof, first, second);
48 }
49
50 void compute(
51 const XiDofVec3d<3>& a1, const XiDofVec3d<3>& a2, const bool first, const bool second,
52 const bool a3_second, XiDofVec3d<3>& a3) {
53 if (a3_second) {
54 compute_d2(a1, a2, first, second, a3);
55 } else {
56 compute_d1(a1, a2, first, second, a3);
57 }
58 }
59
60 void compute_ts(
61 const XiDofVec3d<3>& a1, const XiDofVec3d<3>& a2, const ScalarXiDof<3>& w1,
62 const ScalarXiDof<3>& w2, const bool first, const bool second, const bool a3_second,
63 XiDofVec3d<3>& a3) {
64 if (a3_second) {
65 compute_d2_ts(a1, a2, w1, w2, first, second, a3);
66 } else {
67 compute_d1_ts(a1, a2, w1, w2, first, second, a3);
68 }
69 }
70
71 void compute_an(
72 const XiDofVec3d<3>& a1, const XiDofVec3d<3>& a2, const double w1, const double w2,
73 const bool first, const bool second, const bool a3_second, XiDofVec3d<3>& an) {
74 if (a3_second) {
75 compute_an_d2(a1, a2, w1, w2, first, second, an);
76 } else {
77 compute_an_d1(a1, a2, w1, w2, first, second, an);
78 }
79 }
80
81 void compute_an_d00(
82 const XiDofVec3d<1>& a1, const XiDofVec3d<1>& a2, const double w1, const double w2,
83 const bool first, const bool second, XiDofVec3d<1>& an);
84
85public:
86 void compute_d00(
87 const XiDofVec3d<1>& a1, const XiDofVec3d<1>& a2, const bool first, const bool second,
88 XiDofVec3d<1>& a3);
89
90 void compute_d00_2nd(
91 const XiDofVec3d<1>& a1, const XiDofVec3d<1>& a2, const bool first, const bool second,
92 XiDofVec3d<1>& a3);
93
94 void compute_d00_ts(
95 const XiDofVec3d<1>& a1, const XiDofVec3d<1>& a2, const ScalarXiDof<1>& w1,
96 const ScalarXiDof<1>& w2, const bool first, const bool second, XiDofVec3d<1>& a3);
97
98 void compute_d00_kl(
99 const b2linalg::Matrix<double, b2linalg::Mrectangle>& N, const Vec3d pos[6],
100 DofVec3d& a3) {
101 const unsigned int nnode = (unsigned int)N.size1();
102 const unsigned int ndof = (unsigned int)N.size1() * 3;
103 assert(N.size2() >= 3);
104 if (n.size() != ndof) { n.resize(ndof, true, false); }
105 if (a3.size() != ndof) { a3.resize(ndof, true, false); }
106
107 n[d00].d0 = outer_product(pos[d10], pos[d01]);
108 const double l = norm(n[d00].d0);
109 const double l_inv = 1. / l;
110 a3.d0 = n[d00].d0 * l_inv;
111
112 for (unsigned int i = 0; i != nnode; ++i) {
113 const unsigned int j = 0;
114 const unsigned int ii = i * 3 + j;
115 const double N1 = N[d10][i];
116 const double N2 = N[d01][i];
117 n[d00].d1[ii][0] = 0;
118 n[d00].d1[ii][1] = N1 * (-pos[d01][2]) - N2 * (-pos[d10][2]);
119 n[d00].d1[ii][2] = N1 * (+pos[d01][1]) - N2 * (+pos[d10][1]);
120 }
121 for (unsigned int i = 0; i != nnode; ++i) {
122 const unsigned int j = 1;
123 const unsigned int ii = i * 3 + j;
124 const double N1 = N[d10][i];
125 const double N2 = N[d01][i];
126 n[d00].d1[ii][0] = N1 * (+pos[d01][2]) - N2 * (+pos[d10][2]);
127 n[d00].d1[ii][1] = 0;
128 n[d00].d1[ii][2] = N1 * (-pos[d01][0]) - N2 * (-pos[d10][0]);
129 }
130 for (unsigned int i = 0; i != nnode; ++i) {
131 const unsigned int j = 2;
132 const unsigned int ii = i * 3 + j;
133 const double N1 = N[d10][i];
134 const double N2 = N[d01][i];
135 n[d00].d1[ii][0] = N1 * (-pos[d01][1]) - N2 * (-pos[d10][1]);
136 n[d00].d1[ii][1] = N1 * (+pos[d01][0]) - N2 * (+pos[d10][0]);
137 n[d00].d1[ii][2] = 0;
138 }
139
140 for (unsigned int i = 0; i != ndof; ++i) {
141 a3.d1[i] = l_inv * (n[d00].d1[i] - a3.d0 * dot(a3.d0, n[d00].d1[i]));
142 }
143 }
144
145 void compute_d00_kl_consistent(
146 const b2linalg::Matrix<double, b2linalg::Mrectangle>& N, const Vec3d pos[6],
147 DofVec3d& a3) {
148 const unsigned int nnode = (unsigned int)N.size1();
149 const unsigned int ndof_ = (unsigned int)N.size1() * 3;
150 ndof = ndof_;
151 assert(N.size2() >= 3);
152 if (nn.size() != ndof || !nn.second) {
153 nn.resize(ndof, true, true);
154 nn.set_zero();
155 }
156 if (n.size() != ndof || !n[d00].second) { n.resize(ndof, true, true); }
157 if (a3.size() != ndof || !a3.second) { a3.resize(ndof, true, true); }
158
159 if (0) {
160 if (aa1.size() != ndof) { aa1.resize(ndof, true, false); }
161 if (aa2.size() != ndof) { aa2.resize(ndof, true, false); }
162 aa1.set_zero();
163 aa2.set_zero();
164 aa1[d00].d0 = pos[d10];
165 aa2[d00].d0 = pos[d01];
166 for (unsigned int i = 0; i != nnode; ++i) {
167 const unsigned int ii = 3 * i;
168 for (int dir = 0; dir != 3; ++dir) {
169 aa1[d00].d1[ii + dir][dir] = N[d10][i];
170 aa2[d00].d1[ii + dir][dir] = N[d01][i];
171 }
172 }
173
174 if (c.size() != ndof) { c.resize(ndof, true, true); }
175 if (aa3.size() != ndof) { aa3.resize(ndof, true, true); }
176 compute_d00_2nd(aa1, aa2, true, true, aa3);
177
178 // a3 = aa3[d00];
179 // return;
180 }
181
182 nn.d0 = outer_product(pos[d10], pos[d01]);
183 const double l = norm(nn.d0);
184 const double l_inv = 1. / l;
185 const double l_inv2 = l_inv * l_inv;
186 a3.d0 = nn.d0 * l_inv;
187
188 //
189 // First derivatives.
190 //
191 for (unsigned int i = 0; i != nnode; ++i) {
192 const unsigned int j = 0;
193 const unsigned int ii = i * 3 + j;
194 const double N1 = N[d10][i];
195 const double N2 = N[d01][i];
196 nn.d1[ii][0] = 0;
197 nn.d1[ii][1] = N1 * (-pos[d01][2]) - N2 * (-pos[d10][2]);
198 nn.d1[ii][2] = N1 * (+pos[d01][1]) - N2 * (+pos[d10][1]);
199 }
200 for (unsigned int i = 0; i != nnode; ++i) {
201 const unsigned int j = 1;
202 const unsigned int ii = i * 3 + j;
203 const double N1 = N[d10][i];
204 const double N2 = N[d01][i];
205 nn.d1[ii][0] = N1 * (+pos[d01][2]) - N2 * (+pos[d10][2]);
206 nn.d1[ii][1] = 0;
207 nn.d1[ii][2] = N1 * (-pos[d01][0]) - N2 * (-pos[d10][0]);
208 }
209 for (unsigned int i = 0; i != nnode; ++i) {
210 const unsigned int j = 2;
211 const unsigned int ii = i * 3 + j;
212 const double N1 = N[d10][i];
213 const double N2 = N[d01][i];
214 nn.d1[ii][0] = N1 * (-pos[d01][1]) - N2 * (-pos[d10][1]);
215 nn.d1[ii][1] = N1 * (+pos[d01][0]) - N2 * (+pos[d10][0]);
216 nn.d1[ii][2] = 0;
217 }
218
219 for (unsigned int i = 0; i != ndof_; ++i) {
220 a3.d1[i] = l_inv * (nn.d1[i] - a3.d0 * dot(a3.d0, nn.d1[i]));
221 }
222
223 //
224 // Second derivatives.
225 //
226 {
227 const unsigned int j0 = 0;
228 const unsigned int j1 = ndof_;
229 const unsigned int j2 = 2 * ndof_;
230 for (unsigned int i = 0; i != nnode; ++i) {
231 const unsigned int ii = i * 3 * ndof_;
232 const double N1i = N[d10][i];
233 const double N2i = N[d01][i];
234 for (unsigned int k = 0; k <= i; ++k) {
235 const double N1k = N[d10][k];
236 const double N2k = N[d01][k];
237 const unsigned int m = ii + 3 * k;
238 const double NN = +N1i * N2k - N2i * N1k;
239
240 nn.d2[m + j1 + 2][0] = +NN;
241 nn.d2[m + j2 + 1][0] = -NN;
242
243 nn.d2[m + j0 + 2][1] = -NN;
244 nn.d2[m + j2 + 0][1] = +NN;
245
246 nn.d2[m + j0 + 1][2] = +NN;
247 nn.d2[m + j1 + 0][2] = -NN;
248 }
249 }
250 }
251
252 // Check.
253 if (0) {
254 bool ok = true;
255 for (unsigned int i = 0; i != ndof; ++i) {
256 if (nn.d1[i] != n[d00].d1[i]) {
257 std::cerr << "d1 i=" << i << " nn=" << nn.d1[i] << " n=" << n[d00].d1[i]
258 << std::endl;
259 ok = false;
260 }
261 for (unsigned int j = 0; j <= i; ++j) {
262 if (nn.d2[i * ndof + j] != n[d00].d2[i * ndof + j]) {
263 std::cerr << "i=" << i << " j=" << j << " nn=" << nn.d2[i * ndof + j]
264 << " n=" << n[d00].d2[i * ndof + j] << std::endl;
265 ok = false;
266 }
267 }
268 }
269 assert(ok);
270 }
271
272 if (a3nn.size() != ndof) { a3nn.resize(ndof); }
273 for (unsigned int i = 0; i != ndof_; ++i) { a3nn[i] = dot(a3.d0, nn.d1[i]); }
274
275 if (a3nnl.size() != ndof) { a3nnl.resize(ndof); }
276 for (unsigned int i = 0; i != ndof_; ++i) { a3nnl[i] = l_inv2 * a3nn[i]; }
277
278 const Vec3d a3l = a3.d0 * l_inv;
279
280 {
281 assert(a3.d1.size() == ndof);
282 assert(a3.d2.size1() == ndof && a3.d2.size2() == ndof);
283 assert(nn.d1.size() == ndof);
284 assert(nn.d2.size1() == ndof && nn.d2.size2() == ndof);
285 unsigned int ii = 0;
286 for (unsigned int i = 0; i != ndof_; ++i) {
287 for (unsigned int j = 0; j <= i; ++j) {
288 const unsigned int m = ii + j;
289 // assert(m < ndof * ndof);
290 // a3.d2[m] = l_inv * (nn.d2[m] - a3.d0 * dot(a3.d0, nn.d2[m])) + l_inv * l_inv
291 // * (a3.d0 * (3 * a3nn[i] * a3nn[j]) - a3.d0 * dot(nn.d1[i], nn.d1[j]) -
292 // nn.d1[i] * a3nn[j] - nn.d1[j] * a3nn[i]);
293 // a3.d2[m] = nn.d2[m] * l_inv +
294 // a3l * (
295 // -dot(a3.d0, nn.d2[m]) +
296 // l_inv * (3 * a3nn[i] * a3nn[j] -
297 // dot(nn.d1[i], nn.d1[j]))
298 // ) -
299 // nn.d1[i] * a3nnl[j] -
300 // nn.d1[j] * a3nnl[i];
301 double f = -dot(a3.d0, nn.d2[m])
302 + l_inv * (3 * a3nn[i] * a3nn[j] - dot(nn.d1[i], nn.d1[j]));
303 Vec3d& v = a3.d2[m];
304 // for (int dir = 0; dir != 3; ++dir) {
305 // v[dir] = nn.d2[m][dir] * l_inv + a3l[dir] * f -
306 // nn.d1[i][dir] * a3nnl[j] -
307 // nn.d1[j][dir] * a3nnl[i];
308 // }
309 v[0] = nn.d2[m][0] * l_inv + a3l[0] * f - nn.d1[i][0] * a3nnl[j]
310 - nn.d1[j][0] * a3nnl[i];
311 v[1] = nn.d2[m][1] * l_inv + a3l[1] * f - nn.d1[i][1] * a3nnl[j]
312 - nn.d1[j][1] * a3nnl[i];
313 v[2] = nn.d2[m][2] * l_inv + a3l[2] * f - nn.d1[i][2] * a3nnl[j]
314 - nn.d1[j][2] * a3nnl[i];
315 }
316 ii += ndof_;
317 }
318 }
319
320 // Check.
321 if (0) {
322 bool ok = true;
323 for (unsigned int i = 0; i != ndof; ++i) {
324 for (unsigned int j = 0; j <= i; ++j) {
325 const double l1 = norm(a3.d2[i * ndof + j]);
326 const double l2 = norm(aa3[d00].d2[i * ndof + j]);
327 const double ll = std::max(l1, l2);
328 if (ll < 1e-10) { continue; }
329 const double ll_inv = 1. / ll;
330 const double diff =
331 norm(a3.d2[i * ndof + j] - aa3[d00].d2[i * ndof + j]) * ll_inv;
332 if (diff > 1e-3) {
333 std::cerr << "i=" << i << " j=" << j << " a3=" << a3.d2[i * ndof + j]
334 << " aa3=" << aa3[d00].d2[i * ndof + j] << " diff=" << diff
335 << std::endl;
336 ok = false;
337 }
338 }
339 }
340 assert(ok);
341 }
342 }
343
344 void compute_an_d00_kl(
345 const b2linalg::Matrix<double, b2linalg::Mrectangle>& N, const Vec3d pos[6],
346 const double w1, const double w2, DofVec3d& an) {
347 const unsigned int nnode = (unsigned int)N.size1();
348 const unsigned int ndof = (unsigned int)N.size1() * 3;
349 assert(N.size2() >= 3);
350 if (ann.size() != ndof) { ann.resize(ndof, true, false); }
351 if (an.size() != ndof) { an.resize(ndof, true, false); }
352
353 // Aliases.
354 const Vec3d& a1 = pos[d10];
355 const Vec3d& a2 = pos[d01];
356
357 // Components of the metric tensor.
358 const double a11 = dot(a1, a1);
359 const double a12 = dot(a1, a2);
360 const double a22 = dot(a2, a2);
361
362 // No derivatives.
363 ann.d0 = -(a1 * a12 - a2 * a11) * w1 + (a2 * a12 - a1 * a22) * w2;
364 const double l = norm(ann.d0);
365 const double l_inv = 1. / l;
366 an.d0 = ann.d0 * l_inv;
367
368 // First derivatives of non-normalised vector.
369 for (unsigned int i = 0; i != nnode; ++i) {
370 const double N1 = N[d10][i];
371 const double N2 = N[d01][i];
372 for (int j = 0; j != 3; ++j) {
373 const unsigned int ii = i * 3 + j;
374 ann.d1[ii] =
375 -(a1 * (a1[j] * N2 + a2[j] * N1) - a2 * (a1[j] * N1 + a1[j] * N1)) * w1
376 + (a2 * (a1[j] * N2 + a2[j] * N1) - a1 * (a2[j] * N2 + a2[j] * N2)) * w2;
377 ann.d1[ii][j] += -(N1 * a12 - N2 * a11) * w1 + (N2 * a12 - N1 * a22) * w2;
378 }
379 }
380
381 // First derivatives of normalised vector.
382 for (unsigned int i = 0; i != ndof; ++i) {
383 an.d1[i] = l_inv * (ann.d1[i] - an.d0 * dot(an.d0, ann.d1[i]));
384 }
385 }
386
387 void compute_an_d00_kl_consistent(
388 const b2linalg::Matrix<double, b2linalg::Mrectangle>& N, const Vec3d pos[6],
389 const double w1, const double w2, DofVec3d& an) {
390 const unsigned int nnode = (unsigned int)N.size1();
391 const unsigned int ndof_ = (unsigned int)N.size1() * 3;
392 ndof = ndof_;
393 assert(N.size2() >= 3);
394 if (ann.size() != ndof || !ann.second) {
395 ann.resize(ndof, true, true);
396 ann.set_zero();
397 }
398 if (n.size() != ndof || !n[d00].second) { n.resize(ndof, true, true); }
399 if (an.size() != ndof || !an.second) { an.resize(ndof, true, true); }
400
401 // Aliases.
402 const Vec3d& a1 = pos[d10];
403 const Vec3d& a2 = pos[d01];
404
405 // Components of the metric tensor.
406 const double a11 = dot(a1, a1);
407 const double a12 = dot(a1, a2);
408 const double a22 = dot(a2, a2);
409
410 // No derivatives.
411 ann.d0 = -(a1 * a12 - a2 * a11) * w1 + (a2 * a12 - a1 * a22) * w2;
412 const double l = norm(ann.d0);
413 const double l_inv = 1. / l;
414 const double l_inv2 = l_inv * l_inv;
415 an.d0 = ann.d0 * l_inv;
416
417 // First derivatives of non-normalised vector.
418 for (unsigned int i = 0; i != nnode; ++i) {
419 const double N1 = N[d10][i];
420 const double N2 = N[d01][i];
421 for (int j = 0; j != 3; ++j) {
422 const unsigned int ii = i * 3 + j;
423 ann.d1[ii] =
424 -(a1 * (a1[j] * N2 + a2[j] * N1) - a2 * (a1[j] * N1 + a1[j] * N1)) * w1
425 + (a2 * (a1[j] * N2 + a2[j] * N1) - a1 * (a2[j] * N2 + a2[j] * N2)) * w2;
426 ann.d1[ii][j] += -(N1 * a12 - N2 * a11) * w1 + (N2 * a12 - N1 * a22) * w2;
427 }
428 }
429
430 // First derivatives of normalised vector.
431 for (unsigned int i = 0; i != ndof; ++i) {
432 an.d1[i] = l_inv * (ann.d1[i] - an.d0 * dot(an.d0, ann.d1[i]));
433 }
434
435 //
436 // Second derivatives.
437 //
438 {
439 for (unsigned int ij = 0; ij != ndof_; ++ij) {
440 const unsigned int i = ij / 3;
441 const unsigned int j = ij % 3;
442 const double N1i = N[d10][i];
443 const double N2i = N[d01][i];
444 for (unsigned int kl = 0; kl <= ij; ++kl) {
445 const unsigned int k = kl / 3;
446 const unsigned int l = kl % 3;
447 const unsigned m = ij * ndof_ + kl;
448 const double N1k = N[d10][k];
449 const double N2k = N[d01][k];
450
451 Vec3d& v = ann.d2[m];
452 v.set_zero();
453 v[j] += -(N1i * (N1k * a2[l] + N2k * a1[l]) - 2 * N2i * N1k * a1[l]) * w1
454 + (N2i * (N1k * a2[l] + N2k * a1[l]) - 2 * N1i * N2k * a2[l]) * w2;
455 v[l] += -(N1i * (N1k * a2[j] + N2k * a1[j]) - 2 * N2i * N1k * a1[j]) * w1
456 + (N2i * (N1k * a2[j] + N2k * a1[j]) - 2 * N1i * N2k * a2[j]) * w2;
457 }
458 }
459 }
460
461 if (an_nn.size() != ndof) { an_nn.resize(ndof); }
462 for (unsigned int i = 0; i != ndof_; ++i) { an_nn[i] = dot(an.d0, nn.d1[i]); }
463
464 if (an_nnl.size() != ndof) { an_nnl.resize(ndof); }
465 for (unsigned int i = 0; i != ndof_; ++i) { an_nnl[i] = l_inv2 * an_nn[i]; }
466
467 const Vec3d anl = an.d0 * l_inv;
468
469 {
470 assert(an.d1.size() == ndof);
471 assert(an.d2.size1() == ndof && an.d2.size2() == ndof);
472 assert(nn.d1.size() == ndof);
473 assert(nn.d2.size1() == ndof && nn.d2.size2() == ndof);
474 unsigned int ii = 0;
475 for (unsigned int i = 0; i != ndof_; ++i) {
476 for (unsigned int j = 0; j <= i; ++j) {
477 const unsigned int m = ii + j;
478 double f = -dot(an.d0, nn.d2[m])
479 + l_inv * (3 * an_nn[i] * an_nn[j] - dot(nn.d1[i], nn.d1[j]));
480 Vec3d& v = an.d2[m];
481 v[0] = nn.d2[m][0] * l_inv + anl[0] * f - nn.d1[i][0] * an_nnl[j]
482 - nn.d1[j][0] * an_nnl[i];
483 v[1] = nn.d2[m][1] * l_inv + anl[1] * f - nn.d1[i][1] * an_nnl[j]
484 - nn.d1[j][1] * an_nnl[i];
485 v[2] = nn.d2[m][2] * l_inv + anl[2] * f - nn.d1[i][2] * an_nnl[j]
486 - nn.d1[j][2] * an_nnl[i];
487 }
488 ii += ndof_;
489 }
490 }
491 }
492
493private:
494 void compute_d1(
495 const XiDofVec3d<3>& a1, const XiDofVec3d<3>& a2, const bool first, const bool second,
496 XiDofVec3d<3>& a3);
497
498 void compute_d1_ts(
499 const XiDofVec3d<3>& a1, const XiDofVec3d<3>& a2, const ScalarXiDof<3>& w1,
500 const ScalarXiDof<3>& w2, const bool first, const bool second, XiDofVec3d<3>& a3);
501
502 void compute_d2(
503 const XiDofVec3d<3>& a1, const XiDofVec3d<3>& a2, const bool first, const bool second,
504 XiDofVec3d<3>& a3);
505
506 void compute_d2_ts(
507 const XiDofVec3d<3>& a1, const XiDofVec3d<3>& a2, const ScalarXiDof<3>& w1,
508 const ScalarXiDof<3>& w2, const bool first, const bool second, XiDofVec3d<3>& a3);
509
510 void compute_an_d1(
511 const XiDofVec3d<3>& a1, const XiDofVec3d<3>& a2, const double w1, const double w2,
512 const bool first, const bool second, XiDofVec3d<3>& an);
513
514 void compute_an_d2(
515 const XiDofVec3d<3>& a1, const XiDofVec3d<3>& a2, const double w1, const double w2,
516 const bool first, const bool second, XiDofVec3d<3>& an);
517
518 size_t ndof;
519 ScalarXiDof<3> c;
520 XiDofVec3d<3> n;
521
522 // For the unit side-normal vector.
523 DofVec3d ann;
524 std::vector<double> an_nn;
525 std::vector<double> an_nnl;
526
527 // For now.
528 DofVec3d nn;
529 std::vector<double> a3nn;
530 std::vector<double> a3nnl;
531
532 XiDofVec3d<1> aa1;
533 XiDofVec3d<1> aa2;
534 XiDofVec3d<1> aa3;
535};
536
537class ComputeA3test {
538public:
539 void compute(
540 const XiDofVec3d<3>& a1, const XiDofVec3d<3>& a2, const bool first, const bool second,
541 XiDofVec3d<3>& n, ScalarXiDof<3>& c);
542};
543
544void normalize(const XiDofVec3d<3>& a, XiDofVec3d<3>& n);
545
546void normalize(const XiDofVec3d<1>& a, XiDofVec3d<1>& n);
547
548} // namespace b2000::klt
549
550#endif // B2ELEMENT_KLT_COMPUTE_A3_H_
csda< T > norm(const csda< T > &a)
Definition b2csda.H:343