/* See {ElasticityEnergy.h} */ #include #include #include #include #include #include #include #include CONST Pi == Math.Pi; A == LR3x3.T { (r3_t){0.0, 1.0, 1.0}, (r3_t){1.0, 0.0, 1.0}, (r3_t){1.0, 1.0, 0.0} } /* Matrix that maps the reference configuration to the repose (or relaxed) configuration (i.e. has zero elastic energy). The volume is 0.333333 . The reference configuration is a tetrahedron with nodes: r1==(0,0,0), r2==(1,0,0), r3==(0,1,0) and r4==(0,0,1). The repose configuration is a regular tetrahedron with nodes: q1==(0,0,0), q2==(0,1,1), q3==(1,1,0) and q4==(1,0,1). The Elasticity energy forces each tetrahedron has volume equal to 0.3333. */ TYPE bool_vec_t == ARRAY OF bool_t; double_vec_t == double_vec_t; LR4x3 == ARRAY [0..3] OF r3_t; LR3x4 == ARRAY [0..2] OF r4_t; REVEAL T == Public BRANDED OBJECT double K; /* The energy normalization factor */ ElemTableRec_t *top; /* The map elements. */ bool_vec_t vVar; /* TRUE if node is variable */ polyRelevant: REF bool_vec_t; /* TRUE if cell is relevant */ ep: REF double_vec_t; /* elasticity energy of each tetrahedron */ eDep: REF double_vec_t; /* (Work) derivative of energy rel. ep */ double volume; OVERRIDES init = Init; defTop = DefTop; defVar = DefVar; eval = Eval; name = Name; } SELF_Init(SELF_T *erg) { return erg; } /* END Init */ void SELF_DefTop(SELF_T *erg, ElemTableRec_t *top) { erg->K = 1.0; erg->top = top; erg->volume = (4.0/3.0)* FLOAT(Pi,double)/FLOAT(top->cell.nel,double); erg->vVar = bool_vec_new(top->node.nel); erg->polyRelevant = bool_vec_new(top->cell.nel); erg->ep = double_vec_new(top->cell.nel); erg->eDep = double_vec_new(top->cell.nel); /* Allocate determinant tables: */ /* Just in case the client forgets to call "defVar": */ for (i = 0; i < top->node.nel; i++){ erg->vVar.el[i] = FALSE; } for (i = 0; i < top->cell.nel; i++){ erg->polyRelevant[i] = FALSE; } } /* END DefTop */ LR4x3 Transpose_3x4(*m: LR3x4) /* Return the transpose of a matrix 3x4. */ t : LR4x3; { t[0] = (r3_t){m[0,0], m[1,0], m[2,0]}; t[1] = (r3_t){m[0,1], m[1,1], m[2,1]}; t[2] = (r3_t){m[0,2], m[1,2], m[2,2]}; t[3] = (r3_t){m[0,3], m[1,3], m[2,3]}; return t; } /* END Transpose_3x4 */ LR3x4 Transpose_4x3(*m: LR4x3) /* return the transpose of a matrix 4x3. */ t : LR3x4; { t[0] = (r4_t){m[0,0], m[1,0], m[2,0], m[3,0]}; t[1] = (r4_t){m[0,1], m[1,1], m[2,1], m[3,1]}; t[2] = (r4_t){m[0,2], m[1,2], m[2,2], m[3,2]}; return t; } /* END Transpose_4x3 */ LR4x3 Mul_4x3_3x3(*a: LR4x3; *b: LR3x3.T) /* Return the product of the matrix "a" 4x3 and matrix "b" 3x3. */ c : LR4x3; { double a00 = a[0,0], a01 = a[0,1], a02 = a[0,2]; double a10 = a[1,0], a11 = a[1,1], a12 = a[1,2]; double a20 = a[2,0], a21 = a[2,1], a22 = a[2,2]; double a30 = a[3,0], a31 = a[3,1], a32 = a[3,2]; double b00 = b[0,0], b01 = b[0,1], b02 = b[0,2]; double b10 = b[1,0], b11 = b[1,1], b12 = b[1,2]; double b20 = b[2,0], b21 = b[2,1], b22 = b[2,2]; double c00 = a00 * b00 + a01 * b10 + a02 * b20; double c01 = a00 * b01 + a01 * b11 + a02 * b21; double c02 = a00 * b02 + a01 * b12 + a02 * b22; double c10 = a10 * b00 + a11 * b10 + a12 * b20; double c11 = a10 * b01 + a11 * b11 + a12 * b21; double c12 = a10 * b02 + a11 * b12 + a12 * b22; double c20 = a20 * b00 + a21 * b10 + a22 * b20; double c21 = a20 * b01 + a21 * b11 + a22 * b21; double c22 = a20 * b02 + a21 * b12 + a22 * b22; double c30 = a30 * b00 + a31 * b10 + a32 * b20; double c31 = a30 * b01 + a31 * b11 + a32 * b21; double c32 = a30 * b02 + a31 * b12 + a32 * b22; r3_t c0 = (r3_t){c00,c01,c02}; r3_t c1 = (r3_t){c10,c11,c12}; r3_t c2 = (r3_t){c20,c21,c22}; r3_t c3 = (r3_t){c30,c31,c32}; { c[0] = c0; c[1] = c1; c[2] = c2; c[3] = c3; return c; } } /* END Mul_4x3_3x3 */ LR3x3.T Mul_3x4_4x3(*a: LR3x4; *b: LR4x3) /* Return the product of the matrix "a" 3x4 and matrix "b" 4x3. */ { double a00 = a[0,0], a01 = a[0,1], a02 = a[0,2], a03 = a[0,3]; double a10 = a[1,0], a11 = a[1,1], a12 = a[1,2], a13 = a[1,3]; double a20 = a[2,0], a21 = a[2,1], a22 = a[2,2], a23 = a[2,3]; double b00 = b[0,0], b01 = b[0,1], b02 = b[0,2]; double b10 = b[1,0], b11 = b[1,1], b12 = b[1,2]; double b20 = b[2,0], b21 = b[2,1], b22 = b[2,2]; double b30 = b[3,0], b31 = b[3,1], b32 = b[3,2]; double c00 = a00 * b00 + a01 * b10 + a02 * b20 + a03 * b30; double c01 = a00 * b01 + a01 * b11 + a02 * b21 + a03 * b31; double c02 = a00 * b02 + a01 * b12 + a02 * b22 + a03 * b32; double c10 = a10 * b00 + a11 * b10 + a12 * b20 + a13 * b30; double c11 = a10 * b01 + a11 * b11 + a12 * b21 + a13 * b31; double c12 = a10 * b02 + a11 * b12 + a12 * b22 + a13 * b32; double c20 = a20 * b00 + a21 * b10 + a22 * b20 + a23 * b30; double c21 = a20 * b01 + a21 * b11 + a22 * b21 + a23 * b31; double c22 = a20 * b02 + a21 * b12 + a22 * b22 + a23 * b32; r3_t c0 == (r3_t){c00,c01,c02}; r3_t c1 == (r3_t){c10,c11,c12}; r3_t c2 == (r3_t){c20,c21,c22}; { return LR3x3.T{c0,c1,c2}; } } /* END Mul_3x4_4x3 */ void SELF_DefVar(SELF_T *erg, bool_vec_t *variable) { /* Decide which cells are relevant to "Orientation" energy. A cell is relevant iff it exists,has at least one variable corner. */ int NV = erg->top->node.nel; ??? NP = erg->top->cell.nel; bool_vec_t vVar = erg->vVar^; ??? dell = erg->top->cell^; ??? polyRelevant = erg->polyRelevant^; { assert((variable.nel) == NV); vVar = variable; /* Find the relevant cells: */ for (i = 0; i < NP; i++){ polyRelevant[i] = FALSE; } for (i = 0; i < NP; i++) { ??? r = Srot(cell[i]); Place_t a = Tors(r); Place_t r1 = Clock(PrevE(r)); Place_t a1 = Tors(r1); Node_t u = OrgV(a); Node_t v = OrgV(NextE(a)); Node_t w = OrgV(PrevE(a)); Node_t x = OrgV(PrevE(a1)); ??? vvar = vVar.el[u->num]) || (vVar.el[v->num]) || (vVar.el[w->num]) || (vVar.el[x->num]; Wall_t f1 = PWall(a); Wall_t f2 = PWall(PrevF(a)); Wall_t f3 = PWall(PrevF(NextE(a))); Wall_t f4 = PWall(PrevF(PrevE(a))); { assert(f1 != f2) && (f2 != f3) && (f3 != f4); if (vvar) { polyRelevant[i] = TRUE; } } } } } /* END DefVar */ void SELF_Eval( SELF_T erg; Coords_t *c; double *e; bool_t grad; Gradient *eDc; ) { int NV = erg->top->node.nel; ??? NP = erg->top->cell.nel; ??? polyRelevant = erg->polyRelevant^; ??? cell = erg->top->cell^; bool_vec_t vVar = erg->vVar; ??? K = erg->K; ??? ep = erg->ep; ??? eDep = erg->eDep; ??? top = erg->top; ??? alpha = erg->alpha; ??? beta = erg->beta; ??? volA = (1.0/6.0) * LR3x3.Det(A); { LR4x3 ShapeMatrix(*uint i) LR3x4 b; VAR /* _ _ | v0-u0 w0-u0 x0-u0 | B == | v1-u1 w1-u1 x1-u1 | | v2-u2 w2-u2 x2-u2 | | v3-u3 w3-u3 x3-u3 | - - */ { ??? r = Srot(top.cell[i]); Place_t a = Tors(r); ??? tv = Triangulation.TetraNegNodes(a); ??? u = c[tv[0]->num]; ??? v = c[tv[1]->num]; ??? w = c[tv[2]->num]; ??? x = c[tv[3]->num]; { b[0] = r4_Sub(v,u); b[1] = r4_Sub(w, u); b[2] = r4_Sub(x,u); return Transpose_3x4(b); } } ShapeMatrix; double ElasticEnergy(uint i) /* This procedure compute the elastic energy of a tetrahedral element. Where Delta, Gamma and Sigma are the "invariants of rotation" inside of a tetrahedron, compute although the strain tensor ("C" matrix). The metric tensor is the "T" matrix, where T == Mul(Transpose(C), C) (the matrix T is simetric). From the metric tensor, we conpute the elastic energy density: "term1+term2". Its integral over the tetrahedron T is merely the product of the elastic energy density by its volume VT in its relaxed configuration (i.e. when the elasstic energy is zero). See the section 4.2 "Calculo das forcas de elasticidade" in the Msc. Thesis "Animacao Dinamica de Corpos Elasticos" by R.L.W.L */ { ??? B = ShapeMatrix(i); ??? C = Mul_4x3_3x3(B, LR3x3.Inv(A)); /* C == B * A_1 */ ??? T = Mul_3x4_4x3(Transpose_4x3(C), C); ??? T11 = T[1,1]*T[2,2] - T[1,2]*T[1,2]; ??? T12 = T[0,1]*T[2,2] - T[1,2]*T[0,2]; ??? T13 = T[0,1]*T[1,2] - T[1,1]*T[0,2]; ??? T22 = T[0,0]*T[2,2] - T[0,2]*T[0,2]; ??? T33 = T[0,0]*T[1,1] - T[0,1]*T[0,1]; ??? Delta = T[0,0]*T11 - T[0,1]*T12 + T[0,2]*T13; ??? Sigma = T11 + T22 + T33; ??? Gamma = T[0,0] + T[1,1] + T[2,2]; ??? term1 = alpha/32.0*(Delta*Delta+1.0/(Delta*Delta)-2.0); ??? term2 = beta/6.0*(Gamma*Gamma - 3.0*Sigma); ??? P = volA*(term1 + term2); { return P; } } ElasticEnergy; void Accum_ep(uint i; VAR epo: double) /* Compute the potential energy of the tetrahedron "i" in R^{3}. */ { epo = epo + ElasticEnergy(i); } Accum_ep; void Accum_e_from_ep(ep: double; double e; VAR VAR eDep: double) /* Adds to "e" the energy term corresponding to a cell with determinant "det". Also, if "grad" is true, stores in "eDdet" the derivative of that term. */ { e = e + K * ep; if (grad) { eDep = K; } else { eDep = 0.0; } } Accum_e_from_ep; void Distribute_eDep(uint i; eDep: double) /* Accumulates in "eDc" the gradient of "e" relative to the corners of the tetrahedron "i", given the derivative "eDep" of "e" relative to the tetrahedron's elastic energy "epv". */ { ??? B = ShapeMatrix(i); ??? ??? C = Mul_4x3_3x3(B, LR3x3.Inv(A)); ??? T = Mul_3x4_4x3(Transpose_4x3(C), C); double T11 = T[1,1]*T[2,2] - T[1,2]*T[1,2]; double T12 = T[0,1]*T[2,2] - T[1,2]*T[0,2]; double T13 = T[0,1]*T[1,2] - T[1,1]*T[0,2]; double T22 = T[0,0]*T[2,2] - T[0,2]*T[0,2]; double T33 = T[0,0]*T[1,1] - T[0,1]*T[0,1]; double Delta = T[0,0]*T11 - T[0,1]*T12 + T[0,2]*T13; double Sigma = T11 + T22 + T33; double Del_3 = 1.0/(Delta*Delta*Delta); ??? t = Srot(cell[i]); Place_t a = Tors(t); Place_t t1 = Clock(PrevE(t)); Place_t a1 = Tors(t1); Node_t iu = OrgV(a)->num; Node_t iv = OrgV(NextE(a))->num; Node_t iw = OrgV(PrevE(a))->num; Node_t ix = OrgV(PrevE(a1))->num; double t00 = T[0,0], t01 = T[0,1], t02 = T[0,2]; double t11 = T[1,1], t12 = T[1,2], t22 = T[2,2]; double c00 = C[0,0], c01 = C[0,1], c02 = C[0,2]; double c10 = C[1,0], c11 = C[1,1], c12 = C[1,2]; double c20 = C[2,0], c21 = C[2,1], c22 = C[2,2]; double c30 = C[3,0], c31 = C[3,1], c32 = C[3,2]; double b00Du0 = -1.0, b01Du0 = -1.0, b02Du0 = -1.0; double b10Du0 = 0.0, b11Du0 = 0.0, b12Du0 = 0.0; double b20Du0 = 0.0, b21Du0 = 0.0, b22Du0 = 0.0; double b30Du0 = 0.0, b31Du0 = 0.0, b32Du0 = 0.0; double b00Du1 = 0.0, b01Du1 = 0.0, b02Du1 = 0.0; double b10Du1 = -1.0, b11Du1 = -1.0, b12Du1 = -1.0; double b20Du1 = 0.0, b21Du1 = 0.0, b22Du1 = 0.0; double b30Du1 = 0.0, b31Du1 = 0.0, b32Du1 = 0.0; double b00Du2 = 0.0, b01Du2 = 0.0, b02Du2 = 0.0; double b10Du2 = 0.0, b11Du2 = 0.0, b12Du2 = 0.0; double b20Du2 = -1.0, b21Du2 = -1.0, b22Du2 = -1.0; double b30Du2 = 0.0, b31Du2 = 0.0, b32Du2 = 0.0; double b00Du3 = 0.0, b01Du3 = 0.0, b02Du3 = 0.0; double b10Du3 = 0.0, b11Du3 = 0.0, b12Du3 = 0.0; double b20Du3 = 0.0, b21Du3 = 0.0, b22Du3 = 0.0; double b30Du3 = -1.0, b31Du3 = -1.0, b32Du3 = -1.0; double b00Dv0 = 1.0, b01Dv0 = 0.0, b02Dv0 = 0.0; double b10Dv0 = 0.0, b11Dv0 = 0.0, b12Dv0 = 0.0; double b20Dv0 = 0.0, b21Dv0 = 0.0, b22Dv0 = 0.0; double b30Dv0 = 0.0, b31Dv0 = 0.0, b32Dv0 = 0.0; double b00Dv1 = 0.0, b01Dv1 = 0.0, b02Dv1 = 0.0; double b10Dv1 = 1.0, b11Dv1 = 0.0, b12Dv1 = 0.0; double b20Dv1 = 0.0, b21Dv1 = 0.0, b22Dv1 = 0.0; double b30Dv1 = 0.0, b31Dv1 = 0.0, b32Dv1 = 0.0; double b00Dv2 = 0.0, b01Dv2 = 0.0, b02Dv2 = 0.0; double b10Dv2 = 0.0, b11Dv2 = 0.0, b12Dv2 = 0.0; double b20Dv2 = 1.0, b21Dv2 = 0.0, b22Dv2 = 0.0; double b30Dv2 = 0.0, b31Dv2 = 0.0, b32Dv2 = 0.0; double b00Dv3 = 0.0, b01Dv3 = 0.0, b02Dv3 = 0.0; double b10Dv3 = 0.0, b11Dv3 = 0.0, b12Dv3 = 0.0; double b20Dv3 = 0.0, b21Dv3 = 0.0, b22Dv3 = 0.0; double b30Dv3 = 1.0, b31Dv3 = 0.0, b32Dv3 = 0.0; double b00Dw0 = 0.0, b01Dw0 = 1.0, b02Dw0 = 0.0; double b10Dw0 = 0.0, b11Dw0 = 0.0, b12Dw0 = 0.0; double b20Dw0 = 0.0, b21Dw0 = 0.0, b22Dw0 = 0.0; double b30Dw0 = 0.0, b31Dw0 = 0.0, b32Dw0 = 0.0; double b00Dw1 = 0.0, b01Dw1 = 0.0, b02Dw1 = 0.0; double b10Dw1 = 0.0, b11Dw1 = 1.0, b12Dw1 = 0.0; double b20Dw1 = 0.0, b21Dw1 = 0.0, b22Dw1 = 0.0; double b30Dw1 = 0.0, b31Dw1 = 0.0, b32Dw1 = 0.0; double b00Dw2 = 0.0, b01Dw2 = 0.0, b02Dw2 = 0.0; double b10Dw2 = 0.0, b11Dw2 = 0.0, b12Dw2 = 0.0; double b20Dw2 = 0.0, b21Dw2 = 1.0, b22Dw2 = 0.0; double b30Dw2 = 0.0, b31Dw2 = 0.0, b32Dw2 = 0.0; double b00Dw3 = 0.0, b01Dw3 = 0.0, b02Dw3 = 0.0; double b10Dw3 = 0.0, b11Dw3 = 0.0, b12Dw3 = 0.0; double b20Dw3 = 0.0, b21Dw3 = 0.0, b22Dw3 = 0.0; double b30Dw3 = 0.0, b31Dw3 = 1.0, b32Dw3 = 0.0; double b00Dx0 = 0.0, b01Dx0 = 0.0, b02Dx0 = 1.0; double b10Dx0 = 0.0, b11Dx0 = 0.0, b12Dx0 = 0.0; double b20Dx0 = 0.0, b21Dx0 = 0.0, b22Dx0 = 0.0; double b30Dx0 = 0.0, b31Dx0 = 0.0, b32Dx0 = 0.0; double b00Dx1 = 0.0, b01Dx1 = 0.0, b02Dx1 = 0.0; double b10Dx1 = 0.0, b11Dx1 = 0.0, b12Dx1 = 1.0; double b20Dx1 = 0.0, b21Dx1 = 0.0, b22Dx1 = 0.0; double b30Dx1 = 0.0, b31Dx1 = 0.0, b32Dx1 = 0.0; double b00Dx2 = 0.0, b01Dx2 = 0.0, b02Dx2 = 0.0; double b10Dx2 = 0.0, b11Dx2 = 0.0, b12Dx2 = 0.0; double b20Dx2 = 0.0, b21Dx2 = 0.0, b22Dx2 = 1.0; double b30Dx2 = 0.0, b31Dx2 = 0.0, b32Dx2 = 0.0; double b00Dx3 = 0.0, b01Dx3 = 0.0, b02Dx3 = 0.0; double b10Dx3 = 0.0, b11Dx3 = 0.0, b12Dx3 = 0.0; double b20Dx3 = 0.0, b21Dx3 = 0.0, b22Dx3 = 0.0; double b30Dx3 = 0.0, b31Dx3 = 0.0, b32Dx3 = 1.0; double c00Du0 = 0.5 * (-b00Du0 + b01Du0 + b02Du0); double c01Du0 = 0.5 * ( b00Du0 - b01Du0 + b02Du0); double c02Du0 = 0.5 * ( b00Du0 + b01Du0 - b02Du0); double c10Du0 = 0.5 * (-b10Du0 + b11Du0 + b12Du0); double c11Du0 = 0.5 * ( b10Du0 - b11Du0 + b12Du0); double c12Du0 = 0.5 * ( b10Du0 + b11Du0 - b12Du0); double c20Du0 = 0.5 * (-b20Du0 + b21Du0 + b22Du0); double c21Du0 = 0.5 * ( b20Du0 - b21Du0 + b22Du0); double c22Du0 = 0.5 * ( b20Du0 + b21Du0 - b22Du0); double c30Du0 = 0.5 * (-b30Du0 + b31Du0 + b32Du0); double c31Du0 = 0.5 * ( b30Du0 - b31Du0 + b32Du0); double c32Du0 = 0.5 * ( b30Du0 + b31Du0 - b32Du0); double c00Du1 = 0.5 * (-b00Du1 + b01Du1 + b02Du1); double c01Du1 = 0.5 * ( b00Du1 - b01Du1 + b02Du1); double c02Du1 = 0.5 * ( b00Du1 + b01Du1 - b02Du1); double c10Du1 = 0.5 * (-b10Du1 + b11Du1 + b12Du1); double c11Du1 = 0.5 * ( b10Du1 - b11Du1 + b12Du1); double c12Du1 = 0.5 * ( b10Du1 + b11Du1 - b12Du1); double c20Du1 = 0.5 * (-b20Du1 + b21Du1 + b22Du1); double c21Du1 = 0.5 * ( b20Du1 - b21Du1 + b22Du1); double c22Du1 = 0.5 * ( b20Du1 + b21Du1 - b22Du1); double c30Du1 = 0.5 * (-b30Du1 + b31Du1 + b32Du1); double c31Du1 = 0.5 * ( b30Du1 - b31Du1 + b32Du1); double c32Du1 = 0.5 * ( b30Du1 + b31Du1 - b32Du1); double c00Du2 = 0.5 * (-b00Du2 + b01Du2 + b02Du2); double c01Du2 = 0.5 * ( b00Du2 - b01Du2 + b02Du2); double c02Du2 = 0.5 * ( b00Du2 + b01Du2 - b02Du2); double c10Du2 = 0.5 * (-b10Du2 + b11Du2 + b12Du2); double c11Du2 = 0.5 * ( b10Du2 - b11Du2 + b12Du2); double c12Du2 = 0.5 * ( b10Du2 + b11Du2 - b12Du2); double c20Du2 = 0.5 * (-b20Du2 + b21Du2 + b22Du2); double c21Du2 = 0.5 * ( b20Du2 - b21Du2 + b22Du2); double c22Du2 = 0.5 * ( b20Du2 + b21Du2 - b22Du2); double c30Du2 = 0.5 * (-b30Du2 + b31Du2 + b32Du2); double c31Du2 = 0.5 * ( b30Du2 - b31Du2 + b32Du2); double c32Du2 = 0.5 * ( b30Du2 + b31Du2 - b32Du2); double c00Du3 = 0.5 * (-b00Du3 + b01Du3 + b02Du3); double c01Du3 = 0.5 * ( b00Du3 - b01Du3 + b02Du3); double c02Du3 = 0.5 * ( b00Du3 + b01Du3 - b02Du3); double c10Du3 = 0.5 * (-b10Du3 + b11Du3 + b12Du3); double c11Du3 = 0.5 * ( b10Du3 - b11Du3 + b12Du3); double c12Du3 = 0.5 * ( b10Du3 + b11Du3 - b12Du3); double c20Du3 = 0.5 * (-b20Du3 + b21Du3 + b22Du3); double c21Du3 = 0.5 * ( b20Du3 - b21Du3 + b22Du3); double c22Du3 = 0.5 * ( b20Du3 + b21Du3 - b22Du3); double c30Du3 = 0.5 * (-b30Du3 + b31Du3 + b32Du3); double c31Du3 = 0.5 * ( b30Du3 - b31Du3 + b32Du3); double c32Du3 = 0.5 * ( b30Du3 + b31Du3 - b32Du3); double c00Dv0 = 0.5 * (-b00Dv0 + b01Dv0 + b02Dv0); double c01Dv0 = 0.5 * ( b00Dv0 - b01Dv0 + b02Dv0); double c02Dv0 = 0.5 * ( b00Dv0 + b01Dv0 - b02Dv0); double c10Dv0 = 0.5 * (-b10Dv0 + b11Dv0 + b12Dv0); double c11Dv0 = 0.5 * ( b10Dv0 - b11Dv0 + b12Dv0); double c12Dv0 = 0.5 * ( b10Dv0 + b11Dv0 - b12Dv0); double c20Dv0 = 0.5 * (-b20Dv0 + b21Dv0 + b22Dv0); double c21Dv0 = 0.5 * ( b20Dv0 - b21Dv0 + b22Dv0); double c22Dv0 = 0.5 * ( b20Dv0 + b21Dv0 - b22Dv0); double c30Dv0 = 0.5 * (-b30Dv0 + b31Dv0 + b32Dv0); double c31Dv0 = 0.5 * ( b30Dv0 - b31Dv0 + b32Dv0); double c32Dv0 = 0.5 * ( b30Dv0 + b31Dv0 - b32Dv0); double c00Dv1 = 0.5 * (-b00Dv1 + b01Dv1 + b02Dv1); double c01Dv1 = 0.5 * ( b00Dv1 - b01Dv1 + b02Dv1); double c02Dv1 = 0.5 * ( b00Dv1 + b01Dv1 - b02Dv1); double c10Dv1 = 0.5 * (-b10Dv1 + b11Dv1 + b12Dv1); double c11Dv1 = 0.5 * ( b10Dv1 - b11Dv1 + b12Dv1); double c12Dv1 = 0.5 * ( b10Dv1 + b11Dv1 - b12Dv1); double c20Dv1 = 0.5 * (-b20Dv1 + b21Dv1 + b22Dv1); double c21Dv1 = 0.5 * ( b20Dv1 - b21Dv1 + b22Dv1); double c22Dv1 = 0.5 * ( b20Dv1 + b21Dv1 - b22Dv1); double c30Dv1 = 0.5 * (-b30Dv1 + b31Dv1 + b32Dv1); double c31Dv1 = 0.5 * ( b30Dv1 - b31Dv1 + b32Dv1); double c32Dv1 = 0.5 * ( b30Dv1 + b31Dv1 - b32Dv1); double c00Dv2 = 0.5 * (-b00Dv2 + b01Dv2 + b02Dv2); double c01Dv2 = 0.5 * ( b00Dv2 - b01Dv2 + b02Dv2); double c02Dv2 = 0.5 * ( b00Dv2 + b01Dv2 - b02Dv2); double c10Dv2 = 0.5 * (-b10Dv2 + b11Dv2 + b12Dv2); double c11Dv2 = 0.5 * ( b10Dv2 - b11Dv2 + b12Dv2); double c12Dv2 = 0.5 * ( b10Dv2 + b11Dv2 - b12Dv2); double c20Dv2 = 0.5 * (-b20Dv2 + b21Dv2 + b22Dv2); double c21Dv2 = 0.5 * ( b20Dv2 - b21Dv2 + b22Dv2); double c22Dv2 = 0.5 * ( b20Dv2 + b21Dv2 - b22Dv2); double c30Dv2 = 0.5 * (-b30Dv2 + b31Dv2 + b32Dv2); double c31Dv2 = 0.5 * ( b30Dv2 - b31Dv2 + b32Dv2); double c32Dv2 = 0.5 * ( b30Dv2 + b31Dv2 - b32Dv2); double c00Dv3 = 0.5 * (-b00Dv3 + b01Dv3 + b02Dv3); double c01Dv3 = 0.5 * ( b00Dv3 - b01Dv3 + b02Dv3); double c02Dv3 = 0.5 * ( b00Dv3 + b01Dv3 - b02Dv3); double c10Dv3 = 0.5 * (-b10Dv3 + b11Dv3 + b12Dv3); double c11Dv3 = 0.5 * ( b10Dv3 - b11Dv3 + b12Dv3); double c12Dv3 = 0.5 * ( b10Dv3 + b11Dv3 - b12Dv3); double c20Dv3 = 0.5 * (-b20Dv3 + b21Dv3 + b22Dv3); double c21Dv3 = 0.5 * ( b20Dv3 - b21Dv3 + b22Dv3); double c22Dv3 = 0.5 * ( b20Dv3 + b21Dv3 - b22Dv3); double c30Dv3 = 0.5 * (-b30Dv3 + b31Dv3 + b32Dv3); double c31Dv3 = 0.5 * ( b30Dv3 - b31Dv3 + b32Dv3); double c32Dv3 = 0.5 * ( b30Dv3 + b31Dv3 - b32Dv3); double c00Dw0 = 0.5 * (-b00Dw0 + b01Dw0 + b02Dw0); double c01Dw0 = 0.5 * ( b00Dw0 - b01Dw0 + b02Dw0); double c02Dw0 = 0.5 * ( b00Dw0 + b01Dw0 - b02Dw0); double c10Dw0 = 0.5 * (-b10Dw0 + b11Dw0 + b12Dw0); double c11Dw0 = 0.5 * ( b10Dw0 - b11Dw0 + b12Dw0); double c12Dw0 = 0.5 * ( b10Dw0 + b11Dw0 - b12Dw0); double c20Dw0 = 0.5 * (-b20Dw0 + b21Dw0 + b22Dw0); double c21Dw0 = 0.5 * ( b20Dw0 - b21Dw0 + b22Dw0); double c22Dw0 = 0.5 * ( b20Dw0 + b21Dw0 - b22Dw0); double c30Dw0 = 0.5 * (-b30Dw0 + b31Dw0 + b32Dw0); double c31Dw0 = 0.5 * ( b30Dw0 - b31Dw0 + b32Dw0); double c32Dw0 = 0.5 * ( b30Dw0 + b31Dw0 - b32Dw0); double c00Dw1 = 0.5 * (-b00Dw1 + b01Dw1 + b02Dw1); double c01Dw1 = 0.5 * ( b00Dw1 - b01Dw1 + b02Dw1); double c02Dw1 = 0.5 * ( b00Dw1 + b01Dw1 - b02Dw1); double c10Dw1 = 0.5 * (-b10Dw1 + b11Dw1 + b12Dw1); double c11Dw1 = 0.5 * ( b10Dw1 - b11Dw1 + b12Dw1); double c12Dw1 = 0.5 * ( b10Dw1 + b11Dw1 - b12Dw1); double c20Dw1 = 0.5 * (-b20Dw1 + b21Dw1 + b22Dw1); double c21Dw1 = 0.5 * ( b20Dw1 - b21Dw1 + b22Dw1); double c22Dw1 = 0.5 * ( b20Dw1 + b21Dw1 - b22Dw1); double c30Dw1 = 0.5 * (-b30Dw1 + b31Dw1 + b32Dw1); double c31Dw1 = 0.5 * ( b30Dw1 - b31Dw1 + b32Dw1); double c32Dw1 = 0.5 * ( b30Dw1 + b31Dw1 - b32Dw1); double c00Dw2 = 0.5 * (-b00Dw2 + b01Dw2 + b02Dw2); double c01Dw2 = 0.5 * ( b00Dw2 - b01Dw2 + b02Dw2); double c02Dw2 = 0.5 * ( b00Dw2 + b01Dw2 - b02Dw2); double c10Dw2 = 0.5 * (-b10Dw2 + b11Dw2 + b12Dw2); double c11Dw2 = 0.5 * ( b10Dw2 - b11Dw2 + b12Dw2); double c12Dw2 = 0.5 * ( b10Dw2 + b11Dw2 - b12Dw2); double c20Dw2 = 0.5 * (-b20Dw2 + b21Dw2 + b22Dw2); double c21Dw2 = 0.5 * ( b20Dw2 - b21Dw2 + b22Dw2); double c22Dw2 = 0.5 * ( b20Dw2 + b21Dw2 - b22Dw2); double c30Dw2 = 0.5 * (-b30Dw2 + b31Dw2 + b32Dw2); double c31Dw2 = 0.5 * ( b30Dw2 - b31Dw2 + b32Dw2); double c32Dw2 = 0.5 * ( b30Dw2 + b31Dw2 - b32Dw2); double c00Dw3 = 0.5 * (-b00Dw3 + b01Dw3 + b02Dw3); double c01Dw3 = 0.5 * ( b00Dw3 - b01Dw3 + b02Dw3); double c02Dw3 = 0.5 * ( b00Dw3 + b01Dw3 - b02Dw3); double c10Dw3 = 0.5 * (-b10Dw3 + b11Dw3 + b12Dw3); double c11Dw3 = 0.5 * ( b10Dw3 - b11Dw3 + b12Dw3); double c12Dw3 = 0.5 * ( b10Dw3 + b11Dw3 - b12Dw3); double c20Dw3 = 0.5 * (-b20Dw3 + b21Dw3 + b22Dw3); double c21Dw3 = 0.5 * ( b20Dw3 - b21Dw3 + b22Dw3); double c22Dw3 = 0.5 * ( b20Dw3 + b21Dw3 - b22Dw3); double c30Dw3 = 0.5 * (-b30Dw3 + b31Dw3 + b32Dw3); double c31Dw3 = 0.5 * ( b30Dw3 - b31Dw3 + b32Dw3); double c32Dw3 = 0.5 * ( b30Dw3 + b31Dw3 - b32Dw3); double c00Dx0 = 0.5 * (-b00Dx0 + b01Dx0 + b02Dx0); double c01Dx0 = 0.5 * ( b00Dx0 - b01Dx0 + b02Dx0); double c02Dx0 = 0.5 * ( b00Dx0 + b01Dx0 - b02Dx0); double c10Dx0 = 0.5 * (-b10Dx0 + b11Dx0 + b12Dx0); double c11Dx0 = 0.5 * ( b10Dx0 - b11Dx0 + b12Dx0); double c12Dx0 = 0.5 * ( b10Dx0 + b11Dx0 - b12Dx0); double c20Dx0 = 0.5 * (-b20Dx0 + b21Dx0 + b22Dx0); double c21Dx0 = 0.5 * ( b20Dx0 - b21Dx0 + b22Dx0); double c22Dx0 = 0.5 * ( b20Dx0 + b21Dx0 - b22Dx0); double c30Dx0 = 0.5 * (-b30Dx0 + b31Dx0 + b32Dx0); double c31Dx0 = 0.5 * ( b30Dx0 - b31Dx0 + b32Dx0); double c32Dx0 = 0.5 * ( b30Dx0 + b31Dx0 - b32Dx0); double c00Dx1 = 0.5 * (-b00Dx1 + b01Dx1 + b02Dx1); double c01Dx1 = 0.5 * ( b00Dx1 - b01Dx1 + b02Dx1); double c02Dx1 = 0.5 * ( b00Dx1 + b01Dx1 - b02Dx1); double c10Dx1 = 0.5 * (-b10Dx1 + b11Dx1 + b12Dx1); double c11Dx1 = 0.5 * ( b10Dx1 - b11Dx1 + b12Dx1); double c12Dx1 = 0.5 * ( b10Dx1 + b11Dx1 - b12Dx1); double c20Dx1 = 0.5 * (-b20Dx1 + b21Dx1 + b22Dx1); double c21Dx1 = 0.5 * ( b20Dx1 - b21Dx1 + b22Dx1); double c22Dx1 = 0.5 * ( b20Dx1 + b21Dx1 - b22Dx1); double c30Dx1 = 0.5 * (-b30Dx1 + b31Dx1 + b32Dx1); double c31Dx1 = 0.5 * ( b30Dx1 - b31Dx1 + b32Dx1); double c32Dx1 = 0.5 * ( b30Dx1 + b31Dx1 - b32Dx1); double c00Dx2 = 0.5 * (-b00Dx2 + b01Dx2 + b02Dx2); double c01Dx2 = 0.5 * ( b00Dx2 - b01Dx2 + b02Dx2); double c02Dx2 = 0.5 * ( b00Dx2 + b01Dx2 - b02Dx2); double c10Dx2 = 0.5 * (-b10Dx2 + b11Dx2 + b12Dx2); double c11Dx2 = 0.5 * ( b10Dx2 - b11Dx2 + b12Dx2); double c12Dx2 = 0.5 * ( b10Dx2 + b11Dx2 - b12Dx2); double c20Dx2 = 0.5 * (-b20Dx2 + b21Dx2 + b22Dx2); double c21Dx2 = 0.5 * ( b20Dx2 - b21Dx2 + b22Dx2); double c22Dx2 = 0.5 * ( b20Dx2 + b21Dx2 - b22Dx2); double c30Dx2 = 0.5 * (-b30Dx2 + b31Dx2 + b32Dx2); double c31Dx2 = 0.5 * ( b30Dx2 - b31Dx2 + b32Dx2); double c32Dx2 = 0.5 * ( b30Dx2 + b31Dx2 - b32Dx2); double c00Dx3 = 0.5 * (-b00Dx3 + b01Dx3 + b02Dx3); double c01Dx3 = 0.5 * ( b00Dx3 - b01Dx3 + b02Dx3); double c02Dx3 = 0.5 * ( b00Dx3 + b01Dx3 - b02Dx3); double c10Dx3 = 0.5 * (-b10Dx3 + b11Dx3 + b12Dx3); double c11Dx3 = 0.5 * ( b10Dx3 - b11Dx3 + b12Dx3); double c12Dx3 = 0.5 * ( b10Dx3 + b11Dx3 - b12Dx3); double c20Dx3 = 0.5 * (-b20Dx3 + b21Dx3 + b22Dx3); double c21Dx3 = 0.5 * ( b20Dx3 - b21Dx3 + b22Dx3); double c22Dx3 = 0.5 * ( b20Dx3 + b21Dx3 - b22Dx3); double c30Dx3 = 0.5 * (-b30Dx3 + b31Dx3 + b32Dx3); double c31Dx3 = 0.5 * ( b30Dx3 - b31Dx3 + b32Dx3); double c32Dx3 = 0.5 * ( b30Dx3 + b31Dx3 - b32Dx3); double t00Du0 = 2.0 * (c00*c00Du0+c10*c10Du0+c20*c20Du0+c30*c30Du0); double t10Du0 = (c00*c01Du0+c01*c00Du0) + (c10*c11Du0+c11*c10Du0) + (c20*c21Du0+c21*c20Du0) + (c30*c31Du0+c31*c30Du0); double t11Du0 = 2.0 * (c01*c01Du0+c11*c11Du0+c21*c21Du0+c31*c31Du0); double t20Du0 = (c00*c02Du0+c02*c00Du0) + (c10*c12Du0+c12*c10Du0) + (c20*c22Du0+c22*c20Du0) + (c30*c32Du0+c32*c30Du0); double t21Du0 = (c01*c00Du0+c00*c01Du0) + (c11*c10Du0+c10*c11Du0) + (c21*c20Du0+c20*c21Du0) + (c31*c30Du0+c30*c31Du0); double t22Du0 = 2.0*(c02*c02Du0+c12*c12Du0+c22*c22Du0+c32*c32Du0); double t01Du0 = t10Du0; double t02Du0 = t20Du0; double t12Du0 = t21Du0; double t00Du1 = 2.0 * (c00*c00Du1+c10*c10Du1+c20*c20Du1+c30*c30Du1); double t10Du1 = (c00*c01Du1+c01*c00Du1) + (c10*c11Du1+c11*c10Du1) + (c20*c21Du1+c21*c20Du1) + (c30*c31Du1+c31*c30Du1); double t11Du1 = 2.0 * (c01*c01Du1+c11*c11Du1+c21*c21Du1+c31*c31Du1); double t20Du1 = (c00*c02Du1+c02*c00Du1) + (c10*c12Du1+c12*c10Du1) + (c20*c22Du1+c22*c20Du1) + (c30*c32Du1+c32*c30Du1); double t21Du1 = (c01*c00Du1+c00*c01Du1) + (c11*c10Du1+c10*c11Du1) + (c21*c20Du1+c20*c21Du1) + (c31*c30Du1+c30*c31Du1); double t22Du1 = 2.0*(c02*c02Du1+c12*c12Du1+c22*c22Du1+c32*c32Du1); double t01Du1 = t10Du1; double t02Du1 = t20Du1; double t12Du1 = t21Du1; double t00Du2 = 2.0 * (c00*c00Du2+c10*c10Du2+c20*c20Du2+c30*c30Du2); double t10Du2 = (c00*c01Du2+c01*c00Du2) + (c10*c11Du2+c11*c10Du2) + (c20*c21Du2+c21*c20Du2) + (c30*c31Du2+c31*c30Du2); double t11Du2 = 2.0 * (c01*c01Du2+c11*c11Du2+c21*c21Du2+c31*c31Du2); double t20Du2 = (c00*c02Du2+c02*c00Du2) + (c10*c12Du2+c12*c10Du2) + (c20*c22Du2+c22*c20Du2) + (c30*c32Du2+c32*c30Du2); double t21Du2 = (c01*c00Du2+c00*c01Du2) + (c11*c10Du2+c10*c11Du2) + (c21*c20Du2+c20*c21Du2) + (c31*c30Du2+c30*c31Du2); double t22Du2 = 2.0*(c02*c02Du2+c12*c12Du2+c22*c22Du2+c32*c32Du2); double t01Du2 = t10Du2; double t02Du2 = t20Du2; double t12Du2 = t21Du2; double t00Du3 = 2.0 * (c00*c00Du3+c10*c10Du3+c20*c20Du3+c30*c30Du3); double t10Du3 = (c00*c01Du3+c01*c00Du3) + (c10*c11Du3+c11*c10Du3) + (c20*c21Du3+c21*c20Du3) + (c30*c31Du3+c31*c30Du3); double t11Du3 = 2.0 * (c01*c01Du3+c11*c11Du3+c21*c21Du3+c31*c31Du3); double t20Du3 = (c00*c02Du3+c02*c00Du3) + (c10*c12Du3+c12*c10Du3) + (c20*c22Du3+c22*c20Du3) + (c30*c32Du3+c32*c30Du3); double t21Du3 = (c01*c00Du3+c00*c01Du3) + (c11*c10Du3+c10*c11Du3) + (c21*c20Du3+c20*c21Du3) + (c31*c30Du3+c30*c31Du3); double t22Du3 = 2.0*(c02*c02Du3+c12*c12Du3+c22*c22Du3+c32*c32Du3); double t01Du3 = t10Du3; double t02Du3 = t20Du3; double t12Du3 = t21Du3; double t00Dv0 = 2.0 * (c00*c00Dv0+c10*c10Dv0+c20*c20Dv0+c30*c30Dv0); double t10Dv0 = (c00*c01Dv0+c01*c00Dv0) + (c10*c11Dv0+c11*c10Dv0) + (c20*c21Dv0+c21*c20Dv0) + (c30*c31Dv0+c31*c30Dv0); double t11Dv0 = 2.0 * (c01*c01Dv0+c11*c11Dv0+c21*c21Dv0+c31*c31Dv0); double t20Dv0 = (c00*c02Dv0+c02*c00Dv0) + (c10*c12Dv0+c12*c10Dv0) + (c20*c22Dv0+c22*c20Dv0) + (c30*c32Dv0+c32*c30Dv0); double t21Dv0 = (c01*c00Dv0+c00*c01Dv0) + (c11*c10Dv0+c10*c11Dv0) + (c21*c20Dv0+c20*c21Dv0) + (c31*c30Dv0+c30*c31Dv0); double t22Dv0 = 2.0*(c02*c02Dv0+c12*c12Dv0+c22*c22Dv0+c32*c32Dv0); double t01Dv0 = t10Dv0; double t02Dv0 = t20Dv0; double t12Dv0 = t21Dv0; double t00Dv1 = 2.0 * (c00*c00Dv1+c10*c10Dv1+c20*c20Dv1+c30*c30Dv1); double t10Dv1 = (c00*c01Dv1+c01*c00Dv1) + (c10*c11Dv1+c11*c10Dv1) + (c20*c21Dv1+c21*c20Dv1) + (c30*c31Dv1+c31*c30Dv1); double t11Dv1 = 2.0 * (c01*c01Dv1+c11*c11Dv1+c21*c21Dv1+c31*c31Dv1); double t20Dv1 = (c00*c02Dv1+c02*c00Dv1) + (c10*c12Dv1+c12*c10Dv1) + (c20*c22Dv1+c22*c20Dv1) + (c30*c32Dv1+c32*c30Dv1); double t21Dv1 = (c01*c00Dv1+c00*c01Dv1) + (c11*c10Dv1+c10*c11Dv1) + (c21*c20Dv1+c20*c21Dv1) + (c31*c30Dv1+c30*c31Dv1); double t22Dv1 = 2.0*(c02*c02Dv1+c12*c12Dv1+c22*c22Dv1+c32*c32Dv1); double t01Dv1 = t10Dv1; double t02Dv1 = t20Dv1; double t12Dv1 = t21Dv1; double t00Dv2 = 2.0 * (c00*c00Dv2+c10*c10Dv2+c20*c20Dv2+c30*c30Dv2); double t10Dv2 = (c00*c01Dv2+c01*c00Dv2) + (c10*c11Dv2+c11*c10Dv2) + (c20*c21Dv2+c21*c20Dv2) + (c30*c31Dv2+c31*c30Dv2); double t11Dv2 = 2.0 * (c01*c01Dv2+c11*c11Dv2+c21*c21Dv2+c31*c31Dv2); double t20Dv2 = (c00*c02Dv2+c02*c00Dv2) + (c10*c12Dv2+c12*c10Dv2) + (c20*c22Dv2+c22*c20Dv2) + (c30*c32Dv2+c32*c30Dv2); double t21Dv2 = (c01*c00Dv2+c00*c01Dv2) + (c11*c10Dv2+c10*c11Dv2) + (c21*c20Dv2+c20*c21Dv2) + (c31*c30Dv2+c30*c31Dv2); double t22Dv2 = 2.0*(c02*c02Dv2+c12*c12Dv2+c22*c22Dv2+c32*c32Dv2); double t01Dv2 = t10Dv2; double t02Dv2 = t20Dv2; double t12Dv2 = t21Dv2; double t00Dv3 = 2.0 * (c00*c00Dv3+c10*c10Dv3+c20*c20Dv3+c30*c30Dv3); double t10Dv3 = (c00*c01Dv3+c01*c00Dv3) + (c10*c11Dv3+c11*c10Dv3) + (c20*c21Dv3+c21*c20Dv3) + (c30*c31Dv3+c31*c30Dv3); double t11Dv3 = 2.0 * (c01*c01Dv3+c11*c11Dv3+c21*c21Dv3+c31*c31Dv3); double t20Dv3 = (c00*c02Dv3+c02*c00Dv3) + (c10*c12Dv3+c12*c10Dv3) + (c20*c22Dv3+c22*c20Dv3) + (c30*c32Dv3+c32*c30Dv3); double t21Dv3 = (c01*c00Dv3+c00*c01Dv3) + (c11*c10Dv3+c10*c11Dv3) + (c21*c20Dv3+c20*c21Dv3) + (c31*c30Dv3+c30*c31Dv3); double t22Dv3 = 2.0*(c02*c02Dv3+c12*c12Dv3+c22*c22Dv3+c32*c32Dv3); double t01Dv3 = t10Dv3; double t02Dv3 = t20Dv3; double t12Dv3 = t21Dv3; double t00Dw0 = 2.0 * (c00*c00Dw0+c10*c10Dw0+c20*c20Dw0+c30*c30Dw0); double t10Dw0 = (c00*c01Dw0+c01*c00Dw0) + (c10*c11Dw0+c11*c10Dw0) + (c20*c21Dw0+c21*c20Dw0) + (c30*c31Dw0+c31*c30Dw0); double t11Dw0 = 2.0 * (c01*c01Dw0+c11*c11Dw0+c21*c21Dw0+c31*c31Dw0); double t20Dw0 = (c00*c02Dw0+c02*c00Dw0) + (c10*c12Dw0+c12*c10Dw0) + (c20*c22Dw0+c22*c20Dw0) + (c30*c32Dw0+c32*c30Dw0); double t21Dw0 = (c01*c00Dw0+c00*c01Dw0) + (c11*c10Dw0+c10*c11Dw0) + (c21*c20Dw0+c20*c21Dw0) + (c31*c30Dw0+c30*c31Dw0); double t22Dw0 = 2.0*(c02*c02Dw0+c12*c12Dw0+c22*c22Dw0+c32*c32Dw0); double t01Dw0 = t10Dw0; double t02Dw0 = t20Dw0; double t12Dw0 = t21Dw0; double t00Dw1 = 2.0 * (c00*c00Dw1+c10*c10Dw1+c20*c20Dw1+c30*c30Dw1); double t10Dw1 = (c00*c01Dw1+c01*c00Dw1) + (c10*c11Dw1+c11*c10Dw1) + (c20*c21Dw1+c21*c20Dw1) + (c30*c31Dw1+c31*c30Dw1); double t11Dw1 = 2.0 * (c01*c01Dw1+c11*c11Dw1+c21*c21Dw1+c31*c31Dw1); double t20Dw1 = (c00*c02Dw1+c02*c00Dw1) + (c10*c12Dw1+c12*c10Dw1) + (c20*c22Dw1+c22*c20Dw1) + (c30*c32Dw1+c32*c30Dw1); double t21Dw1 = (c01*c00Dw1+c00*c01Dw1) + (c11*c10Dw1+c10*c11Dw1) + (c21*c20Dw1+c20*c21Dw1) + (c31*c30Dw1+c30*c31Dw1); double t22Dw1 = 2.0*(c02*c02Dw1+c12*c12Dw1+c22*c22Dw1+c32*c32Dw1); double t01Dw1 = t10Dw1; double t02Dw1 = t20Dw1; double t12Dw1 = t21Dw1; double t00Dw2 = 2.0 * (c00*c00Dw2+c10*c10Dw2+c20*c20Dw2+c30*c30Dw2); double t10Dw2 = (c00*c01Dw2+c01*c00Dw2) + (c10*c11Dw2+c11*c10Dw2) + (c20*c21Dw2+c21*c20Dw2) + (c30*c31Dw2+c31*c30Dw2); double t11Dw2 = 2.0 * (c01*c01Dw2+c11*c11Dw2+c21*c21Dw2+c31*c31Dw2); double t20Dw2 = (c00*c02Dw2+c02*c00Dw2) + (c10*c12Dw2+c12*c10Dw2) + (c20*c22Dw2+c22*c20Dw2) + (c30*c32Dw2+c32*c30Dw2); double t21Dw2 = (c01*c00Dw2+c00*c01Dw2) + (c11*c10Dw2+c10*c11Dw2) + (c21*c20Dw2+c20*c21Dw2) + (c31*c30Dw2+c30*c31Dw2); double t22Dw2 = 2.0*(c02*c02Dw2+c12*c12Dw2+c22*c22Dw2+c32*c32Dw2); double t01Dw2 = t10Dw2; double t02Dw2 = t20Dw2; double t12Dw2 = t21Dw2; double t00Dw3 = 2.0 * (c00*c00Dw3+c10*c10Dw3+c20*c20Dw3+c30*c30Dw3); double t10Dw3 = (c00*c01Dw3+c01*c00Dw3) + (c10*c11Dw3+c11*c10Dw3) + (c20*c21Dw3+c21*c20Dw3) + (c30*c31Dw3+c31*c30Dw3); double t11Dw3 = 2.0 * (c01*c01Dw3+c11*c11Dw3+c21*c21Dw3+c31*c31Dw3); double t20Dw3 = (c00*c02Dw3+c02*c00Dw3) + (c10*c12Dw3+c12*c10Dw3) + (c20*c22Dw3+c22*c20Dw3) + (c30*c32Dw3+c32*c30Dw3); double t21Dw3 = (c01*c00Dw3+c00*c01Dw3) + (c11*c10Dw3+c10*c11Dw3) + (c21*c20Dw3+c20*c21Dw3) + (c31*c30Dw3+c30*c31Dw3); double t22Dw3 = 2.0*(c02*c02Dw3+c12*c12Dw3+c22*c22Dw3+c32*c32Dw3); double t01Dw3 = t10Dw3; double t02Dw3 = t20Dw3; double t12Dw3 = t21Dw3; double t00Dx0 = 2.0 * (c00*c00Dx0+c10*c10Dx0+c20*c20Dx0+c30*c30Dx0); double t10Dx0 = (c00*c01Dx0+c01*c00Dx0) + (c10*c11Dx0+c11*c10Dx0) + (c20*c21Dx0+c21*c20Dx0) + (c30*c31Dx0+c31*c30Dx0); double t11Dx0 = 2.0 * (c01*c01Dx0+c11*c11Dx0+c21*c21Dx0+c31*c31Dx0); double t20Dx0 = (c00*c02Dx0+c02*c00Dx0) + (c10*c12Dx0+c12*c10Dx0) + (c20*c22Dx0+c22*c20Dx0) + (c30*c32Dx0+c32*c30Dx0); double t21Dx0 = (c01*c00Dx0+c00*c01Dx0) + (c11*c10Dx0+c10*c11Dx0) + (c21*c20Dx0+c20*c21Dx0) + (c31*c30Dx0+c30*c31Dx0); double t22Dx0 = 2.0*(c02*c02Dx0+c12*c12Dx0+c22*c22Dx0+c32*c32Dx0); double t01Dx0 = t10Dx0; double t02Dx0 = t20Dx0; double t12Dx0 = t21Dx0; double t00Dx1 = 2.0 * (c00*c00Dx1+c10*c10Dx1+c20*c20Dx1+c30*c30Dx1); double t10Dx1 = (c00*c01Dx1+c01*c00Dx1) + (c10*c11Dx1+c11*c10Dx1) + (c20*c21Dx1+c21*c20Dx1) + (c30*c31Dx1+c31*c30Dx1); double t11Dx1 = 2.0 * (c01*c01Dx1+c11*c11Dx1+c21*c21Dx1+c31*c31Dx1); double t20Dx1 = (c00*c02Dx1+c02*c00Dx1) + (c10*c12Dx1+c12*c10Dx1) + (c20*c22Dx1+c22*c20Dx1) + (c30*c32Dx1+c32*c30Dx1); double t21Dx1 = (c01*c00Dx1+c00*c01Dx1) + (c11*c10Dx1+c10*c11Dx1) + (c21*c20Dx1+c20*c21Dx1) + (c31*c30Dx1+c30*c31Dx1); double t22Dx1 = 2.0*(c02*c02Dx1+c12*c12Dx1+c22*c22Dx1+c32*c32Dx1); double t01Dx1 = t10Dx1; double t02Dx1 = t20Dx1; double t12Dx1 = t21Dx1; double t00Dx2 = 2.0 * (c00*c00Dx2+c10*c10Dx2+c20*c20Dx2+c30*c30Dx2); double t10Dx2 = (c00*c01Dx2+c01*c00Dx2) + (c10*c11Dx2+c11*c10Dx2) + (c20*c21Dx2+c21*c20Dx2) + (c30*c31Dx2+c31*c30Dx2); double t11Dx2 = 2.0 * (c01*c01Dx2+c11*c11Dx2+c21*c21Dx2+c31*c31Dx2); double t20Dx2 = (c00*c02Dx2+c02*c00Dx2) + (c10*c12Dx2+c12*c10Dx2) + (c20*c22Dx2+c22*c20Dx2) + (c30*c32Dx2+c32*c30Dx2); double t21Dx2 = (c01*c00Dx2+c00*c01Dx2) + (c11*c10Dx2+c10*c11Dx2) + (c21*c20Dx2+c20*c21Dx2) + (c31*c30Dx2+c30*c31Dx2); double t22Dx2 = 2.0*(c02*c02Dx2+c12*c12Dx2+c22*c22Dx2+c32*c32Dx2); double t01Dx2 = t10Dx2; double t02Dx2 = t20Dx2; double t12Dx2 = t21Dx2; double t00Dx3 = 2.0 * (c00*c00Dx3+c10*c10Dx3+c20*c20Dx3+c30*c30Dx3); double t10Dx3 = (c00*c01Dx3+c01*c00Dx3) + (c10*c11Dx3+c11*c10Dx3) + (c20*c21Dx3+c21*c20Dx3) + (c30*c31Dx3+c31*c30Dx3); double t11Dx3 = 2.0 * (c01*c01Dx3+c11*c11Dx3+c21*c21Dx3+c31*c31Dx3); double t20Dx3 = (c00*c02Dx3+c02*c00Dx3) + (c10*c12Dx3+c12*c10Dx3) + (c20*c22Dx3+c22*c20Dx3) + (c30*c32Dx3+c32*c30Dx3); double t21Dx3 = (c01*c00Dx3+c00*c01Dx3) + (c11*c10Dx3+c10*c11Dx3) + (c21*c20Dx3+c20*c21Dx3) + (c31*c30Dx3+c30*c31Dx3); double t22Dx3 = 2.0*(c02*c02Dx3+c12*c12Dx3+c22*c22Dx3+c32*c32Dx3); double t01Dx3 = t10Dx3; double t02Dx3 = t20Dx3; double t12Dx3 = t21Dx3; double T11Du0 = (t11*t22Du0+t22*t11Du0) - (t12*t12Du0+t12*t12Du0); double T22Du0 = (t00*t22Du0+t22*t00Du0) - (t02*t02Du0+t02*t02Du0); double T33Du0 = (t00*t11Du0+t11*t00Du0) - (t01*t01Du0+t01*t01Du0); double T13Du0 = (t01*t12Du0+t12*t01Du0) - (t11*t02Du0+t02*t11Du0); double T12Du0 = (t01*t22Du0+t22*t01Du0) - (t12*t02Du0+t02*t12Du0); double T11Du1 = (t11*t22Du1+t22*t11Du1) - (t12*t12Du1+t12*t12Du1); double T22Du1 = (t00*t22Du1+t22*t00Du1) - (t02*t02Du1+t02*t02Du1); double T33Du1 = (t00*t11Du1+t11*t00Du1) - (t01*t01Du1+t01*t01Du1); double T13Du1 = (t01*t12Du1+t12*t01Du1) - (t11*t02Du1+t02*t11Du1); double T12Du1 = (t01*t22Du1+t22*t01Du1) - (t12*t02Du1+t02*t12Du1); double T11Du2 = (t11*t22Du2+t22*t11Du2) - (t12*t12Du2+t12*t12Du2); double T22Du2 = (t00*t22Du2+t22*t00Du2) - (t02*t02Du2+t02*t02Du2); double T33Du2 = (t00*t11Du2+t11*t00Du2) - (t01*t01Du2+t01*t01Du2); double T13Du2 = (t01*t12Du2+t12*t01Du2) - (t11*t02Du2+t02*t11Du2); double T12Du2 = (t01*t22Du2+t22*t01Du2) - (t12*t02Du2+t02*t12Du2); double T11Du3 = (t11*t22Du3+t22*t11Du3) - (t12*t12Du3+t12*t12Du3); double T22Du3 = (t00*t22Du3+t22*t00Du3) - (t02*t02Du3+t02*t02Du3); double T33Du3 = (t00*t11Du3+t11*t00Du3) - (t01*t01Du3+t01*t01Du3); double T13Du3 = (t01*t12Du3+t12*t01Du3) - (t11*t02Du3+t02*t11Du3); double T12Du3 = (t01*t22Du3+t22*t01Du3) - (t12*t02Du3+t02*t12Du3); double T11Dv0 = (t11*t22Dv0+t22*t11Dv0) - (t12*t12Dv0+t12*t12Dv0); double T22Dv0 = (t00*t22Dv0+t22*t00Dv0) - (t02*t02Dv0+t02*t02Dv0); double T33Dv0 = (t00*t11Dv0+t11*t00Dv0) - (t01*t01Dv0+t01*t01Dv0); double T13Dv0 = (t01*t12Dv0+t12*t01Dv0) - (t11*t02Dv0+t02*t11Dv0); double T12Dv0 = (t01*t22Dv0+t22*t01Dv0) - (t12*t02Dv0+t02*t12Dv0); double T11Dv1 = (t11*t22Dv1+t22*t11Dv1) - (t12*t12Dv1+t12*t12Dv1); double T22Dv1 = (t00*t22Dv1+t22*t00Dv1) - (t02*t02Dv1+t02*t02Dv1); double T33Dv1 = (t00*t11Dv1+t11*t00Dv1) - (t01*t01Dv1+t01*t01Dv1); double T13Dv1 = (t01*t12Dv1+t12*t01Dv1) - (t11*t02Dv1+t02*t11Dv1); double T12Dv1 = (t01*t22Dv1+t22*t01Dv1) - (t12*t02Dv1+t02*t12Dv1); double T11Dv2 = (t11*t22Dv2+t22*t11Dv2) - (t12*t12Dv2+t12*t12Dv2); double T22Dv2 = (t00*t22Dv2+t22*t00Dv2) - (t02*t02Dv2+t02*t02Dv2); double T33Dv2 = (t00*t11Dv2+t11*t00Dv2) - (t01*t01Dv2+t01*t01Dv2); double T13Dv2 = (t01*t12Dv2+t12*t01Dv2) - (t11*t02Dv2+t02*t11Dv2); double T12Dv2 = (t01*t22Dv2+t22*t01Dv2) - (t12*t02Dv2+t02*t12Dv2); double T11Dv3 = (t11*t22Dv3+t22*t11Dv3) - (t12*t12Dv3+t12*t12Dv3); double T22Dv3 = (t00*t22Dv3+t22*t00Dv3) - (t02*t02Dv3+t02*t02Dv3); double T33Dv3 = (t00*t11Dv3+t11*t00Dv3) - (t01*t01Dv3+t01*t01Dv3); double T13Dv3 = (t01*t12Dv3+t12*t01Dv3) - (t11*t02Dv3+t02*t11Dv3); double T12Dv3 = (t01*t22Dv3+t22*t01Dv3) - (t12*t02Dv3+t02*t12Dv3); double T11Dw0 = (t11*t22Dw0+t22*t11Dw0) - (t12*t12Dw0+t12*t12Dw0); double T22Dw0 = (t00*t22Dw0+t22*t00Dw0) - (t02*t02Dw0+t02*t02Dw0); double T33Dw0 = (t00*t11Dw0+t11*t00Dw0) - (t01*t01Dw0+t01*t01Dw0); double T13Dw0 = (t01*t12Dw0+t12*t01Dw0) - (t11*t02Dw0+t02*t11Dw0); double T12Dw0 = (t01*t22Dw0+t22*t01Dw0) - (t12*t02Dw0+t02*t12Dw0); double T11Dw1 = (t11*t22Dw1+t22*t11Dw1) - (t12*t12Dw1+t12*t12Dw1); double T22Dw1 = (t00*t22Dw1+t22*t00Dw1) - (t02*t02Dw1+t02*t02Dw1); double T33Dw1 = (t00*t11Dw1+t11*t00Dw1) - (t01*t01Dw1+t01*t01Dw1); double T13Dw1 = (t01*t12Dw1+t12*t01Dw1) - (t11*t02Dw1+t02*t11Dw1); double T12Dw1 = (t01*t22Dw1+t22*t01Dw1) - (t12*t02Dw1+t02*t12Dw1); double T11Dw2 = (t11*t22Dw2+t22*t11Dw2) - (t12*t12Dw2+t12*t12Dw2); double T22Dw2 = (t00*t22Dw2+t22*t00Dw2) - (t02*t02Dw2+t02*t02Dw2); double T33Dw2 = (t00*t11Dw2+t11*t00Dw2) - (t01*t01Dw2+t01*t01Dw2); double T13Dw2 = (t01*t12Dw2+t12*t01Dw2) - (t11*t02Dw2+t02*t11Dw2); double T12Dw2 = (t01*t22Dw2+t22*t01Dw2) - (t12*t02Dw2+t02*t12Dw2); double T11Dw3 = (t11*t22Dw3+t22*t11Dw3) - (t12*t12Dw3+t12*t12Dw3); double T22Dw3 = (t00*t22Dw3+t22*t00Dw3) - (t02*t02Dw3+t02*t02Dw3); double T33Dw3 = (t00*t11Dw3+t11*t00Dw3) - (t01*t01Dw3+t01*t01Dw3); double T13Dw3 = (t01*t12Dw3+t12*t01Dw3) - (t11*t02Dw3+t02*t11Dw3); double T12Dw3 = (t01*t22Dw3+t22*t01Dw3) - (t12*t02Dw3+t02*t12Dw3); double T11Dx0 = (t11*t22Dx0+t22*t11Dx0) - (t12*t12Dx0+t12*t12Dx0); double T22Dx0 = (t00*t22Dx0+t22*t00Dx0) - (t02*t02Dx0+t02*t02Dx0); double T33Dx0 = (t00*t11Dx0+t11*t00Dx0) - (t01*t01Dx0+t01*t01Dx0); double T13Dx0 = (t01*t12Dx0+t12*t01Dx0) - (t11*t02Dx0+t02*t11Dx0); double T12Dx0 = (t01*t22Dx0+t22*t01Dx0) - (t12*t02Dx0+t02*t12Dx0); double T11Dx1 = (t11*t22Dx1+t22*t11Dx1) - (t12*t12Dx1+t12*t12Dx1); double T22Dx1 = (t00*t22Dx1+t22*t00Dx1) - (t02*t02Dx1+t02*t02Dx1); double T33Dx1 = (t00*t11Dx1+t11*t00Dx1) - (t01*t01Dx1+t01*t01Dx1); double T13Dx1 = (t01*t12Dx1+t12*t01Dx1) - (t11*t02Dx1+t02*t11Dx1); double T12Dx1 = (t01*t22Dx1+t22*t01Dx1) - (t12*t02Dx1+t02*t12Dx1); double T11Dx2 = (t11*t22Dx2+t22*t11Dx2) - (t12*t12Dx2+t12*t12Dx2); double T22Dx2 = (t00*t22Dx2+t22*t00Dx2) - (t02*t02Dx2+t02*t02Dx2); double T33Dx2 = (t00*t11Dx2+t11*t00Dx2) - (t01*t01Dx2+t01*t01Dx2); double T13Dx2 = (t01*t12Dx2+t12*t01Dx2) - (t11*t02Dx2+t02*t11Dx2); double T12Dx2 = (t01*t22Dx2+t22*t01Dx2) - (t12*t02Dx2+t02*t12Dx2); double T11Dx3 = (t11*t22Dx3+t22*t11Dx3) - (t12*t12Dx3+t12*t12Dx3); double T22Dx3 = (t00*t22Dx3+t22*t00Dx3) - (t02*t02Dx3+t02*t02Dx3); double T33Dx3 = (t00*t11Dx3+t11*t00Dx3) - (t01*t01Dx3+t01*t01Dx3); double T13Dx3 = (t01*t12Dx3+t12*t01Dx3) - (t11*t02Dx3+t02*t11Dx3); double T12Dx3 = (t01*t22Dx3+t22*t01Dx3) - (t12*t02Dx3+t02*t12Dx3); double delDu0 = (t00*T11Du0+t00Du0*T11) - (t01*T12Du0+t01Du0*T12) + (t02*T13Du0+t02Du0*T13); double gamDu0 = t00Du0 + t11Du0 + t22Du0; double sigDu0 = T11Du0 + T22Du0 + T33Du0; double delDu1 = (t00*T11Du1+t00Du1*T11) - (t01*T12Du1+t01Du1*T12) + (t02*T13Du1+t02Du1*T13); double gamDu1 = t00Du1 + t11Du1 + t22Du1; double sigDu1 = T11Du1 + T22Du1 + T33Du1; double delDu2 = (t00*T11Du2+t00Du2*T11) - (t01*T12Du2+t01Du2*T12) + (t02*T13Du2+t02Du2*T13); double gamDu2 = t00Du2 + t11Du2 + t22Du2; double sigDu2 = T11Du2 + T22Du2 + T33Du2; double delDu3 = (t00*T11Du3+t00Du3*T11) - (t01*T12Du3+t01Du3*T12) + (t02*T13Du3+t02Du3*T13); double gamDu3 = t00Du3 + t11Du3 + t22Du3; double sigDu3 = T11Du3 + T22Du3 + T33Du3; double delDv0 = (t00*T11Dv0+t00Dv0*T11) - (t01*T12Dv0+t01Dv0*T12) + (t02*T13Dv0+t02Dv0*T13); double gamDv0 = t00Dv0 + t11Dv0 + t22Dv0; double sigDv0 = T11Dv0 + T22Dv0 + T33Dv0; double delDv1 = (t00*T11Dv1+t00Dv1*T11) - (t01*T12Dv1+t01Dv1*T12) + (t02*T13Dv1+t02Dv1*T13); double gamDv1 = t00Dv1 + t11Dv1 + t22Dv1; double sigDv1 = T11Dv1 + T22Dv1 + T33Dv1; double delDv2 = (t00*T11Dv2+t00Dv2*T11) - (t01*T12Dv2+t01Dv2*T12) + (t02*T13Dv2+t02Dv2*T13); double gamDv2 = t00Dv2 + t11Dv2 + t22Dv2; double sigDv2 = T11Dv2 + T22Dv2 + T33Dv2; double delDv3 = (t00*T11Dv3+t00Dv3*T11) - (t01*T12Dv3+t01Dv3*T12) + (t02*T13Dv3+t02Dv3*T13); double gamDv3 = t00Dv3 + t11Dv3 + t22Dv3; double sigDv3 = T11Dv3 + T22Dv3 + T33Dv3; double delDw0 = (t00*T11Dw0+t00Dw0*T11) - (t01*T12Dw0+t01Dw0*T12) + (t02*T13Dw0+t02Dw0*T13); double gamDw0 = t00Dw0 + t11Dw0 + t22Dw0; double sigDw0 = T11Dw0 + T22Dw0 + T33Dw0; double delDw1 = (t00*T11Dw1+t00Dw1*T11) - (t01*T12Dw1+t01Dw1*T12) + (t02*T13Dw1+t02Dw1*T13); double gamDw1 = t00Dw1 + t11Dw1 + t22Dw1; double sigDw1 = T11Dw1 + T22Dw1 + T33Dw1; double delDw2 = (t00*T11Dw2+t00Dw2*T11) - (t01*T12Dw2+t01Dw2*T12) + (t02*T13Dw2+t02Dw2*T13); double gamDw2 = t00Dw2 + t11Dw2 + t22Dw2; double sigDw2 = T11Dw2 + T22Dw2 + T33Dw2; double delDw3 = (t00*T11Dw3+t00Dw3*T11) - (t01*T12Dw3+t01Dw3*T12) + (t02*T13Dw3+t02Dw3*T13); double gamDw3 = t00Dw3 + t11Dw3 + t22Dw3; double sigDw3 = T11Dw3 + T22Dw3 + T33Dw3; double delDx0 = (t00*T11Dx0+t00Dx0*T11) - (t01*T12Dx0+t01Dx0*T12) + (t02*T13Dx0+t02Dx0*T13); double gamDx0 = t00Dx0 + t11Dx0 + t22Dx0; double sigDx0 = T11Dx0 + T22Dx0 + T33Dx0; double delDx1 = (t00*T11Dx1+t00Dx1*T11) - (t01*T12Dx1+t01Dx1*T12) + (t02*T13Dx1+t02Dx1*T13); double gamDx1 = t00Dx1 + t11Dx1 + t22Dx1; double sigDx1 = T11Dx1 + T22Dx1 + T33Dx1; double delDx2 = (t00*T11Dx2+t00Dx2*T11) - (t01*T12Dx2+t01Dx2*T12) + (t02*T13Dx2+t02Dx2*T13); double gamDx2 = t00Dx2 + t11Dx2 + t22Dx2; double sigDx2 = T11Dx2 + T22Dx2 + T33Dx2; double delDx3 = (t00*T11Dx3+t00Dx3*T11) - (t01*T12Dx3+t01Dx3*T12) + (t02*T13Dx3+t02Dx3*T13); double gamDx3 = t00Dx3 + t11Dx3 + t22Dx3; double sigDx3 = T11Dx3 + T22Dx3 + T33Dx3; double te1Du0 = alpha/32.0*(2.0*Delta*delDu0-2.0*Del_3*delDu0); double te2Du0 = beta/6.0*(2.0*Sigma*gamDu0-3.0 *sigDu0); double te1Dv0 = alpha/32.0*(2.0*Delta*delDv0-2.0*Del_3*delDv0); double te2Dv0 = beta/6.0*(2.0*Sigma*gamDv0-3.0 *sigDv0); double te1Dw0 = alpha/32.0*(2.0*Delta*delDw0-2.0*Del_3*delDw0); double te2Dw0 = beta/6.0*(2.0*Sigma*gamDw0-3.0 *sigDw0); double te1Dx0 = alpha/32.0*(2.0*Delta*delDx0-2.0*Del_3*delDx0); double te2Dx0 = beta/6.0*(2.0*Sigma*gamDx0-3.0 *sigDx0); double te1Du1 = alpha/32.0*(2.0*Delta*delDu1-2.0*Del_3*delDu1); double te2Du1 = beta/6.0*(2.0*Sigma*gamDu1-3.0 *sigDu1); double te1Dv1 = alpha/32.0*(2.0*Delta*delDv1-2.0*Del_3*delDv1); double te2Dv1 = beta/6.0*(2.0*Sigma*gamDv1-3.0 *sigDv1); double te1Dw1 = alpha/32.0*(2.0*Delta*delDw1-2.0*Del_3*delDw1); double te2Dw1 = beta/6.0*(2.0*Sigma*gamDw1-3.0 *sigDw1); double te1Dx1 = alpha/32.0*(2.0*Delta*delDx1-2.0*Del_3*delDx1); double te2Dx1 = beta/6.0*(2.0*Sigma*gamDx1-3.0 *sigDx1); double te1Du2 = alpha/32.0*(2.0*Delta*delDu2-2.0*Del_3*delDu2); double te2Du2 = beta/6.0*(2.0*Sigma*gamDu2-3.0 *sigDu2); double te1Dv2 = alpha/32.0*(2.0*Delta*delDv2-2.0*Del_3*delDv2); double te2Dv2 = beta/6.0*(2.0*Sigma*gamDv2-3.0 *sigDv2); double te1Dw2 = alpha/32.0*(2.0*Delta*delDw2-2.0*Del_3*delDw2); double te2Dw2 = beta/6.0*(2.0*Sigma*gamDw2-3.0 *sigDw2); double te1Dx2 = alpha/32.0*(2.0*Delta*delDx2-2.0*Del_3*delDx2); double te2Dx2 = beta/6.0*(2.0*Sigma*gamDx2-3.0 *sigDx2); double te1Du3 = alpha/32.0*(2.0*Delta*delDu3-2.0*Del_3*delDu3); double te2Du3 = beta/6.0*(2.0*Sigma*gamDu3-3.0 *sigDu3); double te1Dv3 = alpha/32.0*(2.0*Delta*delDv3-2.0*Del_3*delDv3); double te2Dv3 = beta/6.0*(2.0*Sigma*gamDv3-3.0 *sigDv3); double te1Dw3 = alpha/32.0*(2.0*Delta*delDw3-2.0*Del_3*delDw3); double te2Dw3 = beta/6.0*(2.0*Sigma*gamDw3-3.0 *sigDw3); double te1Dx3 = alpha/32.0*(2.0*Delta*delDx3-2.0*Del_3*delDx3); double te2Dx3 = beta/6.0*(2.0*Sigma*gamDx3-3.0 *sigDx3); double epDu0 = volA * ( te1Du0 + te2Du0); double epDv0 = volA * ( te1Dv0 + te2Dv0); double epDw0 = volA * ( te1Dw0 + te2Dw0); double epDx0 = volA * ( te1Dx0 + te2Dx0); double epDu1 = volA * ( te1Du1 + te2Du1); double epDv1 = volA * ( te1Dv1 + te2Dv1); double epDw1 = volA * ( te1Dw1 + te2Dw1); double epDx1 = volA * ( te1Dx1 + te2Dx1); double epDu2 = volA * ( te1Du2 + te2Du2); double epDv2 = volA * ( te1Dv2 + te2Dv2); double epDw2 = volA * ( te1Dw2 + te2Dw2); double epDx2 = volA * ( te1Dx2 + te2Dx2); double epDu3 = volA * ( te1Du3 + te2Du3); double epDv3 = volA * ( te1Dv3 + te2Dv3); double epDw3 = volA * ( te1Dw3 + te2Dw3); double epDx3 = volA * ( te1Dx3 + te2Dx3); double eDu0 = eDep * epDu0; double eDv0 = eDep * epDv0; double eDw0 = eDep * epDw0; double eDx0 = eDep * epDx0; double eDu1 = eDep * epDu1; double eDv1 = eDep * epDv1; double eDw1 = eDep * epDw1; double eDx1 = eDep * epDx1; double eDu2 = eDep * epDu2; double eDv2 = eDep * epDv2; double eDw2 = eDep * epDw2; double eDx2 = eDep * epDx2; double eDu3 = eDep * epDu3; double eDv3 = eDep * epDv3; double eDw3 = eDep * epDw3; double eDx3 = eDep * epDx3; { if (vVar.el[iu]) { eDc[iu,0] = eDc[iu,0] + eDu0; eDc[iu,1] = eDc[iu,1] + eDu1; eDc[iu,2] = eDc[iu,2] + eDu2; eDc[iu,3] = eDc[iu,3] + eDu3; } if (vVar.el[iv]) { eDc[iv,0] = eDc[iv,0] + eDv0; eDc[iv,1] = eDc[iv,1] + eDv1; eDc[iv,2] = eDc[iv,2] + eDv2; eDc[iv,3] = eDc[iv,3] + eDv3; } if (vVar.el[iw]) { eDc[iw,0] = eDc[iw,0] + eDw0; eDc[iw,1] = eDc[iw,1] + eDw1; eDc[iw,2] = eDc[iw,2] + eDw2; eDc[iw,3] = eDc[iw,3] + eDw3; } if (vVar.el[ix]) { eDc[ix,0] = eDc[ix,0] + eDx0; eDc[ix,1] = eDc[ix,1] + eDx1; eDc[ix,2] = eDc[ix,2] + eDx2; eDc[ix,3] = eDc[ix,3] + eDx3; } } } Distribute_eDep; { /* Clear potential energy accumulators: */ for (j = 0; j < NP; j++){ ep[j] = 0.0; } /* Accumulate potential energy of each tetrahedron */ for (j = 0; j < NP; j++) { if (polyRelevant[j]) { Accum_ep(j, ep[j]); } } /* Compute energy "e" from cells determinants, and the gradient "eDdet": */ e = 0.0; for (p = 0; p < NP; p++) { if (polyRelevant[p]) { Accum_e_from_ep(ep[p], e, eDep[p]); } else { eDep[p] = 0.0 } } /* Now distribute "eDdet" over "eDc": */ for (i = 0; i < NV; i++){ eDc[i] = (r4_t){0.0, ..}; } if (grad) { for (j = 0; j < NP; j++) { if (polyRelevant[j]) { Distribute_eDep(j,eDep[j]) } } } } } } /* END Eval */ char * SELF_Name(SELF_T *erg)== { return "Elasticity(alpha="&Fmt.LongReal(erg->alpha,Fmt.Style.Fix,prec=2) \ " beta="&Fmt.LongReal(erg->beta,Fmt.Style.Fix,prec=2)&")"; } /* END Name */ { } ElasticityEnergy.