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