/* See {CmpVolEnergy.h} */ #include #include #include #include #include #include CONST Zero == 0.0; TYPE bool_vec_t == ARRAY OF bool_t; double_vec_t == double_vec_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 */ vp: REF double_vec_t; /* (Work) Volume of each cell */ eDvp: REF double_vec_t; /* (Work) Derivate of energy rel. to volume cells */ 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/((double)top->cell.nel); erg->top = top; erg->vVar = bool_vec_new(top->node.nel); erg->polyRelevant = bool_vec_new(top->cell.nel); /* Allocate volume tables: */ erg->vp = double_vec_new(top->cell.nel); erg->eDvp = double_vec_new(top->cell.nel); /* 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 */ void SELF_DefVar(SELF_T *erg, bool_vec_t *variable) { /* Decide which cells are relevant to volume compression energy. A cell is relevant iff it exists, has exactly four walls being least one existing triangular wall having at least one variable corner. */ int NV = erg->top->node.nel; ??? NP = erg->top->cell.nel; bool_vec_t vVar = erg->vVar^; ??? cell = 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++) { ??? p = cell[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))); bool_t fvar = (f1->exists) || (f2->exists) || (f3->exists) || (f4->exists); { assert(f1 != f2) && (f2 != f3) && (f3 != f4); if ((p->exists) && (fvar) && (vvar)) { /*assert(u->exists) && (v->exists) && (w->exists) && (x->exists);*/ 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; bool_vec_t vVar = erg->vVar^; ??? polyRelevant = erg->polyRelevant^; ??? cell = erg->top->cell; ??? vp = erg->vp^; ??? eDvp = erg->eDvp^; ??? K = erg->K; ??? V = ((double)erg->volume); { void Accum_v(*u, v, w, x: r4_t; VAR vo: double) /* Compute the volume of the tetrahedron "u v w x". Also stores in {n} the cross product "(v - u) X (w - u) X (x - u)". */ { ??? uv = r4_Sub(v, u); ??? uw == r4_Sub(w, u); ??? ux == r4_Sub(x, u); ??? n == r4_cross(uv, uw, ux); { vo = vo + 1.0/6.0 * r4_Norm(n); } } Accum_v; void Accum_e_from_v(v: double; double e; VAR VAR etDv: double) /* Adds to "e" the energy term corresponding to a cell with volume "v". Also, if "grad" is true, stores in "etDv" the derivative of that term. */ { assert(v!=0.0); if (v <= 0.0) { etDv = 0.0 } else { double r = v/V; double s = 1.0/r; double r2 = r*r; double s2 = s*s; double h = r2 + s2 - 2.0; double v2 = v*v; double v3 = v2*v; double V2 == V*V; { e = e + K * h; if (grad) { etDv = 2.0 * K * ( (v/V2) - (V2/v3)); } else { etDv = 0.0; } } } } Accum_e_from_v; void Distribute_eDv(iu, iv, iw, ix: uint; eDvo: double) /* Accumulates in "eDc" the gradient of "e" relative to the corners of the tetrahedron "iu iv iw ix", given the derivative "eDv" of "e" relative to the tetrahedron's volume "v". */ { r4_t u = c[iu], v = c[iv], w = c[iw], x = c[ix]; double u1 = u[0], u2 = u[1], u3 = u[2], u4 = u[3]; double v1 = v[0], v2 = v[1], v3 = v[2], v4 = v[3]; double w1 = w[0], w2 = w[1], w3 = w[2], w4 = w[3]; double x1 = x[0], x2 = x[1], x3 = x[2], x4 = x[3]; double d1 = +u1 * (v2 * (w3-x3) - v3 * (w2-x2) + w2*x3-w3*x2) -u2 * (v1 * (w3-x3) - v3 * (w1-x1) + w1*x3-w3*x1) +u3 * (v1 * (w2-x2) - v2 * (w1-x1) + w1*x2-x1*w2) -(v1*(w2*x3-x2*w3)-v2*(w1*x3-x1*w3)+v3*(w1*x2-x1*w2)), d2 = +u1 * (v2 * (w4-x4) - v4 * (w2-x2) + w2*x4-x2*w4) -u2 * (v1 * (w4-x4) - v4 * (w1-x1) + w1*x4-x1*w4) +u4 * (v1 * (w2-x2) - v2 * (w1-x1) + w1*x2-x1*w2) -(v1*(w2*x4-x2*w4)-v2*(w1*x4-x1*w4)+v4*(w1*x2-x1*w2)), d3 = +u2 * (v3 * (w4-x4) - v4 * (w3-x3) + w3*x4-x3*w4) -u3 * (v2 * (w4-x4) - v4 * (w2-x2) + w2*x4-x2*w4) +u4 * (v2 * (w3-x3) - v3 * (w2-x2) + w2*x3-x2*w3) -(v2*(w3*x4-x3*w4)-v3*(w2*x4-x2*w4)+v4*(w2*x3-x2*w3)), d4 = +u1 * (v3 * (w4-x4) - v4 * (w3-x3) + w3*x4-x3*w4) -u3 * (v1 * (w4-x4) - v4 * (w1-x1) + w1*x4-x1*w4) +u4 * (v1 * (w3-x3) - v3 * (w1-x1) + w1*x3-x1*w3) -(v1*(w3*x4-x3*w4)-v3*(w1*x4-x1*w4)+v4*(w1*x3-x1*w3)); r4_t d = (r4_t){d1,d2,d3,d4}; double sum2 = d1*d1 + d2*d2 + d3*d3 + d4*d4; double term = (1.0/6.0) * (1.0/sqrt(sum2)); double d1Du1 = v2 * (w3-x3) - v3 * (w2-x2) + w2*x3-w3*x2; double d2Du1 = v2 * (w4-x4) - v4 * (w2-x2) + w2*x4-x2*w4; double d3Du1 = Zero; double d4Du1 = v3 * (w4-x4) - v4 * (w3-x3) + w3*x4-x3*w4; r4_t dDu1 = (r4_t){d1Du1, d2Du1, d3Du1, d4Du1}, double vDu1 = term * r4_dot(d,dDu1); double eDu1 = eDvo * vDu1; double d1Du2 = - (v1 * (w3-x3) - v3 * (w1-x1) + w1*x3-w3*x1); double d2Du2 = - (v1 * (w4-x4) - v4 * (w1-x1) + w1*x4-x1*w4); double d3Du2 = v3 * (w4-x4) - v4 * (w3-x3) + w3*x4-x3*w4; double d4Du2 = Zero; r4_t dDu2 = (r4_t){d1Du2, d2Du2, d3Du2, d4Du2}, double vDu2 = term * r4_dot(d,dDu2); double eDu2 = eDvo * vDu2; double d1Du3 = v1 * (w2-x2) - v2 * (w1-x1) + w1*x2-x1*w2; double d2Du3 = Zero; double d3Du3 = - (v2 * (w4-x4) - v4 * (w2-x2) + w2*x4-x2*w4); double d4Du3 = - (v1 * (w4-x4) - v4 * (w1-x1) + w1*x4-x1*w4); r4_t dDu3 = (r4_t){d1Du3, d2Du3, d3Du3, d4Du3}, double vDu3 = term * r4_dot(d,dDu3); double eDu3 = eDvo * vDu3; double d1Du4 = Zero; double d2Du4 = v1 * (w2-x2) - v2 * (w1-x1) + w1*x2-x1*w2; double d3Du4 = v2 * (w3-x3) - v3 * (w2-x2) + w2*x3-x2*w3; double d4Du4 = v1 * (w3-x3) - v3 * (w1-x1) + w1*x3-x1*w3; r4_t dDu4 = (r4_t){d1Du4, d2Du4, d3Du4, d4Du4}, double vDu4 = term * r4_dot(d,dDu4); double eDu4 = eDvo * vDu4; double d1Dv1 = u2 * (x3-w3) + u3 * (w2-x2) - w2*x3+w3*x2; double d2Dv1 = u2 * (x4-w4) + u4 * (w2-x2) - w2*x4+x2*w4; double d3Dv1 = Zero; double d4Dv1 = u3 * (x4-w4) + u4 * (w3-x3) - w3*x4+x3*w4; r4_t dDv1 = (r4_t){d1Dv1, d2Dv1, d3Dv1, d4Dv1}, double vDv1 = term * r4_dot(d,dDv1); double eDv1 = eDvo * vDv1; double d1Dv2 = u1 * (w3-x3) + u3 * (x1-w1) + w1*x3-w3*x1; double d2Dv2 = u1 * (w4-x4) + u4 * (x1-w1) + w1*x4-x1*w4; double d3Dv2 = u3 * (x4-w4) + u4 * (w3-x3) - w3*x4+x3*w4; double d4Dv2 = Zero; r4_t dDv2 = (r4_t){d1Dv2, d2Dv2, d3Dv2, d4Dv2}, double vDv2 = term * r4_dot(d,dDv2); double eDv2 = eDvo * vDv2; double d1Dv3 = u1 * (x2-w2) + u2 * (w1-x1) - w1*x2+x1*w2; double d2Dv3 = Zero; double d3Dv3 = u2 * (w4-x4) + u4 * (x2-w2) + w2*x4-x2*w4; double d4Dv3 = u1 * (w4-x4) + u4 * (x1-w1) + w1*x4-x1*w4; r4_t dDv3 = (r4_t){d1Dv3, d2Dv3, d3Dv3, d4Dv3}, double vDv3 = term * r4_dot(d,dDv3); double eDv3 = eDvo * vDv3; double d1Dv4 = Zero; double d2Dv4 = u1 * (x2-w2) + u2 * (w1-x1) - w1*x2+x1*w2; double d3Dv4 = u2 * (x3-w3) + u3 * (w2-x2) - w2*x3+x2*w3; double d4Dv4 = u1 * (x3-w3) + u3 * (w1-x1) - w1*x3+x1*w3; r4_t dDv4 = (r4_t){d1Dv4, d2Dv4, d3Dv4, d4Dv4}, double vDv4 = term * r4_dot(d,dDv4); double eDv4 = eDvo * vDv4; double d1Dw1 = u2 * (v3-x3) + u3 * (x2-v2) + v2*x3-v3*x2; double d2Dw1 = u2 * (v4-x4) + u4 * (x2-v2) + v2*x4-x2*v4; double d3Dw1 = Zero; double d4Dw1 = u3 * (v4-x4) + u4 * (x3-v3) + v3*x4-x3*v4; double dDw1 = (r4_t){d1Dw1, d2Dw1, d3Dw1, d4Dw1}; double vDw1 = term * r4_dot(d,dDw1); double eDw1 = eDvo * vDw1; double d1Dw2 = u1 * (x3-v3) + u3 * (v1-x1) - v1*x3+v3*x1; double d2Dw2 = u1 * (x4-v4) + u4 * (v1-x1) - v1*x4+x1*v4; double d3Dw2 = u3 * (v4-x4) + u4 * (x3-v3) + v3*x4-x3*v4; double d4Dw2 = Zero; r4_t dDw2 = (r4_t){d1Dw2, d2Dw2, d3Dw2, d4Dw2}, double vDw2 = term * r4_dot(d,dDw2); double eDw2 = eDvo * vDw2; double d1Dw3 = u1 * (v2-x2) + u2 * (x1-v1) + v1*x2-x1*v2; double d2Dw3 = Zero; double d3Dw3 = u2 * (x4-v4) + u4 * (v2-x2) - v2*x4+x2*v4; double d4Dw3 = u1 * (x4-v4) + u4 * (v1-x1) - v1*x4+x1*v4; r4_t dDw3 = (r4_t){d1Dw3, d2Dw3, d3Dw3, d4Dw3}, double vDw3 = term * r4_dot(d,dDw3); double eDw3 = eDvo * vDw3; double d1Dw4 = Zero; double d2Dw4 = u1 * (v2-x2) + u2 * (x1-v1) + v1*x2-x1*v2; double d3Dw4 = u2 * (v3-x3) + u3 * (x2-v2) + v2*x3-x2*v3; double d4Dw4 = u1 * (v3-x3) + u3 * (x1-v1) + v1*x3-x1*v3; r4_t dDw4 = (r4_t){d1Dw4, d2Dw4, d3Dw4, d4Dw4}, double vDw4 = term * r4_dot(d,dDw4); double eDw4 = eDvo * vDw4; double d1Dx1 = u2 * (w3-v3) + u3 * (v2-w2) - v2*w3+v3*w2; double d2Dx1 = u2 * (w4-v4) + u4 * (v2-w2) - v2*w4+w2*v4; double d3Dx1 = Zero; double d4Dx1 = u3 * (w4-v4) + u4 * (v3-w3) - v3*w4+w3*v4; r4_t dDx1 = (r4_t){d1Dx1, d2Dx1, d3Dx1, d4Dx1}, double vDx1 = term * r4_dot(d,dDx1); double eDx1 = eDvo * vDx1; double d1Dx2 = u1 * (v3-w3) + u3 * (w1-v1) + v1*w3-v3*w1; double d2Dx2 = u1 * (v4-w4) + u4 * (w1-v1) + v1*w4-w1*v4; double d3Dx2 = u3 * (w4-v4) + u4 * (v3-w3) + v4*w3-v3*w4; double d4Dx2 = Zero; r4_t dDx2 = (r4_t){d1Dx2, d2Dx2, d3Dx2, d4Dx2}, double vDx2 = term * r4_dot(d,dDx2); double eDx2 = eDvo * vDx2; double d1Dx3 = u1 * (w2-v2) + u2 * (v1-w1) - v1*w2+w1*v2; double d2Dx3 = Zero; double d3Dx3 = u2 * (v4-w4) + u4 * (w2-v2) + v2*w4-w2*v4; double d4Dx3 = u1 * (v4-w4) + u4 * (w1-v1) + v1*w4-w1*v4; r4_t dDx3 = (r4_t){d1Dx3, d2Dx3, d3Dx3, d4Dx3}, double vDx3 = term * r4_dot(d,dDx3); double eDx3 = eDvo * vDx3; double d1Dx4 = Zero; double d2Dx4 = u1 * (w2-v2) + u2 * (v1-w1) - v1*w2+w1*v2; double d3Dx4 = u2 * (w3-v3) + u3 * (v2-w2) - v2*w3+w2*v3; double d4Dx4 = u1 * (w3-v3) + u3 * (v1-w1) - v1*w3+w1*v3; r4_t dDx4 = (r4_t){d1Dx4, d2Dx4, d3Dx4, d4Dx4}, double vDx4 = term * r4_dot(d,dDx4); eDx4 = eDvo * vDx4 ){ if (vVar.el[iu]) { eDc[iu,0] = eDc[iu,0] + eDu1; eDc[iu,1] = eDc[iu,1] + eDu2; eDc[iu,2] = eDc[iu,2] + eDu3; eDc[iu,3] = eDc[iu,3] + eDu4; } if (vVar.el[iv]) { eDc[iv,0] = eDc[iv,0] + eDv1; eDc[iv,1] = eDc[iv,1] + eDv2; eDc[iv,2] = eDc[iv,2] + eDv3; eDc[iv,3] = eDc[iv,3] + eDv4; } if (vVar.el[iw]) { eDc[iw,0] = eDc[iw,0] + eDw1; eDc[iw,1] = eDc[iw,1] + eDw2; eDc[iw,2] = eDc[iw,2] + eDw3; eDc[iw,3] = eDc[iw,3] + eDw4; } if (vVar.el[ix]) { eDc[ix,0] = eDc[ix,0] + eDx1; eDc[ix,1] = eDc[ix,1] + eDx2; eDc[ix,2] = eDc[ix,2] + eDx3; eDc[ix,3] = eDc[ix,3] + eDx4; } } } Distribute_eDv; { /* Clear volume accumulators: */ for (j = 0; j < NP; j++){ vp[j] = 0.0; } /* Enumerate cells and accumulate cell volumes: */ for (j = 0; j < NP; j++) { if (polyRelevant[j]) { ??? t = Srot(cell[j]); Place_t a = Tors(t); Place_t t1 = Clock(PrevE(t)); Place_t a1 = Tors(t1); Node_t un = OrgV(a)->num; Node_t vn = OrgV(NextE(a))->num; Node_t wn = OrgV(PrevE(a))->num; Node_t xn = OrgV(PrevE(a1))->num; { Accum_v(c[un], c[vn], c[wn], c[xn], vp[j]) } } } /* Compute energy "e" from cells volumes, and the gradient "eDvp": */ e = 0.0; for (p = 0; p < NP; p++) { if (polyRelevant[p]) { Accum_e_from_v(vp[p], e, eDvp[p]) } else { eDvp[p] = 0.0 } } /* Now distribute "eDvp" 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]) { ??? t = Srot(cell[j]); Place_t a = Tors(t); Place_t t1 = Clock(PrevE(t)); Place_t a1 = Tors(t1); Node_t un = OrgV(a)->num; Node_t vn = OrgV(NextE(a))->num; Node_t wn = OrgV(PrevE(a))->num; Node_t xn = OrgV(PrevE(a1))->num; { Distribute_eDv(un, vn, wn, xn, eDvp[j]) } } } } } } } /* END Eval */ char * SELF_Name(SELF_T *erg)== { return "Compr(iVol = " & Fmt.Real(erg->volume, Fmt.Style.Fix, prec = 3) & ")" } /* END Name */