21#ifndef B2ELEMENT_KLT_COMPUTE_A3_H_
22#define B2ELEMENT_KLT_COMPUTE_A3_H_
24#include "b2element_klt_util.H"
25#include "utils/b2linear_algebra.H"
30void compute_A3(
const double A1[3][3],
const double A2[3][3],
double A3[3][3]);
32void compute_A3(
const Vec3d A1[3],
const Vec3d A2[3], Vec3d A3[3]);
38 ComputeA3d(
const size_t ndof_,
const bool first,
const bool second) {
39 resize(ndof_, first, second);
42 size_t get_ndof()
const {
return ndof; }
44 void resize(
const size_t ndof_,
const bool first,
const bool second) {
46 c.resize(ndof, first, second);
47 n.resize(ndof, first, second);
51 const XiDofVec3d<3>& a1,
const XiDofVec3d<3>& a2,
const bool first,
const bool second,
52 const bool a3_second, XiDofVec3d<3>& a3) {
54 compute_d2(a1, a2, first, second, a3);
56 compute_d1(a1, a2, first, second, a3);
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,
65 compute_d2_ts(a1, a2, w1, w2, first, second, a3);
67 compute_d1_ts(a1, a2, w1, w2, first, second, a3);
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) {
75 compute_an_d2(a1, a2, w1, w2, first, second, an);
77 compute_an_d1(a1, a2, w1, w2, first, second, an);
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);
87 const XiDofVec3d<1>& a1,
const XiDofVec3d<1>& a2,
const bool first,
const bool second,
91 const XiDofVec3d<1>& a1,
const XiDofVec3d<1>& a2,
const bool first,
const bool second,
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);
99 const b2linalg::Matrix<double, b2linalg::Mrectangle>& N,
const Vec3d pos[6],
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); }
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;
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]);
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]);
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;
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]));
145 void compute_d00_kl_consistent(
146 const b2linalg::Matrix<double, b2linalg::Mrectangle>& N,
const Vec3d pos[6],
148 const unsigned int nnode = (
unsigned int)N.size1();
149 const unsigned int ndof_ = (
unsigned int)N.size1() * 3;
151 assert(N.size2() >= 3);
152 if (nn.size() != ndof || !nn.second) {
153 nn.resize(ndof,
true,
true);
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); }
160 if (aa1.size() != ndof) { aa1.resize(ndof,
true,
false); }
161 if (aa2.size() != ndof) { aa2.resize(ndof,
true,
false); }
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];
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);
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;
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];
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]);
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]);
207 nn.d1[ii][2] = N1 * (-pos[d01][0]) - N2 * (-pos[d10][0]);
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]);
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]));
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;
240 nn.d2[m + j1 + 2][0] = +NN;
241 nn.d2[m + j2 + 1][0] = -NN;
243 nn.d2[m + j0 + 2][1] = -NN;
244 nn.d2[m + j2 + 0][1] = +NN;
246 nn.d2[m + j0 + 1][2] = +NN;
247 nn.d2[m + j1 + 0][2] = -NN;
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]
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;
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]); }
275 if (a3nnl.size() != ndof) { a3nnl.resize(ndof); }
276 for (
unsigned int i = 0; i != ndof_; ++i) { a3nnl[i] = l_inv2 * a3nn[i]; }
278 const Vec3d a3l = a3.d0 * l_inv;
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);
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;
301 double f = -dot(a3.d0, nn.d2[m])
302 + l_inv * (3 * a3nn[i] * a3nn[j] - dot(nn.d1[i], nn.d1[j]));
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];
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;
331 norm(a3.d2[i * ndof + j] - aa3[d00].d2[i * ndof + j]) * ll_inv;
333 std::cerr <<
"i=" << i <<
" j=" << j <<
" a3=" << a3.d2[i * ndof + j]
334 <<
" aa3=" << aa3[d00].d2[i * ndof + j] <<
" diff=" << diff
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); }
354 const Vec3d& a1 = pos[d10];
355 const Vec3d& a2 = pos[d01];
358 const double a11 = dot(a1, a1);
359 const double a12 = dot(a1, a2);
360 const double a22 = dot(a2, a2);
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;
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;
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;
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]));
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;
393 assert(N.size2() >= 3);
394 if (ann.size() != ndof || !ann.second) {
395 ann.resize(ndof,
true,
true);
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); }
402 const Vec3d& a1 = pos[d10];
403 const Vec3d& a2 = pos[d01];
406 const double a11 = dot(a1, a1);
407 const double a12 = dot(a1, a2);
408 const double a22 = dot(a2, a2);
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;
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;
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;
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]));
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];
451 Vec3d& v = ann.d2[m];
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;
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]); }
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]; }
467 const Vec3d anl = an.d0 * l_inv;
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);
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]));
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];
495 const XiDofVec3d<3>& a1,
const XiDofVec3d<3>& a2,
const bool first,
const bool second,
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);
503 const XiDofVec3d<3>& a1,
const XiDofVec3d<3>& a2,
const bool first,
const bool second,
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);
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);
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);
524 std::vector<double> an_nn;
525 std::vector<double> an_nnl;
529 std::vector<double> a3nn;
530 std::vector<double> a3nnl;
540 const XiDofVec3d<3>& a1,
const XiDofVec3d<3>& a2,
const bool first,
const bool second,
541 XiDofVec3d<3>& n, ScalarXiDof<3>& c);
544void normalize(
const XiDofVec3d<3>& a, XiDofVec3d<3>& n);
546void normalize(
const XiDofVec3d<1>& a, XiDofVec3d<1>& n);
csda< T > norm(const csda< T > &a)
Definition b2csda.H:343