/* See {UneqVolEnergy.h} */ #include #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 */ bool_vec_t wallRelevant; /* TRUE if wall is relevant */ vr: REF double_vec_t; /* (Work) Volume ratio of two cells incident to wall relevant */ eDvr: REF double_vec_t; /* (Work) Derivate of energy rel. to volume ratio */ 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->wall.nel); erg->top = top; erg->vVar = bool_vec_new(top->node.nel); erg->wallRelevant = bool_vec_new(top->wall.nel); /* Allocate volume ratio tables: */ erg->vr = double_vec_new(top->wall.nel); erg->eDvr = double_vec_new(top->wall.nel); /* Just in case the client forgets to call "defVar": */ for (i = 0; i < top->node.nel; i++){ erg->vVar^[i] = FALSE; } for (i = 0; i < top->wall.nel; i++){ erg->wallRelevant^[i] = FALSE; } } /* END DefTop */ void SELF_DefVar(SELF_T *erg, bool_vec_t *variable) VAR p1,p2: Cell; { /* Decide which wall are relevant to unequal volume energy. A wall is relevant iff it exists, has exactly two cell incident to it being least one existing cell and the wall have at least one variable corner. */ int NV = erg->top->node.nel; int NF = erg->top->wall.nel; bool_vec_t vVar = erg->vVar^; ??? wall = erg->top->wall^; ??? wallRelevant = erg->wallRelevant^; { assert((variable.nel) == NV); vVar = variable; /* Find the relevant walls: */ for (i = 0; i < NF; i++){ wallRelevant.el[i] = FALSE; } for (i = 0; i < NF; i++) { ??? f = wall[i]; ??? a = f.pa^; { p1 = PnegP(a); p2 = PposP(a); if ((p1!=NULL) && (p2!=NULL)) { ??? pvar = p1->exists) || (p2->exists; Node_t u = OrgV(a); Node_t v = OrgV(NextE(a)); Node_t w = OrgV(PrevE(a)); ??? vvar = vVar.el[u->num]) || (vVar.el[v->num]) || (vVar.el[w->num]; { if ((f->exists) && (pvar) && (vvar)) { assert(u->exists) && (v->exists) && (w->exists); wallRelevant.el[i] = TRUE; } } } } } } } /* END DefVar */ void SELF_Eval( SELF_T erg; Coords_t *c; double *e; bool_t grad; Gradient *eDc; ) { ??? K = erg->K; int NV = erg->top->node.nel; int NF = erg->top->wall.nel; ??? wall = erg->top->wall^; bool_vec_t vVar = erg->vVar^; ??? wallRelevant = erg->wallRelevant^; ??? vr = erg->vr^; ??? eDvr = erg->eDvr^; { void Accum_vr(*u, v, w, x, y: r4_t; VAR vro: double) /* Compute the volume ratio of tetrahedrons "u w v x" (PposP) and "u w y v" (PnegP). "n1" is the cross product "(w - u)X(v - u)X(x - u)" and "n2" is the cross product (w - u)X(y - u)X(v - u)". The coordinates "u v w" defined the common wall of two tetrahedrons where "x" is the node extremus of cell PposP and "y" is the node extremus of cell PnegP. Note that the rigth-hand rule is valid in R^{4} for determining in which of the two possible direc- tions "a x b x c" points. */ { ??? uv = r4_Sub(v, u); ??? uw = r4_Sub(w, u); ??? ux = r4_Sub(x, u); ??? uy = r4_Sub(y, u); ??? n1 = r4_cross(uw, uv, ux); ??? n2 = r4_cross(uw, uy, uv); ??? v1 = 1.0/6.0 * r4_Norm(n1); ??? v2 = 1.0/6.0 * r4_Norm(n2); { vro = vro + v1/v2; } } Accum_vr; void Accum_e_from_vr(vr: double; double e; VAR VAR etDvr: double) /* Adds to "e" the energy term corresponding to a wall with volume ratio "vr" of cells incidents to it. Also, if "grad" is true, stores in "etDvr" the derivative of that term. */ { assert(vr!=0.0); if (vr <= 0.0) { etDvr = 0.0 } else { ??? d = vr - 1.0; ??? d2 = d * d; { e = e + K * d2; if (grad) { etDvr = 2.0 * K * d; } else { etDvr = 0.0; } } } } Accum_e_from_vr; void Distribute_eDvr(iu, iv, iw, ix, iy: uint; eDvr: double) /* Accumulates in "eDc" the gradient of "e" relative to the corners of the two tetrahedrons "iu iv iw ix", and "iu iv iw iy" given the derivative "eDv" of "e" relative to volume ratio of tetrahedrons "vr". */ void CalcPartialDerivatives( uint *u,v,w,x; double sum2; VAR VAR voDu,voDv,voDw,voDx: r4_t; ) double *term; { ??? u = c[u], v = c[v], w = c[w], x = c[x]; ??? u1 = u[0], u2 = u[1], u3 = u[2], u4 = u[3]; ??? v1 = v[0], v2 = v[1], v3 = v[2], v4 = v[3]; ??? w1 = w[0], w2 = w[1], w3 = w[2], w4 = w[3]; ??? 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)); double 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)); double 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)); double 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)); double d = (r4_t){d1,d2,d3,d4}; 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; double dDu1 = (r4_t){d1Du1, d2Du1, d3Du1, d4Du1}; 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; double dDu2 = (r4_t){d1Du2, d2Du2, d3Du2, d4Du2}; 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); double dDu3 = (r4_t){d1Du3, d2Du3, d3Du3, d4Du3}; 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; double dDu4 = (r4_t){d1Du4, d2Du4, d3Du4, d4Du4}; 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; double dDv1 = (r4_t){d1Dv1, d2Dv1, d3Dv1, d4Dv1}; 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; double dDv2 = (r4_t){d1Dv2, d2Dv2, d3Dv2, d4Dv2}; 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; double dDv3 = (r4_t){d1Dv3, d2Dv3, d3Dv3, d4Dv3}; 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; double dDv4 = (r4_t){d1Dv4, d2Dv4, d3Dv4, d4Dv4}; 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 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; double dDw2 = (r4_t){d1Dw2, d2Dw2, d3Dw2, d4Dw2}; 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; double dDw3 = (r4_t){d1Dw3, d2Dw3, d3Dw3, d4Dw3}; 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; double dDw4 = (r4_t){d1Dw4, d2Dw4, d3Dw4, d4Dw4}; 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; double dDx1 = (r4_t){d1Dx1, d2Dx1, d3Dx1, d4Dx1}; 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; double dDx2 = (r4_t){d1Dx2, d2Dx2, d3Dx2, d4Dx2}; 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; double dDx3 = (r4_t){d1Dx3, d2Dx3, d3Dx3, d4Dx3}; 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; double dDx4 = (r4_t){d1Dx4, d2Dx4, d3Dx4, d4Dx4}; { sum2 = d1*d1 + d2*d2 + d3*d3 + d4*d4; term = 1.0/sqrt(sum2); voDu[0] = term * r4_dot(d,dDu1); voDu[1] = term * r4_dot(d,dDu2); voDu[2] = term * r4_dot(d,dDu3); voDu[3] = term * r4_dot(d,dDu4); voDv[0] = term * r4_dot(d,dDv1); voDv[1] = term * r4_dot(d,dDv2); voDv[2] = term * r4_dot(d,dDv3); voDv[3] = term * r4_dot(d,dDv4); voDw[0] = term * r4_dot(d,dDw1); voDw[1] = term * r4_dot(d,dDw2); voDw[2] = term * r4_dot(d,dDw3); voDw[3] = term * r4_dot(d,dDw4); voDx[0] = term * r4_dot(d,dDx1); voDx[1] = term * r4_dot(d,dDx2); voDx[2] = term * r4_dot(d,dDx3); voDx[3] = term * r4_dot(d,dDx4); voDu = (r4_t){voDu[0], voDu[1], voDu[2], voDu[3]}; voDv = (r4_t){voDv[0], voDv[1], voDv[2], voDv[3]}; voDw = (r4_t){voDw[0], voDw[1], voDw[2], voDw[3]}; voDx = (r4_t){voDx[0], voDx[1], voDx[2], voDx[3]}; } } CalcPartialDerivatives; VAR vpDu,vpDv,vpDw,vpDx,vnDu,vnDv,vnDw,vnDy : r4_t; sum1, sum2 : double; { CalcPartialDerivatives(iu,iv,iw,ix,sum1,vpDu,vpDv,vpDw,vpDx); CalcPartialDerivatives(iu,iv,iw,iy,sum2,vnDu,vnDv,vnDw,vnDy); ??? vrDu1 = (sqrt(sum2) * vpDu[0] - sqrt(sum1) * vnDu[0])/sum2; ??? vrDu2 = (sqrt(sum2) * vpDu[1] - sqrt(sum1) * vnDu[1])/sum2; ??? vrDu3 = (sqrt(sum2) * vpDu[2] - sqrt(sum1) * vnDu[2])/sum2; ??? vrDu4 = (sqrt(sum2) * vpDu[3] - sqrt(sum1) * vnDu[3])/sum2; double vrDv1 = (sqrt(sum2) * vpDv[0] - sqrt(sum1) * vnDv[0])/sum2; double vrDv2 = (sqrt(sum2) * vpDv[1] - sqrt(sum1) * vnDv[1])/sum2; double vrDv3 = (sqrt(sum2) * vpDv[2] - sqrt(sum1) * vnDv[2])/sum2; double vrDv4 = (sqrt(sum2) * vpDv[3] - sqrt(sum1) * vnDv[3])/sum2; double vrDw1 = (sqrt(sum2) * vpDw[0] - sqrt(sum1) * vnDw[0])/sum2; double vrDw2 = (sqrt(sum2) * vpDw[1] - sqrt(sum1) * vnDw[1])/sum2; double vrDw3 = (sqrt(sum2) * vpDw[2] - sqrt(sum1) * vnDw[2])/sum2; double vrDw4 = (sqrt(sum2) * vpDw[3] - sqrt(sum1) * vnDw[3])/sum2; double vrDx1 = (sqrt(sum2) * vpDx[0])/sum2; double vrDx2 = (sqrt(sum2) * vpDx[1])/sum2; double vrDx3 = (sqrt(sum2) * vpDx[2])/sum2; double vrDx4 = (sqrt(sum2) * vpDx[3])/sum2; double vrDy1 = -(sqrt(sum1) * vnDy[0])/sum2; double vrDy2 = -(sqrt(sum1) * vnDy[1])/sum2; double vrDy3 = -(sqrt(sum1) * vnDy[2])/sum2; double vrDy4 = -(sqrt(sum1) * vnDy[3])/sum2; double eDu1 = eDvr * vrDu1; double eDu2 = eDvr * vrDu2; double eDu3 = eDvr * vrDu3; double eDu4 = eDvr * vrDu4; double eDv1 = eDvr * vrDv1; double eDv2 = eDvr * vrDv2; double eDv3 = eDvr * vrDv3; double eDv4 = eDvr * vrDv4; double eDw1 = eDvr * vrDw1; double eDw2 = eDvr * vrDw2; double eDw3 = eDvr * vrDw3; double eDw4 = eDvr * vrDw4; double eDx1 = eDvr * vrDx1; double eDx2 = eDvr * vrDx2; double eDx3 = eDvr * vrDx3; double eDx4 = eDvr * vrDx4; double eDy1 = eDvr * vrDy1; double eDy2 = eDvr * vrDy2; double eDy3 = eDvr * vrDy3; double eDy4 = eDvr * vrDy4 { 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; } if (vVar.el[iy]) { eDc[iy,0] = eDc[iy,0] + eDy1; eDc[iy,1] = eDc[iy,1] + eDy2; eDc[iy,2] = eDc[iy,2] + eDy3; eDc[iy,3] = eDc[iy,3] + eDy4; } } } Distribute_eDvr; { /* Clear volume ratio accumulators: */ for (j = 0; j < NF; j++){ vr[j] = 0.0; } /* Enumerate walls and accumulate volume ratio of cells incidents. */ for (j = 0; j < NF; j++) { if (wallRelevant.el[j]) { ??? f = wall[j]; ??? a = f.pa^; 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(NextF(a)))->num; Node_t yn = OrgV(PrevE(PrevF(a)))->num; { Accum_vr(c[un], c[vn], c[wn], c[xn], c[yn], vr[j]) } } } /* Compute energy "e" from volumes ratios, and the gradient "eDvr": */ e = 0.0; for (j = 0; j < NF; j++) { if (wallRelevant.el[j]) { Accum_e_from_vr(vr[j], e, eDvr[j]) } else { eDvr[j] = 0.0 } } /* Now distribute "eDvr" over "eDc": */ for (i = 0; i < NV; i++){ eDc[i] = (r4_t){0.0, ..}; } if (grad) { for (j = 0; j < NF; j++) { if (wallRelevant.el[j]) { ??? f = wall[j]; ??? a = f.pa^; 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(NextF(a)))->num; Node_t yn = OrgV(PrevE(PrevF(a)))->num; { Distribute_eDvr(un, vn, wn, xn, yn, eDvr[j]) } } } } } } } /* END Eval */ PROCEDURE Name(<*UNUSED); erg: T): char *== { return "Unequal()"; } /* END Name */ { } UneqVolEnergy.