#define PROG_NAME "BuildRefinement" #define PROG_DESC "???" #define PROG_VERS "1.0" #define BuildRefinement_C_COPYRIGHT \ "" #define PROG_INFO \ "" \ " " /* This module builds a refinement of a triangulation $T$ exchanging each single tetrahedron of $T$ by a $K$-refined tetrahedron as created by the "Refine" procedure on the library "libm3subi". */ #include #define _GNU_SOURCE #include #include #include #include #define _GNU_SOURCE #include #include #include TYPE typedef struct Options_t { char *inFileTp; /* Input file name (minus extensions) */ char *inFileSt; /* Input file name (minus extensions) */ char *outFile; /* Output file name prefix */ uint order; /* order of the refinement */ bool_t detail; /* print some detail information */ bool_t fixed; /* retain the previus geometry */ bool_t net; /* drawing each existing wall as a net with small spheres and thin cylinders */ } double Corner = ARRAY [0..23] OF Place_t; double Report = RECORD bool_t bo; INTEGER tu[1+1]; } @PLACE == Place_vec_t; double CENTER = Node_vec_t; PROCEDURE MyAssert(*p: r4_t) { assert(!( (p[0] == 0.0) AND (p[1] == 0.0) AND (p[2] == 0.0) AND (p[3] == 0.0) ) ); } /* END MyAssert */ bool_t AreGlued(Place_t @p, Place_t b) /* This procedure return TRUE iff the corners (@places) "a" e "b" are glued in the original triangulation $T$; otherwise return FALSE. */ bool_t CellsInfo() /* First step: We need to compare the cells information. */ { if (((PposP(a) == PnegP(b))) && ((PposP(b) == PnegP(a)))) { return TRUE; } else { return FALSE } } CellsInfo; bool_t NodesInfo() /* Second step: We need to compare the nodes information. */ { if (((OrgV(a) == OrgV(b)) AND (OrgV(NextE(a)) == OrgV(NextE(b))) AND (OrgV(PrevE(a)) == OrgV(PrevE(b))))){ return TRUE; } else { return FALSE } } NodesInfo; { if ((CellsInfo())){ if ((NodesInfo())) { assert(a == Spin(b)) && (b == Spin(a)); return TRUE; } } return FALSE; } /* END AreGlued */ Node MidOcta(Place_t @p) /* Return a node that is the medial node of an triangulated octahedron. This node has a "VP" label as type. */ { Node_t v = OrgV(PrevE(PrevF(PrevF(a)))); { v.label = "VP"; return v; } } /* END MidOcta */ Node InterNode(Place_t @p) /* Return a node that is a internal to the macro-tetrahedron. This node is located on a @{edge->?} and wall and present only in order greater than four. This node has a "VP" label as type. */ { Node_t v = OrgV(PrevE(PrevF(a))); { v.label = "FE"; return v; } } /* END InterNode */ Report AreGluedCorners(*ca,cb: Corner) /* This procedure return TRUE iff the tetrahedra with corners "ca" and "cb" must be glued. */ Report *re; { for (i = 0; i <= (23); i++){ for (j = 0; j <= (23); j++) { if ((AreGlued(ca.el[i],cb.el[j]))) { re.bo = TRUE; re.tu[0] = i; re.tu[1] = j; } } } return re; } /* END AreGluedCorners */ Options_t *GetOptions(int argc, char **argv); int main(int argc, char **argv) VAR ElemTableRec_t newtop; vme : REF ARRAY OF CENTER; /* Array of mid node of octahedra */ vin : REF ARRAY OF CENTER; /* An array of internal node */ nvv,nve,nvf,nvp,nfe: uint = 0; /* numbers nodes of each type */ nov,noe,nof,nop : uint = 0; { Options_t *o = GetOptions(argc, argv); char *topo_cmt = jsprintf("Created by %s on %s", PROG_NAME, Today()); /* Random_t coins = MakeRandomSource(4615); */ ??? tc = Triangulation.ReadToMa(o->inFileTp); ??? top = tc.top; ??? co = NEW(REF ARRAY OF Corner, top->cell.nel); ??? nco = NEW(REF ARRAY OF Refine.Corner, top->cell.nel); ??? rc = Triangulation.ReadState(o->inFileSt), ??? c = rc^ { /* Assertion */ if (top->der!=3) { fprintf(stderr,"\nThis topology isn't a triangulation\n"); Process.Exit¦(1); } fprintf(stderr,"\n" & Fmt.Int(o->order) & "-Refining from: " & o->inFileTp & ".tp\n"); /* compute Corners table */ for (i = 0; i < top->cell.nel; i++) { ??? da = Srot(top.cell[i]); Place_t pa = Tors(da); { co.el[i] = BuildCornerTable(pa); nco.el[i] = SubdivideTetrahedron(o->order, pa, o->net, top); } if (o->detail) { fprintf(stderr, "tetrahedron " & Fmt.Int(i) & ":\n"); PrintCornerTable(co.el[i]); } } /* save memory for the medial and internal nodes .*/ if (o->order == 1) { vme = NEW(REF ARRAY OF CENTER, top->cell.nel, 0); } else if (o->order == 2){ vme = NEW(REF ARRAY OF CENTER, top->cell.nel, 1); } else if (o->order == 3){ vme = NEW(REF ARRAY OF CENTER, top->cell.nel, 4); } else if (o->order == 4){ vme = NEW(REF ARRAY OF CENTER, top->cell.nel, 10); vin = NEW(REF ARRAY OF CENTER, top->cell.nel, 1); } else if (o->order == 5){ vme = NEW(REF ARRAY OF CENTER, top->cell.nel, 20); vin = NEW(REF ARRAY OF CENTER, top->cell.nel, 4); } for (i = 0; i < top->cell.nel; i++) { if (o->order == 2) { vme[i][0] = MidOcta(nco.el[i].front[0]); } else if (o->order == 3){ vme[i][ 0] = MidOcta(nco.el[i].front[ 0]); vme[i][ 1] = MidOcta(nco.el[i].front[ 1]); vme[i][ 2] = MidOcta(nco.el[i].front[ 3]); vme[i][ 3] = MidOcta(nco.el[i].back [ 0]); } else if (o->order == 4){ vme[i][ 0] = MidOcta(nco.el[i].front[ 0]); vme[i][ 1] = MidOcta(nco.el[i].front[ 1]); vme[i][ 2] = MidOcta(nco.el[i].front[ 3]); vme[i][ 3] = MidOcta(nco.el[i].front[ 4]); vme[i][ 4] = MidOcta(nco.el[i].front[ 6]); vme[i][ 5] = MidOcta(nco.el[i].front[ 8]); vme[i][ 6] = MidOcta(nco.el[i].back [ 0]); vme[i][ 7] = MidOcta(nco.el[i].right[ 6]); vme[i][ 8] = MidOcta(nco.el[i].right[13]); vme[i][ 9] = MidOcta(nco.el[i].left [13]); /* */ vin[i][ 0] = InterNode(nco.el[i].front[6]); } else if (o->order == 5){ vme[i][ 0] = MidOcta(nco.el[i].front[ 0]); vme[i][ 1] = MidOcta(nco.el[i].front[ 1]); vme[i][ 2] = MidOcta(nco.el[i].front[ 3]); vme[i][ 3] = MidOcta(nco.el[i].front[ 4]); vme[i][ 4] = MidOcta(nco.el[i].front[ 6]); vme[i][ 5] = MidOcta(nco.el[i].front[ 8]); vme[i][ 6] = MidOcta(nco.el[i].front[ 9]); vme[i][ 7] = MidOcta(nco.el[i].front[11]); vme[i][ 8] = MidOcta(nco.el[i].front[13]); vme[i][ 9] = MidOcta(nco.el[i].front[15]); vme[i][10] = MidOcta(nco.el[i].back [ 0]); vme[i][11] = MidOcta(nco.el[i].right[11]); vme[i][12] = MidOcta(nco.el[i].right[20]); vme[i][13] = MidOcta(nco.el[i].right[ 6]); vme[i][14] = MidOcta(nco.el[i].right[13]); vme[i][15] = MidOcta(nco.el[i].right[22]); vme[i][16] = MidOcta(nco.el[i].left [20]); vme[i][17] = MidOcta(nco.el[i].left [13]); vme[i][18] = MidOcta(nco.el[i].left [22]); vme[i][19] = MidOcta(nco.el[i].back [ 6]); /* */ vin[i][ 0] = InterNode(nco.el[i].front[ 6]); vin[i][ 1] = InterNode(nco.el[i].front[11]); vin[i][ 2] = InterNode(nco.el[i].front[13]); vin[i][ 3] = InterNode(nco.el[i].right[11]); } } /* compute that corners that must be glued */ for (k = 0; k < top->wall.nel; k++) { ??? f = top->wall[k].pa; ??? i = Triangulation.PposP(f)->num; ??? j = Triangulation.PnegP(f)->num; { if (((Triangulation.PposP(f)!=NULL)) && ((Triangulation.PnegP(f)!=NULL))) { ??? re = AreGluedCorners(co.el[i],co.el[j]); { if (re.bo) { if (o->detail) { fprintf(stderr, "New Corners that must be glue :" \ " NC[" & Pad(Int(i),3) &"," & Pad(Int(re.tu[0]),3)&"] and" \ " NC[" & Pad(Int(j),3) &"," & Pad(Int(re.tu[1]),3)&"]\n"); } ??? ar = Refine.Place_tsOnFrontier(nco.el[i],re.tu[0], o->order); double br = Refine.Place_tsOnFrontier(nco.el[j], re.tu[1], o->order) { GlueMacroTetrahedra(ar,br); } } } } } } for (i = 0; i < top->node.nel; i++) { ??? v = NARROW(top->node[i], Node); { double vl = v.label ){ if (( Text.Equal(vl,"VV"))){ nov++; } else if (0 == strcmp(vl,"VE"))){ noe++; } else if (0 == strcmp(vl,"VF"))){ nof++; } else if (0 == strcmp(vl,"VP"))){ nop++; } } } newtop = MakeElemTable(nco.el[0].right[0],1); for (i = 0; i < newtop.node.nel; i++) { ??? v = NARROW(newtop.node[i], Node); { double vl = v.label ){ if (( Text.Equal(vl,"VV"))){ nvv++; } else if (0 == strcmp(vl,"VE"))){ nve++; } else if (0 == strcmp(vl,"VF"))){ nvf++; } else if (0 == strcmp(vl,"VP"))){ nvp++; } else if (0 == strcmp(vl,"FE"))){ nfe++; } } } assert(newtop.node.nel == nvv + nve + nvf + nvp + nfe); assert(nvv == nov); if (o->order >= 2) { assert(nve == top->NE * (o->order-1) + noe); } if ( o->order == 3){ assert(nvf == top->wall.nel * 1 + nof); } else if (o->order == 4){ assert(nvf == top->wall.nel * 3 + nof); } else if (o->order == 5){ assert(nvf == top->wall.nel * 6 + nof); } /* useful assertion */ ??? o1 = o->order; ??? o2 = o1*o1; ??? o3 = o2*o1; ??? cte = (5*o3-2*o1) DIV 3; { assert(newtop.cell.nel == top->cell.nel * cte); if ( o->order == 1){ assert(nvp == top->cell.nel * 0 + nop); } else if (o->order == 2){ assert(nvp == top->cell.nel * 1 + nop); } else if (o->order == 3){ assert(nvp == top->cell.nel * 4 + nop); } else if (o->order == 4 ){ assert(nvp == top->cell.nel * 10 + nop); assert(nfe == top->cell.nel * 1); } else if (o->order == 5 ){ assert(nvp == top->cell.nel * 20 + nop); assert(nfe == top->cell.nel * 4); } } ??? com = Fmt.Int(o->order) & "-Refined from: " & o->inFileTp & ".tp\n" & "Created by BuildRefinement: " & o->outFile & ".tp"; ??? ncord = SettingCoords(top,newtop,c,co,nco,o->order,vme,vin,o->net); double rcord = GenCoords(newtop)^; { WriteTopology (o->outFile, newtop, com); /*WriteTable (o->outFile, newtop, com);*/ WriteMaterials(o->outFile, newtop, com); if (o->fixed) { for (i = 0 TO LAST(ncord^)-1) { ncord[i] = r4_Scale(((double)1.0), ncord[i]); } WriteState(o->outFile&"-Fix",newtop,ncord^,com&": Fixed Geometry\n"); } else { WriteState(o->outFile, newtop, rcord, com & ": Random Geometry\n"); } } } } /* END DoIt */ Corner BuildCornerTable(Place_t @p) /* This procedure builds an array of 24 wedge @places, such all elements on the array have the same negative cell "PnegP". The @place "a" in the argument must be a @place that represent one tetrahe- dron of $T$ (in the primal space). */ Corner *c; { ??? a0 = a; Place_t a1 = NextE(a0); Place_t a2 = NextE(a1); Place_t b0 = Spin(PrevF(a0)); Place_t c0 = Spin(PrevF(a1)); Place_t d0 = Spin(PrevF(a2)); { /* the 24 elements on the array must have the same negative cell. */ c[0 ] = a0; c[1 ] = NextE(a0); c[2 ] = PrevE(a0); c[3 ] = b0; c[4 ] = NextE(b0); c[5 ] = PrevE(b0); c[6 ] = c0; c[7 ] = NextE(c0); c[8 ] = PrevE(c0); c[9 ] = d0; c[10] = NextE(d0); c[11] = PrevE(d0); assert(c[0] == a); /* basically the last 12 @places in the "corner" array are the "spined- clocked" version of the firsts 12 places. */ for (j = 1; j <= (12); j++) { c[j+11] = Clock(Spin(c[j-1])); } /* for assert that the 24 @places wedge have the same negative cell. */ for (m = 0; m <= (23); m++) { for (n = m+1; n <= (23); n++) { assert(PnegP(c[m]) == PnegP(c[n])); } } } return c; } /* END BuildCornerTable */ PROCEDURE PrintCornerTable(*ac : Corner) /* Prints detail information about the macro-tetrahedra that must be glued. */ { for (j = 0; j < 12; j++){ fprintf(stderr,Fmt.Pad(Fmt.Int(j),6) & " "); } fprintf(stderr,"\n--------------------------------------------------" \ "---------------------------------------\n"); for (j = 0; j < 12; j++) { Octf.PrintPlace(stderr, ac[j], 3, FALSE); } fprintf(stderr, "\ng\n"); for (j = 12; j <= (23); j++) { fprintf(stderr,Fmt.Pad(Fmt.Int(j),6) & " "); } fprintf(stderr,"\n--------------------------------------------------" \ "---------------------------------------\n"); for (j = 12; j <= (23); j++) { Octf.PrintPlace(stderr, ac[j], 3, FALSE); } fprintf(stderr, "\n\n"); } /* END PrintCornerTable */ Refine.Corner SubdivideTetrahedron( uint order; Place_t @p; bool_t net; ElemTableRec_t *top) /* This procedure rescues the propiertes of the nodes, @{edge->?}s and walls belonging to the original tetrahedron represented by the @place "a", and expands this propiertes to the macro-tetrahedron. */ void InheritNodes(v,V: Node) /* Inherits node propiertes.*/ { V->exists = v->exists; V.fixed = v.fixed; V.color = v.color; V.transp = v.transp; V.radius = v.radius; V.label = v.label; } InheritNodes; PROCEDURE Inherit@{Edge->?}s(e,E: @{Edge->?}) == /* Inherits @{edge->?} propiertes.*/ { ??? ed = top.edge[e->num]; { E->exists = ed->exists; E.color = ed.color; E.transp = ed.transp; E.radius = ed.radius; E->root = ed->root; } } Inherit@{Edge->?}s; void InheritWalls(f,F: Wall) /* Inherits wall propiertes.*/ { ??? fd = top.wall[f->num]; { F->exists = fd->exists; F.color = fd.color; F.transp = fd.transp; F->root = fd->root; } } InheritWalls; co : Refine.Corner; { if ( order == 1){ co = MakeTetra(1,net); } else if (order == 2){ co = MakeTetra(2,net); } else if (order == 3){ co = MakeTetra(3,net); } else if (order == 4){ co = MakeTetra(4,net); } else if (order == 5){ co = MakeTetra(5,net); } else { fprintf(stderr,"BuildRefinement: Order must be less equal to 5\n"); assert(order <= 5); } ??? o2 = order*order; ??? ve = Triangulation.TetraNegNodes(a); ??? ed = Triangulation.TetraEdges(a); ??? fa = Triangulation.TetraWalls(a); ??? v0 = ve[0]; ??? v1 = ve[1]; ??? v2 = ve[2]; ??? v3 = ve[3]; ??? e0 = ed[0]; ??? e1 = ed[1]; ??? e2 = ed[2]; ??? e3 = ed[3]; ??? e4 = ed[4]; ??? e5 = ed[5]; ??? f0 = fa[0]; ??? f1 = fa[1]; ??? f2 = fa[2]; ??? f3 = fa[3]; Node_t V0 = OrgV(PrevE(co.back [0])); ??? V1 = OrgV(PrevE(co.front[0])); ??? V2 = OrgV(NextE (co.front[o2-1])); ??? V3 = OrgV(PrevE(co.left [o2-1])); { /* nodes */ InheritNodes(v0,V0); InheritNodes(v1,V1); InheritNodes(v2,V2); InheritNodes(v3,V3); /* walls */ for (i = 0; i < o2; i++) { InheritWalls(f0, PWedge(co.right[i]).wall); InheritWalls(f1, PWedge(co.left [i]).wall); InheritWalls(f2, PWedge(co.front[i]).wall); InheritWalls(f3, PWedge(co.back [i]).wall); } if (net) { /* iff the wall not exists then non exists the elements on the net. */ ??? f0e = top.wall[f0->num]->exists; ??? f1e = top.wall[f1->num]->exists; ??? f2e = top.wall[f2->num]->exists; ??? f3e = top.wall[f3->num]->exists; { ModSetGrade(order,co,f0e,f1e,f2e,f3e); } } /* @{edge->?}s : e0 */ for (i = 1; i <= (order); i++) { Inherit@{Edge->?}s(e0, PWedge(co.right[(i-1)*(i-1)]).edge); } if (! top.edge[e0->num]->exists) { for (i = 1; i < order; i++) { OrgV(co.right[(i-1)*(i-1)])->exists = FALSE; } } ??? cor = top.edge[e0->num].color; { for (i = 1; i < order; i++) { OrgV(co.right[(i-1)*(i-1)]).color = cor; } } /* e1 */ for (i = 1; i <= (order); i++) { Inherit@{Edge->?}s(e1,PWedge( NextE(co.right[i*i-1])).edge); } if (! top.edge[e1->num]->exists) { for (i = 1; i < order; i++) { OrgV(PrevE(co.right[i*i-1]))->exists = FALSE; } } ??? cor = top.edge[e1->num].color; { for (i = 1; i < order; i++) { OrgV(PrevE(co.right[i*i-1])).color = cor; } } /* e3 */ for (i = 1; i <= (order); i++) { Inherit@{Edge->?}s(e3, PWedge(NextE(co.left[i*i-1])).edge); } if (! top.edge[e3->num]->exists) { for (i = 1; i < order; i++) { OrgV(PrevE(co.left[i*i-1]))->exists = FALSE; } } ??? cor = top.edge[e3->num].color; { for (i = 1; i < order; i++) { OrgV(PrevE(co.left[i*i-1])).color = cor; } } /* e4 */ for (i = 1; i <= (order); i++) { Inherit@{Edge->?}s(e4, PWedge(PrevE(co.back[(i-1)*(i-1)])).edge); } if (! top.edge[e4->num]->exists) { for (i = 1; i < order; i++) { OrgV(co.back[(i-1)*(i-1)])->exists = FALSE; } } ??? cor = top.edge[e4->num].color; { for (i = 1; i < order; i++) { OrgV(co.back[(i-1)*(i-1)]).color = cor; } } /* e2 */ for (i = 1; i <= (order); i++) { Inherit@{Edge->?}s(e2, PWedge(NextE(co.back[i*i-1])).edge); } if (! top.edge[e2->num]->exists) { for (i = 1; i < order; i++) { OrgV(NextE(co.back[i*i-1]))->exists = FALSE; } } ??? cor = top.edge[e2->num].color; { for (i = 1; i < order; i++) { OrgV(NextE(co.back[i*i-1])).color = cor; } } /* e5 */ INTEGER count = o2-1; VAR { Inherit@{Edge->?}s(e5, PWedge(co.front[count]).edge); for (j = 1; j < order; j++) { Inherit@{Edge->?}s(e5, PWedge(co.front[count-2]).edge); count = count-2; }; count = o2-1; if (! top.edge[e5->num]->exists) { for (j = 1; j < order; j++) { OrgV(co.front[count])->exists = FALSE; count = count-2; } }; count = o2-1; ??? cor = top.edge[e5->num].color; { for (j = 1; j < order; j++) { OrgV(co.front[count]).color = cor; count = count-2; } } } return co; } } /* END SubdivideTetrahedron */ PROCEDURE SettingCoords( ElemTableRec_t *oldtop; ElemTableRec_t *newtop; *oldco : Triangulation.Coords_t; *cor : REF ARRAY OF Corner; *ncor : REF ARRAY OF Refine.Corner; uint order; mocta : REF ARRAY OF CENTER; inter : REF ARRAY OF CENTER; bool_t net; ) : REF Triangulation.Coords_t == uint NewNumVer(Node_t v) /*Return the new number node of the node "v" on the "newtp" table.*/ VAR uint n = 0; { for (i = 0; i < newtp.NV; i++) { ??? w = newtp.node[i]; { if (w == v){ n = i; EXIT; } } } return n; } NewNumVer; bool_t PresentNode(Node_t v) /* Return TRUE iff the node "v" is present on the "newtp" table.*/ { for (i = 0; i < newtp.NV; i++) { ??? w = newtp.node[i]; { if (w == v){ return TRUE; }; } } return FALSE; } PresentNode; VAR newco : REF Triangulation.Coords_t; /* nodes */ v0,v1,v2,v3 : uint; /* @{edge->?}s */ e0 = NEW(REF uint_vec_t , order-1); e1 = NEW(REF uint_vec_t , order-1); e2 = NEW(REF uint_vec_t , order-1); e3 = NEW(REF uint_vec_t , order-1); e4 = NEW(REF uint_vec_t , order-1); e5 = NEW(REF uint_vec_t , order-1); /* walls */ f0 : REF uint_vec_t; f1 : REF uint_vec_t; f2 : REF uint_vec_t; f3 : REF uint_vec_t; { if (order == 3){ f0 = NEW(REF uint_vec_t , 1); f1 = NEW(REF uint_vec_t , 1); f2 = NEW(REF uint_vec_t , 1); f3 = NEW(REF uint_vec_t , 1); } else if (order == 4){ f0 = NEW(REF uint_vec_t , 3); f1 = NEW(REF uint_vec_t , 3); f2 = NEW(REF uint_vec_t , 3); f3 = NEW(REF uint_vec_t , 3); } else if (order == 5){ f0 = NEW(REF uint_vec_t , 6); f1 = NEW(REF uint_vec_t , 6); f2 = NEW(REF uint_vec_t , 6); f3 = NEW(REF uint_vec_t , 6); } newco = NEW(REF Triangulation.Coords_t, newtp.NV); /* original nodes */ for (i = 0; i < oldtp.NP; i++) { ??? o2 = order*order; Node_t V0 = OrgV(PrevE(ncor[i].back [0])); ??? V1 = OrgV(PrevE(ncor[i].front[0])); ??? V2 = OrgV(NextE (ncor[i].front[o2-1])); ??? V3 = OrgV(PrevE(ncor[i].left [o2-1])); { v0 = OrgV(cor[i,0])->num; v1 = OrgV(cor[i,1])->num; v2 = OrgV(cor[i,2])->num; v3 = OrgV(cor[i,5])->num; if ((PresentNode(V0))) { newco.el[NewNumVer(V0)] = oldco.el[v0]; } if ((PresentNode(V1))) { newco.el[NewNumVer(V1)] = oldco.el[v1]; } if ((PresentNode(V2))) { newco.el[NewNumVer(V2)] = oldco.el[v2]; } if ((PresentNode(V3))) { newco.el[NewNumVer(V3)] = oldco.el[v3]; } } /* nodes on the subdivided @{edge->?}s */ /* e0 : tips nodes V0 and v1 */ ??? dx = (oldco.el[v1][0] - oldco.el[v0][0])/((double)order); ??? dy = (oldco.el[v1][1] - oldco.el[v0][1])/((double)order); ??? dz = (oldco.el[v1][2] - oldco.el[v0][2])/((double)order); ??? dw = (oldco.el[v1][3] - oldco.el[v0][3])/((double)order); { for (j = 1; j < order; j++) { Node_t ve = OrgV(ncor[i].right[(j-1)*(j-1)]); { if ((PresentNode(ve))) { e0[j-1] = NewNumVer(ve); newco.el[e0[j-1]][0] = oldco.el[v1][0] - ((double)j) * dx; newco.el[e0[j-1]][1] = oldco.el[v1][1] - ((double)j) * dy; newco.el[e0[j-1]][2] = oldco.el[v1][2] - ((double)j) * dz; newco.el[e0[j-1]][3] = oldco.el[v1][3] - ((double)j) * dw; } } } } /* end e0 @{edge->?} */ /* e1 : tips nodes V1 and v2 */ ??? dx = (oldco.el[v1][0] - oldco.el[v2][0])/((double)order); ??? dy = (oldco.el[v1][1] - oldco.el[v2][1])/((double)order); ??? dz = (oldco.el[v1][2] - oldco.el[v2][2])/((double)order); ??? dw = (oldco.el[v1][3] - oldco.el[v2][3])/((double)order); { for (j = 1; j < order; j++) { Node_t ve = OrgV(PrevE(ncor[i].right[j*j-1])); { if ((PresentNode(ve))) { e1[j-1] = NewNumVer(ve); newco.el[e1[j-1]][0] = oldco.el[v1][0] - ((double)j) * dx; newco.el[e1[j-1]][1] = oldco.el[v1][1] - ((double)j) * dy; newco.el[e1[j-1]][2] = oldco.el[v1][2] - ((double)j) * dz; newco.el[e1[j-1]][3] = oldco.el[v1][3] - ((double)j) * dw; } } } } /* end e1 @{edge->?} */ /* e3 : tips nodes V1 and v3 */ ??? dx = (oldco.el[v1][0] - oldco.el[v3][0])/((double)order); ??? dy = (oldco.el[v1][1] - oldco.el[v3][1])/((double)order); ??? dz = (oldco.el[v1][2] - oldco.el[v3][2])/((double)order); ??? dw = (oldco.el[v1][3] - oldco.el[v3][3])/((double)order); { for (j = 1; j < order; j++) { Node_t ve = OrgV(PrevE(ncor[i].left[j*j-1])); { if ((PresentNode(ve))) { e3[j-1] = NewNumVer(ve); newco.el[e3[j-1]][0] = oldco.el[v1][0] - ((double)j) * dx; newco.el[e3[j-1]][1] = oldco.el[v1][1] - ((double)j) * dy; newco.el[e3[j-1]][2] = oldco.el[v1][2] - ((double)j) * dz; newco.el[e3[j-1]][3] = oldco.el[v1][3] - ((double)j) * dw; } } } } /* end e3 @{edge->?} */ /* e4 : tips nodes V0 and v3 */ ??? dx = (oldco.el[v0][0] - oldco.el[v3][0])/((double)order); ??? dy = (oldco.el[v0][1] - oldco.el[v3][1])/((double)order); ??? dz = (oldco.el[v0][2] - oldco.el[v3][2])/((double)order); ??? dw = (oldco.el[v0][3] - oldco.el[v3][3])/((double)order); { for (j = 1; j < order; j++) { Node_t ve = OrgV(ncor[i].back[(j-1)*(j-1)]); { if ((PresentNode(ve))) { e4[j-1] = NewNumVer(ve); newco.el[e4[j-1]][0] = oldco.el[v0][0] - ((double)j) * dx; newco.el[e4[j-1]][1] = oldco.el[v0][1] - ((double)j) * dy; newco.el[e4[j-1]][2] = oldco.el[v0][2] - ((double)j) * dz; newco.el[e4[j-1]][3] = oldco.el[v0][3] - ((double)j) * dw; } } } } /* end e4 @{edge->?} */ /* e2 : tips nodes V0 and v2 */ ??? dx = (oldco.el[v0][0] - oldco.el[v2][0])/((double)order); ??? dy = (oldco.el[v0][1] - oldco.el[v2][1])/((double)order); ??? dz = (oldco.el[v0][2] - oldco.el[v2][2])/((double)order); ??? dw = (oldco.el[v0][3] - oldco.el[v2][3])/((double)order); { for (j = 1; j < order; j++) { Node_t ve = OrgV(NextE(ncor[i].back[j*j-1])); { if ((PresentNode(ve))) { e2[j-1] = NewNumVer(ve); newco.el[e2[j-1]][0] = oldco.el[v0][0] - ((double)j) * dx; newco.el[e2[j-1]][1] = oldco.el[v0][1] - ((double)j) * dy; newco.el[e2[j-1]][2] = oldco.el[v0][2] - ((double)j) * dz; newco.el[e2[j-1]][3] = oldco.el[v0][3] - ((double)j) * dw; } } } } /* end e2 @{edge->?} */ /* e5 : tips nodes V2 and v3 */ ??? dx = (oldco.el[v2][0] - oldco.el[v3][0])/((double)order); ??? dy = (oldco.el[v2][1] - oldco.el[v3][1])/((double)order); ??? dz = (oldco.el[v2][2] - oldco.el[v3][2])/((double)order); ??? dw = (oldco.el[v2][3] - oldco.el[v3][3])/((double)order); { VAR uint count; { count = order*order-1; for (j = 1; j < order; j++) { Node_t ve = OrgV(ncor[i].front[count]); { count = count-2; if ((PresentNode(ve))) { e5[j-1] = NewNumVer(ve); newco.el[e5[j-1]][0] = oldco.el[v2][0] - ((double)j) * dx; newco.el[e5[j-1]][1] = oldco.el[v2][1] - ((double)j) * dy; newco.el[e5[j-1]][2] = oldco.el[v2][2] - ((double)j) * dz; newco.el[e5[j-1]][3] = oldco.el[v2][3] - ((double)j) * dw; } } } } } /* end e5 @{edge->?} */ /* nodes on the triangulated wall. */ if (order == 3) { Node_t F0 = OrgV(ncor[i].right[3]); ??? F1 = OrgV(ncor[i].left [3]); ??? F3 = OrgV(ncor[i].back [3]); ??? F2 = OrgV(ncor[i].front[3]); { ??? dx = (newco.el[e1[0]][0] + newco.el[e2[0]][0])/2.0; ??? dy = (newco.el[e1[0]][1] + newco.el[e2[0]][1])/2.0; ??? dz = (newco.el[e1[0]][2] + newco.el[e2[0]][2])/2.0; ??? dw = (newco.el[e1[0]][3] + newco.el[e2[0]][3])/2.0; { if ((PresentNode(F0))) { if (! net){ assert(! F0->exists); } f0[0] = NewNumVer(F0); newco.el[f0[0]][0] = dx; newco.el[f0[0]][1] = dy; newco.el[f0[0]][2] = dz; newco.el[f0[0]][3] = dw; } } ??? dx = (newco.el[e0[0]][0] + newco.el[e4[1]][0])/((double)order-1); ??? dy = (newco.el[e0[0]][1] + newco.el[e4[1]][1])/((double)order-1); ??? dz = (newco.el[e0[0]][2] + newco.el[e4[1]][2])/((double)order-1); ??? dw = (newco.el[e0[0]][3] + newco.el[e4[1]][3])/((double)order-1); { if ((PresentNode(F1))) { if (! net){ assert(! F1->exists); } f1[0] = NewNumVer(F1); newco.el[f1[0]][0] = dx; newco.el[f1[0]][1] = dy; newco.el[f1[0]][2] = dz; newco.el[f1[0]][3] = dw; } } ??? dx = (newco.el[e3[0]][0] + newco.el[e5[0]][0])/((double)order-1); ??? dy = (newco.el[e3[0]][1] + newco.el[e5[0]][1])/((double)order-1); ??? dz = (newco.el[e3[0]][2] + newco.el[e5[0]][2])/((double)order-1); ??? dw = (newco.el[e3[0]][3] + newco.el[e5[0]][3])/((double)order-1); { if ((PresentNode(F2))) { if (! net){ assert(! F2->exists); } f2[0] = NewNumVer(F2); newco.el[f2[0]][0] = dx; newco.el[f2[0]][1] = dy; newco.el[f2[0]][2] = dz; newco.el[f2[0]][3] = dw; } } ??? dx = (newco.el[e4[0]][0] + newco.el[e5[0]][0])/((double)order-1); ??? dy = (newco.el[e4[0]][1] + newco.el[e5[0]][1])/((double)order-1); ??? dz = (newco.el[e4[0]][2] + newco.el[e5[0]][2])/((double)order-1); ??? dw = (newco.el[e4[0]][3] + newco.el[e5[0]][3])/((double)order-1); { if ((PresentNode(F3))) { if (! net){ assert(! F3->exists); } f3[0] = NewNumVer(F3); newco.el[f3[0]][0] = dx; newco.el[f3[0]][1] = dy; newco.el[f3[0]][2] = dz; newco.el[f3[0]][3] = dw; } } } } else if (order == 4){ Node_t F00 = OrgV(ncor[i].right[3]); ??? F01 = OrgV(ncor[i].right[8]); ??? F02 = OrgV(ncor[i].right[6]); ??? F10 = OrgV(ncor[i].left[3]); ??? F11 = OrgV(ncor[i].left[8]); ??? F12 = OrgV(ncor[i].left[6]); ??? F30 = OrgV(ncor[i].back[3]); ??? F31 = OrgV(ncor[i].back[8]); ??? F32 = OrgV(ncor[i].back[6]); ??? F20 = OrgV(ncor[i].front[3]); ??? F21 = OrgV(ncor[i].front[8]); ??? F22 = OrgV(ncor[i].front[6]); { assert0 == strcmp(F00.label,"VF")); assert0 == strcmp(F01.label,"VF")); assert0 == strcmp(F02.label,"VF")); assert0 == strcmp(F10.label,"VF")); assert0 == strcmp(F11.label,"VF")); assert0 == strcmp(F12.label,"VF")); assert0 == strcmp(F30.label,"VF")); assert0 == strcmp(F31.label,"VF")); assert0 == strcmp(F32.label,"VF")); assert0 == strcmp(F20.label,"VF")); assert0 == strcmp(F21.label,"VF")); assert0 == strcmp(F22.label,"VF")); /* Right */ ??? dx = (newco.el[e1[0]][0] - newco.el[e2[0]][0])/((double)order-1); ??? dy = (newco.el[e1[0]][1] - newco.el[e2[0]][1])/((double)order-1); ??? dz = (newco.el[e1[0]][2] - newco.el[e2[0]][2])/((double)order-1); ??? dw = (newco.el[e1[0]][3] - newco.el[e2[0]][3])/((double)order-1); { if ((PresentNode(F00))) { if (! net){ assert(! F00->exists); } f0[0] = NewNumVer(F00); newco.el[f0[0]][0] = newco.el[e1[0]][0] - 1.0 * dx; newco.el[f0[0]][1] = newco.el[e1[0]][1] - 1.0 * dy; newco.el[f0[0]][2] = newco.el[e1[0]][2] - 1.0 * dz; newco.el[f0[0]][3] = newco.el[e1[0]][3] - 1.0 * dw; } if ((PresentNode(F02))) { if (! net){ assert(! F02->exists); } f0[2] = NewNumVer(F02); newco.el[f0[2]][0] = newco.el[e1[0]][0] - 2.0 * dx; newco.el[f0[2]][1] = newco.el[e1[0]][1] - 2.0 * dy; newco.el[f0[2]][2] = newco.el[e1[0]][2] - 2.0 * dz; newco.el[f0[2]][3] = newco.el[e1[0]][3] - 2.0 * dw; } } ??? dx = (newco.el[e1[1]][0] - newco.el[e2[1]][0])/((double)order-2); ??? dy = (newco.el[e1[1]][1] - newco.el[e2[1]][1])/((double)order-2); ??? dz = (newco.el[e1[1]][2] - newco.el[e2[1]][2])/((double)order-2); ??? dw = (newco.el[e1[1]][3] - newco.el[e2[1]][3])/((double)order-2); { if ((PresentNode(F01))) { if (! net){ assert(! F01->exists); } f0[1] = NewNumVer(F01); newco.el[f0[1]][0] = newco.el[e1[1]][0] - 1.0 * dx; newco.el[f0[1]][1] = newco.el[e1[1]][1] - 1.0 * dy; newco.el[f0[1]][2] = newco.el[e1[1]][2] - 1.0 * dz; newco.el[f0[1]][3] = newco.el[e1[1]][3] - 1.0 * dw; } } /* Left */ ??? dx = (newco.el[e0[0]][0] - newco.el[e4[2]][0])/((double)order-1); ??? dy = (newco.el[e0[0]][1] - newco.el[e4[2]][1])/((double)order-1); ??? dz = (newco.el[e0[0]][2] - newco.el[e4[2]][2])/((double)order-1); ??? dw = (newco.el[e0[0]][3] - newco.el[e4[2]][3])/((double)order-1); { if ((PresentNode(F10))) { if (! net){ assert(! F10->exists); } f1[0] = NewNumVer(F10); newco.el[f1[0]][0] = newco.el[e0[0]][0] - 1.0 * dx; newco.el[f1[0]][1] = newco.el[e0[0]][1] - 1.0 * dy; newco.el[f1[0]][2] = newco.el[e0[0]][2] - 1.0 * dz; newco.el[f1[0]][3] = newco.el[e0[0]][3] - 1.0 * dw; } if ((PresentNode(F11))) { if (! net){ assert(! F11->exists); } f1[1] = NewNumVer(F11); newco.el[f1[1]][0] = newco.el[e0[0]][0] - 2.0 * dx; newco.el[f1[1]][1] = newco.el[e0[0]][1] - 2.0 * dy; newco.el[f1[1]][2] = newco.el[e0[0]][2] - 2.0 * dz; newco.el[f1[1]][3] = newco.el[e0[0]][3] - 2.0 * dw; } } ??? dx = (newco.el[e0[1]][0] - newco.el[e4[1]][0])/((double)order-2); ??? dy = (newco.el[e0[1]][1] - newco.el[e4[1]][1])/((double)order-2); ??? dz = (newco.el[e0[1]][2] - newco.el[e4[1]][2])/((double)order-2); ??? dw = (newco.el[e0[1]][3] - newco.el[e4[1]][3])/((double)order-2); { if ((PresentNode(F12))) { if (! net){ assert(! F12->exists); } f1[2] = NewNumVer(F12); newco.el[f1[2]][0] = newco.el[e0[1]][0] - 1.0 * dx; newco.el[f1[2]][1] = newco.el[e0[1]][1] - 1.0 * dy; newco.el[f1[2]][2] = newco.el[e0[1]][2] - 1.0 * dz; newco.el[f1[2]][3] = newco.el[e0[1]][3] - 1.0 * dw; } } /* Front */ ??? dx = (newco.el[e3[0]][0] - newco.el[e5[0]][0])/((double)order-1); ??? dy = (newco.el[e3[0]][1] - newco.el[e5[0]][1])/((double)order-1); ??? dz = (newco.el[e3[0]][2] - newco.el[e5[0]][2])/((double)order-1); ??? dw = (newco.el[e3[0]][3] - newco.el[e5[0]][3])/((double)order-1); { if ((PresentNode(F20))) { if (! net){ assert(! F20->exists); } f2[0] = NewNumVer(F20); newco.el[f2[0]][0] = newco.el[e3[0]][0] - 1.0 * dx; newco.el[f2[0]][1] = newco.el[e3[0]][1] - 1.0 * dy; newco.el[f2[0]][2] = newco.el[e3[0]][2] - 1.0 * dz; newco.el[f2[0]][3] = newco.el[e3[0]][3] - 1.0 * dw; } if ((PresentNode(F21))) { if (! net){ assert(! F21->exists); } f2[1] = NewNumVer(F21); newco.el[f2[1]][0] = newco.el[e3[0]][0] - 2.0 * dx; newco.el[f2[1]][1] = newco.el[e3[0]][1] - 2.0 * dy; newco.el[f2[1]][2] = newco.el[e3[0]][2] - 2.0 * dz; newco.el[f2[1]][3] = newco.el[e3[0]][3] - 2.0 * dw; } } ??? dx = (newco.el[e3[1]][0] - newco.el[e5[1]][0])/((double)order-2); ??? dy = (newco.el[e3[1]][1] - newco.el[e5[1]][1])/((double)order-2); ??? dz = (newco.el[e3[1]][2] - newco.el[e5[1]][2])/((double)order-2); ??? dw = (newco.el[e3[1]][3] - newco.el[e5[1]][3])/((double)order-2); { if ((PresentNode(F22))) { if (! net){ assert(! F22->exists); } f2[2] = NewNumVer(F22); newco.el[f2[2]][0] = newco.el[e3[1]][0] - 1.0 * dx; newco.el[f2[2]][1] = newco.el[e3[1]][1] - 1.0 * dy; newco.el[f2[2]][2] = newco.el[e3[1]][2] - 1.0 * dz; newco.el[f2[2]][3] = newco.el[e3[1]][3] - 1.0 * dw; } } /* Back */ ??? dx = (newco.el[e4[0]][0] - newco.el[e5[0]][0])/((double)order-1); ??? dy = (newco.el[e4[0]][1] - newco.el[e5[0]][1])/((double)order-1); ??? dz = (newco.el[e4[0]][2] - newco.el[e5[0]][2])/((double)order-1); ??? dw = (newco.el[e4[0]][3] - newco.el[e5[0]][3])/((double)order-1); { if ((PresentNode(F30))) { if (! net){ assert(! F30->exists); } f3[0] = NewNumVer(F30); newco.el[f3[0]][0] = newco.el[e4[0]][0] - 1.0 * dx; newco.el[f3[0]][1] = newco.el[e4[0]][1] - 1.0 * dy; newco.el[f3[0]][2] = newco.el[e4[0]][2] - 1.0 * dz; newco.el[f3[0]][3] = newco.el[e4[0]][3] - 1.0 * dw; } if ((PresentNode(F31))) { if (! net){ assert(! F31->exists); } f3[1] = NewNumVer(F31); newco.el[f3[1]][0] = newco.el[e4[0]][0] - 2.0 * dx; newco.el[f3[1]][1] = newco.el[e4[0]][1] - 2.0 * dy; newco.el[f3[1]][2] = newco.el[e4[0]][2] - 2.0 * dz; newco.el[f3[1]][3] = newco.el[e4[0]][3] - 2.0 * dw; } } ??? dx = (newco.el[e4[1]][0] - newco.el[e5[1]][0])/((double)order-2); ??? dy = (newco.el[e4[1]][1] - newco.el[e5[1]][1])/((double)order-2); ??? dz = (newco.el[e4[1]][2] - newco.el[e5[1]][2])/((double)order-2); ??? dw = (newco.el[e4[1]][3] - newco.el[e5[1]][3])/((double)order-2); { if ((PresentNode(F32))) { if (! net){ assert(! F32->exists); } f3[2] = NewNumVer(F32); newco.el[f3[2]][0] = newco.el[e4[1]][0] - 1.0 * dx; newco.el[f3[2]][1] = newco.el[e4[1]][1] - 1.0 * dy; newco.el[f3[2]][2] = newco.el[e4[1]][2] - 1.0 * dz; newco.el[f3[2]][3] = newco.el[e4[1]][3] - 1.0 * dw; } } } } else if (order == 5){ Node_t F00 = OrgV(ncor[i].right[ 3]); ??? F01 = OrgV(ncor[i].right[ 8]); ??? F02 = OrgV(ncor[i].right[ 6]); ??? F03 = OrgV(ncor[i].right[15]); ??? F04 = OrgV(ncor[i].right[13]); ??? F05 = OrgV(ncor[i].right[11]); ??? F10 = OrgV(ncor[i].left [ 3]); ??? F11 = OrgV(ncor[i].left [ 8]); ??? F12 = OrgV(ncor[i].left [ 6]); ??? F13 = OrgV(ncor[i].left [15]); ??? F14 = OrgV(ncor[i].left [13]); ??? F15 = OrgV(ncor[i].left [11]); ??? F30 = OrgV(ncor[i].back [ 3]); ??? F31 = OrgV(ncor[i].back [ 8]); ??? F32 = OrgV(ncor[i].back [ 6]); ??? F33 = OrgV(ncor[i].back [15]); ??? F34 = OrgV(ncor[i].back [13]); ??? F35 = OrgV(ncor[i].back [11]); ??? F20 = OrgV(ncor[i].front[ 3]); ??? F21 = OrgV(ncor[i].front[ 8]); ??? F22 = OrgV(ncor[i].front[ 6]); ??? F23 = OrgV(ncor[i].front[15]); ??? F24 = OrgV(ncor[i].front[13]); ??? F25 = OrgV(ncor[i].front[11]); { /* Right */ ??? dx = (newco.el[e1[0]][0] - newco.el[e2[0]][0])/((double)order-1); ??? dy = (newco.el[e1[0]][1] - newco.el[e2[0]][1])/((double)order-1); ??? dz = (newco.el[e1[0]][2] - newco.el[e2[0]][2])/((double)order-1); ??? dw = (newco.el[e1[0]][3] - newco.el[e2[0]][3])/((double)order-1); { if ((PresentNode(F00))) { if (! net){ assert(! F00->exists); } f0[0] = NewNumVer(F00); newco.el[f0[0]][0] = newco.el[e1[0]][0] - 1.0 * dx; newco.el[f0[0]][1] = newco.el[e1[0]][1] - 1.0 * dy; newco.el[f0[0]][2] = newco.el[e1[0]][2] - 1.0 * dz; newco.el[f0[0]][3] = newco.el[e1[0]][3] - 1.0 * dw; } if ((PresentNode(F02))) { if (! net){ assert(! F02->exists); } f0[2] = NewNumVer(F02); newco.el[f0[2]][0] = newco.el[e1[0]][0] - 2.0 * dx; newco.el[f0[2]][1] = newco.el[e1[0]][1] - 2.0 * dy; newco.el[f0[2]][2] = newco.el[e1[0]][2] - 2.0 * dz; newco.el[f0[2]][3] = newco.el[e1[0]][3] - 2.0 * dw; } if ((PresentNode(F05))) { if (! net){ assert(! F05->exists); } f0[5] = NewNumVer(F05); newco.el[f0[5]][0] = newco.el[e1[0]][0] - 3.0 * dx; newco.el[f0[5]][1] = newco.el[e1[0]][1] - 3.0 * dy; newco.el[f0[5]][2] = newco.el[e1[0]][2] - 3.0 * dz; newco.el[f0[5]][3] = newco.el[e1[0]][3] - 3.0 * dw; } } ??? dx = (newco.el[e1[1]][0] - newco.el[e2[1]][0])/((double)order-2); ??? dy = (newco.el[e1[1]][1] - newco.el[e2[1]][1])/((double)order-2); ??? dz = (newco.el[e1[1]][2] - newco.el[e2[1]][2])/((double)order-2); ??? dw = (newco.el[e1[1]][3] - newco.el[e2[1]][3])/((double)order-2); { if ((PresentNode(F01))) { if (! net){ assert(! F01->exists); } f0[1] = NewNumVer(F01); newco.el[f0[1]][0] = newco.el[e1[1]][0] - 1.0 * dx; newco.el[f0[1]][1] = newco.el[e1[1]][1] - 1.0 * dy; newco.el[f0[1]][2] = newco.el[e1[1]][2] - 1.0 * dz; newco.el[f0[1]][3] = newco.el[e1[1]][3] - 1.0 * dw; } if ((PresentNode(F04))) { if (! net){ assert(! F04->exists); } f0[4] = NewNumVer(F04); newco.el[f0[4]][0] = newco.el[e1[1]][0] - 2.0 * dx; newco.el[f0[4]][1] = newco.el[e1[1]][1] - 2.0 * dy; newco.el[f0[4]][2] = newco.el[e1[1]][2] - 2.0 * dz; newco.el[f0[4]][3] = newco.el[e1[1]][3] - 2.0 * dw; } } ??? dx = (newco.el[e1[2]][0] - newco.el[e2[2]][0])/((double)order-3); ??? dy = (newco.el[e1[2]][1] - newco.el[e2[2]][1])/((double)order-3); ??? dz = (newco.el[e1[2]][2] - newco.el[e2[2]][2])/((double)order-3); ??? dw = (newco.el[e1[2]][3] - newco.el[e2[2]][3])/((double)order-3); { if ((PresentNode(F03))) { if (! net){ assert(! F03->exists); } f0[3] = NewNumVer(F03); newco.el[f0[3]][0] = newco.el[e1[2]][0] - 1.0 * dx; newco.el[f0[3]][1] = newco.el[e1[2]][1] - 1.0 * dy; newco.el[f0[3]][2] = newco.el[e1[2]][2] - 1.0 * dz; newco.el[f0[3]][3] = newco.el[e1[2]][3] - 1.0 * dw; } } /* Left */ ??? dx = (newco.el[e0[0]][0] - newco.el[e4[3]][0])/((double)order-1); ??? dy = (newco.el[e0[0]][1] - newco.el[e4[3]][1])/((double)order-1); ??? dz = (newco.el[e0[0]][2] - newco.el[e4[3]][2])/((double)order-1); ??? dw = (newco.el[e0[0]][3] - newco.el[e4[3]][3])/((double)order-1); { if ((PresentNode(F10))) { if (! net){ assert(! F10->exists); } f1[0] = NewNumVer(F10); newco.el[f1[0]][0] = newco.el[e0[0]][0] - 1.0 * dx; newco.el[f1[0]][1] = newco.el[e0[0]][1] - 1.0 * dy; newco.el[f1[0]][2] = newco.el[e0[0]][2] - 1.0 * dz; newco.el[f1[0]][3] = newco.el[e0[0]][3] - 1.0 * dw; } if ((PresentNode(F11))) { if (! net){ assert(! F11->exists); } f1[1] = NewNumVer(F11); newco.el[f1[1]][0] = newco.el[e0[0]][0] - 2.0 * dx; newco.el[f1[1]][1] = newco.el[e0[0]][1] - 2.0 * dy; newco.el[f1[1]][2] = newco.el[e0[0]][2] - 2.0 * dz; newco.el[f1[1]][3] = newco.el[e0[0]][3] - 2.0 * dw; } if ((PresentNode(F13))) { if (! net){ assert(! F13->exists); } f1[3] = NewNumVer(F13); newco.el[f1[3]][0] = newco.el[e0[0]][0] - 3.0 * dx; newco.el[f1[3]][1] = newco.el[e0[0]][1] - 3.0 * dy; newco.el[f1[3]][2] = newco.el[e0[0]][2] - 3.0 * dz; newco.el[f1[3]][3] = newco.el[e0[0]][3] - 3.0 * dw; } } ??? dx = (newco.el[e0[1]][0] - newco.el[e4[2]][0])/((double)order-2); ??? dy = (newco.el[e0[1]][1] - newco.el[e4[2]][1])/((double)order-2); ??? dz = (newco.el[e0[1]][2] - newco.el[e4[2]][2])/((double)order-2); ??? dw = (newco.el[e0[1]][3] - newco.el[e4[2]][3])/((double)order-2); { if ((PresentNode(F12))) { if (! net){ assert(! F12->exists); } f1[2] = NewNumVer(F12); newco.el[f1[2]][0] = newco.el[e0[1]][0] - 1.0 * dx; newco.el[f1[2]][1] = newco.el[e0[1]][1] - 1.0 * dy; newco.el[f1[2]][2] = newco.el[e0[1]][2] - 1.0 * dz; newco.el[f1[2]][3] = newco.el[e0[1]][3] - 1.0 * dw; } if ((PresentNode(F14))) { if (! net){ assert(! F14->exists); } f1[4] = NewNumVer(F14); newco.el[f1[4]][0] = newco.el[e0[1]][0] - 2.0 * dx; newco.el[f1[4]][1] = newco.el[e0[1]][1] - 2.0 * dy; newco.el[f1[4]][2] = newco.el[e0[1]][2] - 2.0 * dz; newco.el[f1[4]][3] = newco.el[e0[1]][3] - 2.0 * dw; } } ??? dx = (newco.el[e0[2]][0] - newco.el[e4[1]][0])/((double)order-3); ??? dy = (newco.el[e0[2]][1] - newco.el[e4[1]][1])/((double)order-3); ??? dz = (newco.el[e0[2]][2] - newco.el[e4[1]][2])/((double)order-3); ??? dw = (newco.el[e0[2]][3] - newco.el[e4[1]][3])/((double)order-3); { if ((PresentNode(F15))) { if (! net){ assert(! F15->exists); } f1[5] = NewNumVer(F15); newco.el[f1[5]][0] = newco.el[e0[2]][0] - 1.0 * dx; newco.el[f1[5]][1] = newco.el[e0[2]][1] - 1.0 * dy; newco.el[f1[5]][2] = newco.el[e0[2]][2] - 1.0 * dz; newco.el[f1[5]][3] = newco.el[e0[2]][3] - 1.0 * dw; } } /* Front */ ??? dx = (newco.el[e3[0]][0] - newco.el[e5[0]][0])/((double)order-1); ??? dy = (newco.el[e3[0]][1] - newco.el[e5[0]][1])/((double)order-1); ??? dz = (newco.el[e3[0]][2] - newco.el[e5[0]][2])/((double)order-1); ??? dw = (newco.el[e3[0]][3] - newco.el[e5[0]][3])/((double)order-1); { if ((PresentNode(F20))) { if (! net){ assert(! F20->exists); } f2[0] = NewNumVer(F20); newco.el[f2[0]][0] = newco.el[e3[0]][0] - 1.0 * dx; newco.el[f2[0]][1] = newco.el[e3[0]][1] - 1.0 * dy; newco.el[f2[0]][2] = newco.el[e3[0]][2] - 1.0 * dz; newco.el[f2[0]][3] = newco.el[e3[0]][3] - 1.0 * dw; } if ((PresentNode(F21))) { if (! net){ assert(! F21->exists); } f2[1] = NewNumVer(F21); newco.el[f2[1]][0] = newco.el[e3[0]][0] - 2.0 * dx; newco.el[f2[1]][1] = newco.el[e3[0]][1] - 2.0 * dy; newco.el[f2[1]][2] = newco.el[e3[0]][2] - 2.0 * dz; newco.el[f2[1]][3] = newco.el[e3[0]][3] - 2.0 * dw; } if ((PresentNode(F23))) { if (! net){ assert(! F23->exists); } f2[3] = NewNumVer(F23); newco.el[f2[3]][0] = newco.el[e3[0]][0] - 3.0 * dx; newco.el[f2[3]][1] = newco.el[e3[0]][1] - 3.0 * dy; newco.el[f2[3]][2] = newco.el[e3[0]][2] - 3.0 * dz; newco.el[f2[3]][3] = newco.el[e3[0]][3] - 3.0 * dw; } } ??? dx = (newco.el[e3[1]][0] - newco.el[e5[1]][0])/((double)order-2); ??? dy = (newco.el[e3[1]][1] - newco.el[e5[1]][1])/((double)order-2); ??? dz = (newco.el[e3[1]][2] - newco.el[e5[1]][2])/((double)order-2); ??? dw = (newco.el[e3[1]][3] - newco.el[e5[1]][3])/((double)order-2); { if ((PresentNode(F22))) { if (! net){ assert(! F22->exists); } f2[2] = NewNumVer(F22); newco.el[f2[2]][0] = newco.el[e3[1]][0] - 1.0 * dx; newco.el[f2[2]][1] = newco.el[e3[1]][1] - 1.0 * dy; newco.el[f2[2]][2] = newco.el[e3[1]][2] - 1.0 * dz; newco.el[f2[2]][3] = newco.el[e3[1]][3] - 1.0 * dw; } if ((PresentNode(F24))) { if (! net){ assert(! F24->exists); } f2[4] = NewNumVer(F24); newco.el[f2[4]][0] = newco.el[e3[1]][0] - 2.0 * dx; newco.el[f2[4]][1] = newco.el[e3[1]][1] - 2.0 * dy; newco.el[f2[4]][2] = newco.el[e3[1]][2] - 2.0 * dz; newco.el[f2[4]][3] = newco.el[e3[1]][3] - 2.0 * dw; } } ??? dx = (newco.el[e3[2]][0] - newco.el[e5[2]][0])/((double)order-3); ??? dy = (newco.el[e3[2]][1] - newco.el[e5[2]][1])/((double)order-3); ??? dz = (newco.el[e3[2]][2] - newco.el[e5[2]][2])/((double)order-3); ??? dw = (newco.el[e3[2]][3] - newco.el[e5[2]][3])/((double)order-3); { if ((PresentNode(F25))) { if (! net){ assert(! F25->exists); } f2[5] = NewNumVer(F25); newco.el[f2[5]][0] = newco.el[e3[2]][0] - 1.0 * dx; newco.el[f2[5]][1] = newco.el[e3[2]][1] - 1.0 * dy; newco.el[f2[5]][2] = newco.el[e3[2]][2] - 1.0 * dz; newco.el[f2[5]][3] = newco.el[e3[2]][3] - 1.0 * dw; } } /* Back */ ??? dx = (newco.el[e4[0]][0] - newco.el[e5[0]][0])/((double)order-1); ??? dy = (newco.el[e4[0]][1] - newco.el[e5[0]][1])/((double)order-1); ??? dz = (newco.el[e4[0]][2] - newco.el[e5[0]][2])/((double)order-1); ??? dw = (newco.el[e4[0]][3] - newco.el[e5[0]][3])/((double)order-1); { if ((PresentNode(F30))) { if (! net){ assert(! F30->exists); } f3[0] = NewNumVer(F30); newco.el[f3[0]][0] = newco.el[e4[0]][0] - 1.0 * dx; newco.el[f3[0]][1] = newco.el[e4[0]][1] - 1.0 * dy; newco.el[f3[0]][2] = newco.el[e4[0]][2] - 1.0 * dz; newco.el[f3[0]][3] = newco.el[e4[0]][3] - 1.0 * dw; } if ((PresentNode(F31))) { if (! net){ assert(! F31->exists); } f3[1] = NewNumVer(F31); newco.el[f3[1]][0] = newco.el[e4[0]][0] - 2.0 * dx; newco.el[f3[1]][1] = newco.el[e4[0]][1] - 2.0 * dy; newco.el[f3[1]][2] = newco.el[e4[0]][2] - 2.0 * dz; newco.el[f3[1]][3] = newco.el[e4[0]][3] - 2.0 * dw; } if ((PresentNode(F33))) { if (! net){ assert(! F33->exists); } f3[3] = NewNumVer(F33); newco.el[f3[3]][0] = newco.el[e4[0]][0] - 3.0 * dx; newco.el[f3[3]][1] = newco.el[e4[0]][1] - 3.0 * dy; newco.el[f3[3]][2] = newco.el[e4[0]][2] - 3.0 * dz; newco.el[f3[3]][3] = newco.el[e4[0]][3] - 3.0 * dw; } } ??? dx = (newco.el[e4[1]][0] - newco.el[e5[1]][0])/((double)order-2); ??? dy = (newco.el[e4[1]][1] - newco.el[e5[1]][1])/((double)order-2); ??? dz = (newco.el[e4[1]][2] - newco.el[e5[1]][2])/((double)order-2); ??? dw = (newco.el[e4[1]][3] - newco.el[e5[1]][3])/((double)order-2); { if ((PresentNode(F32))) { if (! net){ assert(! F32->exists); } f3[2] = NewNumVer(F32); newco.el[f3[2]][0] = newco.el[e4[1]][0] - 1.0 * dx; newco.el[f3[2]][1] = newco.el[e4[1]][1] - 1.0 * dy; newco.el[f3[2]][2] = newco.el[e4[1]][2] - 1.0 * dz; newco.el[f3[2]][3] = newco.el[e4[1]][3] - 1.0 * dw; } if ((PresentNode(F34))) { if (! net){ assert(! F34->exists); } f3[4] = NewNumVer(F34); newco.el[f3[4]][0] = newco.el[e4[1]][0] - 2.0 * dx; newco.el[f3[4]][1] = newco.el[e4[1]][1] - 2.0 * dy; newco.el[f3[4]][2] = newco.el[e4[1]][2] - 2.0 * dz; newco.el[f3[4]][3] = newco.el[e4[1]][3] - 2.0 * dw; } } ??? dx = (newco.el[e4[2]][0] - newco.el[e5[2]][0])/((double)order-3); ??? dy = (newco.el[e4[2]][1] - newco.el[e5[2]][1])/((double)order-3); ??? dz = (newco.el[e4[2]][2] - newco.el[e5[2]][2])/((double)order-3); ??? dw = (newco.el[e4[2]][3] - newco.el[e5[2]][3])/((double)order-3); { if ((PresentNode(F35))) { if (! net){ assert(! F35->exists); } f3[5] = NewNumVer(F35); newco.el[f3[5]][0] = newco.el[e4[2]][0] - 1.0 * dx; newco.el[f3[5]][1] = newco.el[e4[2]][1] - 1.0 * dy; newco.el[f3[5]][2] = newco.el[e4[2]][2] - 1.0 * dz; newco.el[f3[5]][3] = newco.el[e4[2]][3] - 1.0 * dw; } } } } /* nodes on the octahedron barycenter. */ if (order == 2) { if ((PresentNode(mocta[i][0]))) { if (! net){ assert(! mocta[i][0]->exists); } ??? n = NewNumVer(mocta[i][0]); ??? b0 = Add(newco.el[e3[0]],newco.el[e1[0]]); ??? b1 = Add(newco.el[e4[0]],newco.el[e2[0]]); ??? b = Add(b0,b1); { newco.el[n] = Scale(1.0/((double)4), b); } } } else if (order == 3){ if ((PresentNode(mocta[i][0]))) { if (! net){ assert(! mocta[i][0]->exists); } ??? n = NewNumVer(mocta[i][0]); ??? b0 = Add(Add(newco.el[f2[0]]; with (newco.el[e3[0]]),newco.el[f1[0]]), ??? b1 = Add(Add(newco.el[e1[0]],newco.el[e0[0]]),newco.el[f0[0]]); ??? b = Add(b0,b1); { newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][1]))) { if (! net){ assert(! mocta[i][1]->exists); } ??? n = NewNumVer(mocta[i][1]); ??? b0 = Add(Add(newco.el[f2[0]],newco.el[f3[0]]),newco.el[e5[1]]); ??? b1 = Add(Add(newco.el[f1[0]],newco.el[e4[1]]),newco.el[e3[1]]); ??? b = Add(b0,b1); { MyAssert(newco.el[f2[0]]); MyAssert(newco.el[f3[0]]); MyAssert(newco.el[e5[1]]); MyAssert(newco.el[f1[0]]); MyAssert(newco.el[e4[1]]); MyAssert(newco.el[e3[1]]); newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][2]))) { if (! net){ assert(! mocta[i][2]->exists); } ??? n = NewNumVer(mocta[i][2]); ??? b0 = Add(Add(newco.el[f2[0]]; with (newco.el[f3[0]]),newco.el[e5[0]]), double b1 = Add(Add(newco.el[e1[1]],newco.el[f0[0]]),newco.el[e2[1]]); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][3]))) { if (! net){ assert(! mocta[i][3]->exists); } ??? n = NewNumVer(mocta[i][3]); ??? b0 = Add(Add(newco.el[f3[0]]; with (newco.el[f1[0]]),newco.el[e4[0]]), double b1 = Add(Add(newco.el[f0[0]],newco.el[e0[1]]),newco.el[e2[0]]); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } } else if (order == 4){ if ((PresentNode(mocta[i][0]))) { if (! net){ assert(! mocta[i][0]->exists); } ??? n = NewNumVer(mocta[i][0]); ??? b0 = Add(Add(newco.el[e1[0]]; with (newco.el[e0[0]]),newco.el[f0[0]]), double b1 = Add(Add(newco.el[f1[0]],newco.el[e3[0]]),newco.el[f2[0]]); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][3]))) { if (! net){ assert(! mocta[i][3]->exists); } ??? n = NewNumVer(mocta[i][3]); ??? b0 = Add(Add(newco.el[e3[2]]; with (newco.el[f2[2]]),newco.el[e5[2]]), double b1 = Add(Add(newco.el[f1[1]],newco.el[e4[2]]),newco.el[f3[2]]); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][5]))) { if (! net){ assert(! mocta[i][5]->exists); } ??? n = NewNumVer(mocta[i][5]); ??? b0 = Add(Add(newco.el[f2[1]]; with (newco.el[e1[2]]),newco.el[e5[0]]), double b1 = Add(Add(newco.el[f3[1]],newco.el[e2[2]]),newco.el[f0[1]]); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][6]))) { if (! net){ assert(! mocta[i][6]->exists); } ??? n = NewNumVer(mocta[i][6]); ??? b0 = Add(Add(newco.el[f0[2]]; with (newco.el[e0[2]]),newco.el[e2[0]]), double b1 = Add(Add(newco.el[f1[2]],newco.el[e4[0]]),newco.el[f3[0]]); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } ??? dx = (newco.el[f1[1]][0] - newco.el[f0[1]][0])/((double)order-2); ??? dy = (newco.el[f1[1]][1] - newco.el[f0[1]][1])/((double)order-2); ??? dz = (newco.el[f1[1]][2] - newco.el[f0[1]][2])/((double)order-2); ??? dw = (newco.el[f1[1]][3] - newco.el[f0[1]][3])/((double)order-2); ??? p0 = newco.el[f1[1]][0] - 1.0 * dx; ??? p1 = newco.el[f1[1]][1] - 1.0 * dy; ??? p2 = newco.el[f1[1]][2] - 1.0 * dz; ??? p3 = newco.el[f1[1]][3] - 1.0 * dw; with ( p == (r4_t){p0,p1,p2,p3} ){ if ((PresentNode(inter[i][0]))) { assert( NOT inter[i][0]->exists); ??? n = NewNumVer(inter[i][0]); { newco.el[n] = p; } } if ((PresentNode(mocta[i][1]))) { if (! net){ assert(! mocta[i][1]->exists); } ??? n = NewNumVer(mocta[i][1]); ??? b0 = Add(Add(newco.el[e3[1]]; with (newco.el[f2[0]]),newco.el[f2[2]]), double b1 = Add(Add(newco.el[f1[0]],newco.el[f1[1]]),p); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][2]))) { if (! net){ assert(! mocta[i][2]->exists); } ??? n = NewNumVer(mocta[i][2]); ??? b0 = Add(Add(newco.el[f2[0]]; with (newco.el[e1[1]]),newco.el[f2[1]]), double b1 = Add(Add(newco.el[f0[0]],newco.el[f0[1]]),p); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][4]))) { if (! net){ assert(! mocta[i][4]->exists); } ??? n = NewNumVer(mocta[i][4]); ??? b0 = Add(Add(newco.el[f2[2]]; with (newco.el[f2[1]]),newco.el[e5[1]]), double b1 = Add(Add(newco.el[f3[2]],newco.el[f3[1]]),p); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][7]))) { if (! net){ assert(! mocta[i][7]->exists); } ??? n = NewNumVer(mocta[i][7]); ??? b0 = Add(Add(newco.el[f0[0]]; with (newco.el[e0[1]]),newco.el[f0[2]]), double b1 = Add(Add(newco.el[f1[0]],newco.el[f1[2]]),p); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][8]))) { if (! net){ assert(! mocta[i][8]->exists); } ??? n = NewNumVer(mocta[i][8]); ??? b0 = Add(Add(newco.el[f0[1]]; with (newco.el[f0[2]]),newco.el[e2[1]]), double b1 = Add(Add(newco.el[f3[0]],newco.el[f3[1]]),p); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][9]))) { if (! net){ assert(! mocta[i][9]->exists); } ??? n = NewNumVer(mocta[i][9]); ??? b0 = Add(Add(newco.el[f1[2]]; with (newco.el[f1[1]]),newco.el[e4[1]]), double b1 = Add(Add(newco.el[f3[0]],newco.el[f3[2]]),p); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } } } else if (order == 5){ if ((PresentNode(mocta[i][0]))) { if (! net){ assert(! mocta[i][0]->exists); } ??? n = NewNumVer(mocta[i][0]); ??? b0 = Add(Add(newco.el[e1[0]]; with (newco.el[e0[0]]),newco.el[f0[0]]), double b1 = Add(Add(newco.el[f1[0]],newco.el[e3[0]]),newco.el[f2[0]]); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][6]))) { if (! net){ assert(! mocta[i][6]->exists); } ??? n = NewNumVer(mocta[i][6]); ??? b0 = Add(Add(newco.el[f1[3]]; with (newco.el[e3[3]]),newco.el[e4[3]]), double b1 = Add(Add(newco.el[f3[5]],newco.el[e5[3]]),newco.el[f2[5]]); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][9]))) { if (! net){ assert(! mocta[i][9]->exists); } ??? n = NewNumVer(mocta[i][9]); ??? b0 = Add(Add(newco.el[e1[3]]; with (newco.el[f0[3]]),newco.el[e2[3]]), double b1 = Add(Add(newco.el[f2[3]],newco.el[e5[0]]),newco.el[f3[3]]); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][10]))) { if (! net){ assert(! mocta[i][10]->exists); } ??? n = NewNumVer(mocta[i][10]); ??? b0 = Add(Add(newco.el[f0[5]]; with (newco.el[e0[3]]),newco.el[e2[0]]), double b1 = Add(Add(newco.el[f1[5]],newco.el[e4[0]]),newco.el[f3[0]]); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } ??? dx = (newco.el[f1[2]][0] - newco.el[f2[1]][0])/((double)order-3); ??? dy = (newco.el[f1[2]][1] - newco.el[f2[1]][1])/((double)order-3); ??? dz = (newco.el[f1[2]][2] - newco.el[f2[1]][2])/((double)order-3); ??? dw = (newco.el[f1[2]][3] - newco.el[f2[1]][3])/((double)order-3); ??? Q0 = newco.el[f1[2]][0] - 1.0 * dx; ??? Q1 = newco.el[f1[2]][1] - 1.0 * dy; ??? Q2 = newco.el[f1[2]][2] - 1.0 * dz; ??? Q3 = newco.el[f1[2]][3] - 1.0 * dw; with ( double q1 = (r4_t){Q0,Q1,Q2,Q3} ){ if ((PresentNode(mocta[i][1]))) { if (! net){ assert(! mocta[i][1]->exists); } ??? n = NewNumVer(mocta[i][1]); ??? b0 = Add(Add(newco.el[e3[1]]; with (newco.el[f2[0]]),newco.el[f2[2]]), double b1 = Add(Add(newco.el[f1[0]],newco.el[f1[1]]),q1); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][2]))) { if (! net){ assert(! mocta[i][2]->exists); } ??? n = NewNumVer(mocta[i][2]); ??? b0 = Add(Add(newco.el[e1[1]]; with (newco.el[f0[0]]),newco.el[f0[1]]), double b1 = Add(Add(newco.el[f2[0]],newco.el[f2[1]]),q1); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][13]))) { if (! net){ assert(! mocta[i][13]->exists); } ??? n = NewNumVer(mocta[i][13]); ??? b0 = Add(Add(newco.el[f0[0]]; with (newco.el[e0[1]]),newco.el[f0[2]]), double b1 = Add(Add(newco.el[f1[0]],newco.el[f1[2]]),q1); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } /* first internal node: q1 */ if ((PresentNode(inter[i][0]))) { assert( NOT inter[i][0]->exists); ??? n = NewNumVer(inter[i][0]); { newco.el[n] = q1; } } ??? dx = (newco.el[f1[5]][0] - newco.el[f2[3]][0])/((double)order-2); ??? dy = (newco.el[f1[5]][1] - newco.el[f2[3]][1])/((double)order-2); ??? dz = (newco.el[f1[5]][2] - newco.el[f2[3]][2])/((double)order-2); ??? dw = (newco.el[f1[5]][3] - newco.el[f2[3]][3])/((double)order-2); ??? p0 = newco.el[f1[5]][0] - 1.0 * dx; ??? p1 = newco.el[f1[5]][1] - 1.0 * dy; ??? p2 = newco.el[f1[5]][2] - 1.0 * dz; ??? p3 = newco.el[f1[5]][3] - 1.0 * dw; with ( double q2 = (r4_t){p0,p1,p2,p3}; double r0 = newco.el[f1[5]][0] - 2.0 * dx; double r1 = newco.el[f1[5]][1] - 2.0 * dy; double r2 = newco.el[f1[5]][2] - 2.0 * dz; double r3 = newco.el[f1[5]][3] - 2.0 * dw; double q3 = (r4_t){r0,r1,r2,r3}; double dx = (newco.el[f1[4]][0] - newco.el[f2[4]][0])/((double)order-3); double dy = (newco.el[f1[4]][1] - newco.el[f2[4]][1])/((double)order-3); double dz = (newco.el[f1[4]][2] - newco.el[f2[4]][2])/((double)order-3); double dw = (newco.el[f1[4]][3] - newco.el[f2[4]][3])/((double)order-3); double t0 = newco.el[f1[4]][0] - 1.0 * dx; double t1 = newco.el[f1[4]][1] - 1.0 * dy; double t2 = newco.el[f1[4]][2] - 1.0 * dz; double t3 = newco.el[f1[4]][3] - 1.0 * dw; double q4 = (r4_t){t0,t1,t2,t3} ){ /* internal nodes: q2, q3, q4 */ if ((PresentNode(inter[i][1]))) { assert(! inter[i][1]->exists); ??? n = NewNumVer(inter[i][1]); { newco.el[n] = q4; } } if ((PresentNode(inter[i][2]))) { assert(! inter[i][2]->exists); ??? n = NewNumVer(inter[i][2]); { newco.el[n] = q3; } } if ((PresentNode(inter[i][3]))) { assert(! inter[i][3]->exists); ??? n = NewNumVer(inter[i][3]); { newco.el[n] = q2; } } if ((PresentNode(mocta[i][7]))) { if (! net){ assert(! mocta[i][7]->exists); } ??? n = NewNumVer(mocta[i][7]); ??? b0 = Add(Add(newco.el[f2[5]]; with (newco.el[f2[4]]),newco.el[e5[2]]), double b1 = Add(Add(newco.el[f3[5]],newco.el[f3[4]]),q4); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][8]))) { if (! net){ assert(! mocta[i][8]->exists); } ??? n = NewNumVer(mocta[i][8]); ??? b0 = Add(Add(newco.el[f2[4]]; with (newco.el[f2[3]]),newco.el[e5[1]]), double b1 = Add(Add(newco.el[f3[4]],newco.el[f3[3]]),q3); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][3]))) { if (! net){ assert(! mocta[i][3]->exists); } ??? n = NewNumVer(mocta[i][3]); ??? b0 = Add(Add(newco.el[e3[2]]; with (newco.el[f2[2]]),newco.el[f2[5]]), double b1 = Add(Add(newco.el[f1[1]],newco.el[f1[3]]),q4); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][4]))) { if (! net){ assert(! mocta[i][4]->exists); } ??? n = NewNumVer(mocta[i][4]); ??? b0 = Add(Add(newco.el[f2[2]]; with (newco.el[f2[1]]),newco.el[f2[4]]), double b1 = Add(Add(q1,q3),q4); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][5]))) { if (! net){ assert(! mocta[i][5]->exists); } ??? n = NewNumVer(mocta[i][5]); ??? b0 = Add(Add(newco.el[f2[1]]; with (newco.el[e1[2]]),newco.el[f2[3]]), double b1 = Add(Add(newco.el[f0[1]],newco.el[f0[3]]),q3); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][12]))) { if (! net){ assert(! mocta[i][12]->exists); } ??? n = NewNumVer(mocta[i][12]); ??? b0 = Add(Add(newco.el[f3[0]]; with (newco.el[f3[1]]),newco.el[e2[1]]), double b1 = Add(Add(newco.el[f0[4]],newco.el[f0[5]]),q2); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][16]))) { if (! net){ assert(! mocta[i][16]->exists); } ??? n = NewNumVer(mocta[i][16]); ??? b0 = Add(Add(newco.el[e4[1]]; with (newco.el[f3[0]]),newco.el[f3[2]]), double b1 = Add(Add(newco.el[f1[5]],newco.el[f1[4]]),q2); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][18]))) { if (! net){ assert(! mocta[i][18]->exists); } ??? n = NewNumVer(mocta[i][18]); ??? b0 = Add(Add(newco.el[e4[2]]; with (newco.el[f3[2]]),newco.el[f3[5]]), double b1 = Add(Add(newco.el[f1[4]],newco.el[f1[3]]),q4); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][19]))) { if (! net){ assert(! mocta[i][19]->exists); } ??? n = NewNumVer(mocta[i][19]); ??? b0 = Add(Add(newco.el[f3[2]]; with (newco.el[f3[1]]),newco.el[f3[4]]), double b1 = Add(Add(q2,q3),q4); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][15]))) { if (! net){ assert(! mocta[i][15]->exists); } ??? n = NewNumVer(mocta[i][15]); ??? b0 = Add(Add(newco.el[f3[1]]; with (newco.el[e2[2]]),newco.el[f3[3]]), double b1 = Add(Add(newco.el[f0[3]],newco.el[f0[4]]),q3); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][11]))) { if (! net){ assert(! mocta[i][11]->exists); } ??? n = NewNumVer(mocta[i][11]); ??? b0 = Add(Add(newco.el[f0[2]]; with (newco.el[f0[5]]),q2), double b1 = Add(Add(newco.el[e0[2]],newco.el[f1[2]]), newco.el[f1[5]]); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][14]))) { if (! net){ assert(! mocta[i][14]->exists); } ??? n = NewNumVer(mocta[i][14]); ??? b0 = Add(Add(newco.el[f0[1]]; with (newco.el[f0[2]]),newco.el[f0[4]]), double b1 = Add(Add(q1,q2),q3); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } if ((PresentNode(mocta[i][17]))) { if (! net){ assert(! mocta[i][17]->exists); } ??? n = NewNumVer(mocta[i][17]); ??? b0 = Add(Add(newco.el[f1[2]]; with (newco.el[f1[1]]),newco.el[f1[4]]), double b1 = Add(Add(q1,q2), q4); b == Add(b0,b1) ){ newco.el[n] = Scale(1.0/((double)6), b); } } } } } } return newco; } /* END SettingCoords */ PROCEDURE GlueTetrahedra(Place_t @p, Place_t b) /* This procedure realizes the gluing of two single tetrahedra. */ { if ((PnegP(Spin(a)) == NULL) && (PnegP(b)!=NULL)){ EVAL Triangulation.Glue(Spin(a),b, 1, TRUE); } else if ((PnegP(a) == NULL) && (PnegP(b) == NULL)){ EVAL Triangulation.Glue(a,Spin(b), 1, TRUE); } else { fprintf(stderr,"BuildRefinement: Not gluing some triangular walls\n"); Process.Exit¦(1); } } /* END GlueTetrahedra */ PROCEDURE GlueMacroTetrahedra(ar,br: REF Place_vec_t) /* This procedure realizes the gluing of two macro-tetrahedra trough the triangulated macro-walls. */ { for (i = 0 TO ((ar^).nel)-1){ GlueTetrahedra(ar[i],br[i]); } } /* END GlueMacroTetrahedra */ Options_t GetOptions () { Options_t *o = (Options_t *)malloc(sizeof(Options_t)); argparser_t *pp = argparser_new(stderr, argc, argv); argparser_set_help(pp, PROG_NAME " version " PROG_VERS ", usage:\n" PROG_HELP); argparser_set_info(pp, PROG_INFO); argparser_process_help_info_options(pp); { ¦TRY argparser_get_keyword(pp, "-inFileTp"); o->inFileTp = argparser_get_next(pp); argparser_get_keyword(pp, "-inFileSt"); o->inFileSt = argparser_get_next(pp); argparser_get_keyword(pp, "-outFile"); o->outFile = argparser_get_next(pp); argparser_get_keyword(pp, "-order"); o->order = argparser_get_next_int(pp, 1, 5); o->detail = argparser_keyword_present(pp, "-detail"); o->fixed = argparser_keyword_present(pp, "-fixed"); o->net = argparser_keyword_present(pp, "-net"); argparser_finish(pp); ----------------------------------- #define _HELP \ fprintf(stderr, "Usage: BuildRefinement \\\n" \ " -inFileTp \\\n" \ " -inFileSt \\\n"); fprintf(stderr, " -outFile \\\n" \ " -order [ -detail ] \\\n" \ " [ -fixed ] [ -net ]\n"); END¦ } } return o; } /* END GetOptions */ PROCEDURE ModSetGrade( uint ord; rco: Refine.Corner; bool_t Right,Left,Front,Back; ) /* Set @{edge->?}s and spheres underling on the boundary of a tetrahedron as thin cylinders and small spheres. */ void ResetNodeonNet(Place_t @p) { Node_t vn = OrgV(a); { vn->exists = FALSE; } } ResetNodeonNet; PROCEDURE Reset@{Edge->?}OnNet(Place_t @p) == { Edge_t en = PEdge(a); { en->exists = FALSE; } } Reset@{Edge->?}OnNet; void SetGrade1() { /* nothing */ } SetGrade1; void SetGrade2() { /* 3 @{edge->?}s and 0 nodes for each wall of the refined tetrahedron. */ /* On the right side. */ if (! Right) { Reset@{Edge->?}OnNet(PrevE(rco.right[0])); Reset@{Edge->?}OnNet( rco.right[3]); Reset@{Edge->?}OnNet(NextE (rco.right[2])); } /* On the left side. */ if (! Left) { Reset@{Edge->?}OnNet(PrevE(rco.left[0])); Reset@{Edge->?}OnNet( rco.left[2]); Reset@{Edge->?}OnNet(NextE (rco.left[1])); } /* On the front side. */ if (! Front) { Reset@{Edge->?}OnNet( rco.front[0]); Reset@{Edge->?}OnNet(NextE (rco.front[1])); Reset@{Edge->?}OnNet(NextE (rco.front[2])); } /* On the back side. */ if (! Back) { Reset@{Edge->?}OnNet( rco.back [0]); Reset@{Edge->?}OnNet(NextE (rco.back [1])); Reset@{Edge->?}OnNet(PrevE(rco.back [2])); } } SetGrade2; void SetGrade3() { /* 9 @{edge->?}s and 1 nodes for each wall of the refined tetrahedron. */ SetGrade2(); /* On the right side. */ if (! Right) { Reset@{Edge->?}OnNet(PrevE(rco.right[3])); Reset@{Edge->?}OnNet(PrevE(rco.right[1])); Reset@{Edge->?}OnNet( rco.right[8]); Reset@{Edge->?}OnNet(NextE (rco.right[7])); Reset@{Edge->?}OnNet( rco.right[6]); Reset@{Edge->?}OnNet(NextE (rco.right[5])); ResetNodeonNet(rco.right[3]); } /* On the left side. */ if (! Left) { Reset@{Edge->?}OnNet(PrevE(rco.left[1])); Reset@{Edge->?}OnNet(PrevE(rco.left[3])); Reset@{Edge->?}OnNet(NextE (rco.left[4])); Reset@{Edge->?}OnNet( rco.left[5]); Reset@{Edge->?}OnNet(NextE (rco.left[6])); Reset@{Edge->?}OnNet( rco.left[7]); ResetNodeonNet(rco.left[3]); } /* On the front side. */ if (! Front) { Reset@{Edge->?}OnNet( rco.front[1]); Reset@{Edge->?}OnNet( rco.front[3]); Reset@{Edge->?}OnNet(NextE (rco.front[4])); Reset@{Edge->?}OnNet(NextE (rco.front[5])); Reset@{Edge->?}OnNet(NextE (rco.front[6])); Reset@{Edge->?}OnNet(NextE (rco.front[7])); ResetNodeonNet(rco.front[3]); } /* On the back side. */ if (! Back) { Reset@{Edge->?}OnNet( rco.back [1]); Reset@{Edge->?}OnNet( rco.back [3]); Reset@{Edge->?}OnNet(NextE (rco.back [4])); Reset@{Edge->?}OnNet(PrevE(rco.back [5])); Reset@{Edge->?}OnNet(NextE (rco.back [6])); Reset@{Edge->?}OnNet(PrevE(rco.back [7])); ResetNodeonNet(rco.back[3]); } } SetGrade3; void SetGrade4() { /* 18 @{edge->?}s and 3 nodes for each wall of the refined tetrahedron. */ SetGrade3(); /* On the right side. */ if (! Right) { Reset@{Edge->?}OnNet(PrevE(rco.right[8])); Reset@{Edge->?}OnNet(PrevE(rco.right[6])); Reset@{Edge->?}OnNet(PrevE(rco.right[4])); Reset@{Edge->?}OnNet( rco.right[15]); Reset@{Edge->?}OnNet(NextE (rco.right[14])); Reset@{Edge->?}OnNet( rco.right[13]); Reset@{Edge->?}OnNet(NextE (rco.right[12])); Reset@{Edge->?}OnNet( rco.right[11]); Reset@{Edge->?}OnNet(NextE (rco.right[10])); ResetNodeonNet(rco.right[8]); ResetNodeonNet(rco.right[6]); } /* On the left side. */ if (! Left) { Reset@{Edge->?}OnNet(PrevE(rco.left[4])); Reset@{Edge->?}OnNet(PrevE(rco.left[6])); Reset@{Edge->?}OnNet(PrevE(rco.left[8])); Reset@{Edge->?}OnNet(NextE (rco.left[9 ])); Reset@{Edge->?}OnNet( rco.left[10]); Reset@{Edge->?}OnNet(NextE (rco.left[11])); Reset@{Edge->?}OnNet( rco.left[12]); Reset@{Edge->?}OnNet(NextE (rco.left[13])); Reset@{Edge->?}OnNet( rco.left[14]); ResetNodeonNet(rco.left[8]); ResetNodeonNet(rco.left[6]); } /* On the front side. */ if (! Front) { Reset@{Edge->?}OnNet( rco.front[4]); Reset@{Edge->?}OnNet( rco.front[6]); Reset@{Edge->?}OnNet( rco.front[8]); Reset@{Edge->?}OnNet(NextE (rco.front[ 9])); Reset@{Edge->?}OnNet(NextE (rco.front[10])); Reset@{Edge->?}OnNet(NextE (rco.front[11])); Reset@{Edge->?}OnNet(NextE (rco.front[12])); Reset@{Edge->?}OnNet(NextE (rco.front[13])); Reset@{Edge->?}OnNet(NextE (rco.front[14])); ResetNodeonNet(rco.front[8]); ResetNodeonNet(rco.front[6]); } /* On the back side. */ if (! Back) { Reset@{Edge->?}OnNet( rco.back [4]); Reset@{Edge->?}OnNet( rco.back [6]); Reset@{Edge->?}OnNet( rco.back [8]); Reset@{Edge->?}OnNet(NextE (rco.back [ 9])); Reset@{Edge->?}OnNet(PrevE(rco.back [10])); Reset@{Edge->?}OnNet(NextE (rco.back [11])); Reset@{Edge->?}OnNet(PrevE(rco.back [12])); Reset@{Edge->?}OnNet(NextE (rco.back [13])); Reset@{Edge->?}OnNet(PrevE(rco.back [14])); ResetNodeonNet(rco.back[8]); ResetNodeonNet(rco.back[6]); } } SetGrade4; void SetGrade5() { /* 10 @{edge->?}s and 6 nodes for each wall of the refined tetrahedron. */ SetGrade4(); /* On the right side. */ if (! Right) { Reset@{Edge->?}OnNet(PrevE(rco.right[15])); Reset@{Edge->?}OnNet(PrevE(rco.right[13])); Reset@{Edge->?}OnNet(PrevE(rco.right[11])); Reset@{Edge->?}OnNet(PrevE(rco.right[ 9])); Reset@{Edge->?}OnNet( rco.right[24]); Reset@{Edge->?}OnNet(NextE (rco.right[23])); Reset@{Edge->?}OnNet( rco.right[22]); Reset@{Edge->?}OnNet(NextE (rco.right[21])); Reset@{Edge->?}OnNet( rco.right[20]); Reset@{Edge->?}OnNet(NextE (rco.right[19])); Reset@{Edge->?}OnNet( rco.right[18]); Reset@{Edge->?}OnNet(NextE (rco.right[17])); ResetNodeonNet(rco.right[15]); ResetNodeonNet(rco.right[13]); ResetNodeonNet(rco.right[11]); } /* On the left side. */ if (! Left) { Reset@{Edge->?}OnNet(PrevE(rco.left[ 9])); Reset@{Edge->?}OnNet(PrevE(rco.left[11])); Reset@{Edge->?}OnNet(PrevE(rco.left[13])); Reset@{Edge->?}OnNet(PrevE(rco.left[15])); Reset@{Edge->?}OnNet(NextE (rco.left[16])); Reset@{Edge->?}OnNet( rco.left[17]); Reset@{Edge->?}OnNet(NextE (rco.left[18])); Reset@{Edge->?}OnNet( rco.left[19]); Reset@{Edge->?}OnNet(NextE (rco.left[20])); Reset@{Edge->?}OnNet( rco.left[21]); Reset@{Edge->?}OnNet(NextE (rco.left[22])); Reset@{Edge->?}OnNet( rco.left[23]); ResetNodeonNet(rco.left[15]); ResetNodeonNet(rco.left[13]); ResetNodeonNet(rco.left[11]); } /* On the front side. */ if (! Front) { Reset@{Edge->?}OnNet( rco.front[ 9]); Reset@{Edge->?}OnNet( rco.front[11]); Reset@{Edge->?}OnNet( rco.front[13]); Reset@{Edge->?}OnNet( rco.front[15]); Reset@{Edge->?}OnNet(NextE (rco.front[16])); Reset@{Edge->?}OnNet(NextE (rco.front[17])); Reset@{Edge->?}OnNet(NextE (rco.front[18])); Reset@{Edge->?}OnNet(NextE (rco.front[19])); Reset@{Edge->?}OnNet(NextE (rco.front[20])); Reset@{Edge->?}OnNet(NextE (rco.front[21])); Reset@{Edge->?}OnNet(NextE (rco.front[22])); Reset@{Edge->?}OnNet(NextE (rco.front[23])); ResetNodeonNet(rco.front[15]); ResetNodeonNet(rco.front[13]); ResetNodeonNet(rco.front[11]); } /* On the back side. */ if (! Back) { Reset@{Edge->?}OnNet( rco.back [ 9]); Reset@{Edge->?}OnNet( rco.back [11]); Reset@{Edge->?}OnNet( rco.back [13]); Reset@{Edge->?}OnNet( rco.back [15]); Reset@{Edge->?}OnNet(NextE (rco.back [16])); Reset@{Edge->?}OnNet(PrevE(rco.back [17])); Reset@{Edge->?}OnNet(NextE (rco.back [18])); Reset@{Edge->?}OnNet(PrevE(rco.back [19])); Reset@{Edge->?}OnNet(NextE (rco.back [20])); Reset@{Edge->?}OnNet(PrevE(rco.back [21])); Reset@{Edge->?}OnNet(NextE (rco.back [22])); Reset@{Edge->?}OnNet(PrevE(rco.back [23])); ResetNodeonNet(rco.back[15]); ResetNodeonNet(rco.back[13]); ResetNodeonNet(rco.back[11]); } } SetGrade5; { if ( ord == 1){ SetGrade1(); } else if (ord == 2){ SetGrade2(); } else if (ord == 3){ SetGrade3(); } else if (ord == 4){ SetGrade4(); } else if (ord == 5){ SetGrade5(); } else { fprintf(stderr,"Order must be less equal to 5\n"); assert(ord <= 5); } } /* END ModSetGrade */ { DoIt() } BuildRefinement. /* Copyright © 2001 Universidade Estadual de Campinas (UNICAMP) */