#define PROG_NAME "JSRefineTriang" #define PROG_DESC "???" #define PROG_VERS "1.0" #define JSRefineTriang_C_COPYRIGHT \ "" #define PROG_INFO \ "" \ " " /* Refines (or duplicate) a given triangulation (".tp" file). See notice of copyright at the end of this file. Revisions: 21-05-2000: Fixed and removed a bug that increases the number of original nodes of type \"VV\". 31-05-2000: Added the heredity of the "root" attributes for edges and walls. 11-06-2000: Hide the new elements insert in the refiment process: @{edge->?}s and nodes are set with the attibute "exists== FALSE" 25-11-2000: Added the option "net" for simulate textures on existing walls with thin cylinders and small spheres. */ #include #include #include #include #include #define _GNU_SOURCE #include REF Node_vec_t x ; NNV: uint = 0; TYPE typedef struct Options_t { char *inFileTp; /* Input file name (minus ".tp" extension) */ char *inFileSt; /* Input file name (minus ".st" extension) */ char *outFile; /* Output file name prefix */ bool_t fixOri; /* TRUE to fix original nodes */ bool_t assert; /* make some strong assertions */ bool_t net; /* simulate a net as small spheres and thin cylinders. */ } double RowT = ARRAY [0..3] OF Place_t; CONST double Thin@{Edge->?}Exists = FALSE; Options_t *GetOptions(int argc, char **argv); int main(int argc, char **argv) REF ARRAY OF Place_vec_t gi ; n{ ElemTableRec_t top; 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.ReadToTaMa(o->inFileTp; with (FALSE), /* TRUE for indicate root of tetrahedra */ top == tc.top, rc == Triangulation.ReadState(o->inFileSt), c == rc^, double half = Place_vec_new(2*top->wall.nelE)^; double vnew = Node_vec_new(top->node.nel)^ ){ fprintf(stderr, "Refining from: " & o->inFileTp & ".tp\n"); gi = NEW(REF ARRAY OF Place_vec_t, top->cell.nel, 4); Node_vec_t x = Node_vec_new(top->cell.nel); Place_t Half(Place_t @p) { ??? na = PWedge(a)->num; ??? oa = OrientationBit(a); ??? sa = SpinBit(a); ??? h = half[2*na + oa]; { if (sa == 0){ return h }else{ return Spin(h); } } } Half; /* ====== Refine the Tetrahedral Cell ====== */ void RefineCell(g: RowT; in: uint) { /* To insert four walls obliques "i", "j", "k", "l" */ ??? i = MakeTriangle(FALSE); ??? j = MakeTriangle(FALSE); ??? k = MakeTriangle(FALSE); ??? l = MakeTriangle(FALSE); { /* the wall "i" will be splice with g[1], g[2], g[3] */ SetNextF(i,Clock(g[1])); SetNextF(NextE(i),NextE(g[3])); SetNextF(PrevE(i),Clock(g[2])); /* Update the component @{edge->?} */ SetRingEdgeInfo(i, PEdge(g[1])); SetRingEdgeInfo(NextE(i), PEdge(NextE(g[3]))); SetRingEdgeInfo(PrevE(i), PEdge(g[2])); /* Update the origins */ SetOrg(i, OrgV(g[2])); SetOrg(Clock(i), OrgV(g[1])); SetOrg(NextE(i), OrgV(Clock(i))); SetOrg(Clock(NextE(i)), OrgV(NextE(g[2]))); SetOrg(PrevE(i), OrgV(Clock(NextE(i)))); SetOrg(Clock(PrevE(i)), OrgV(i)); /* the wall "j" will be splice with: g[0], g[1], g[2] */ SetNextF(j,NextE(g[0])); SetNextF(NextE(j),Clock(NextE(g[1]))); SetNextF(PrevE(j),Clock(PrevE(g[2]))); /* Update the component @{edge->?} */ SetRingEdgeInfo(j, PEdge(NextE(g[0]))); SetRingEdgeInfo(NextE(j), PEdge(NextE(g[1]))); SetRingEdgeInfo(PrevE(j), PEdge(PrevE(g[2]))); /* Update the origins */ SetOrg(j, OrgV(NextE(g[0]))); SetOrg(Clock(j), OrgV(PrevE(g[0]))); SetOrg(NextE(j), OrgV(Clock(j))); SetOrg(Clock(NextE(j)), OrgV(g[2])); SetOrg(PrevE(j), OrgV(Clock(NextE(j)))); SetOrg(Clock(PrevE(j)), OrgV(j)); /* the wall "k" will be splice with: g[0], g[1], g[3]) */ SetNextF(k,Clock(PrevF(PrevE(g[0])))); SetNextF(NextE(k),NextF(PrevE(g[1]))); SetNextF(PrevE(k),Clock(PrevF(g[3]))); /* Update the component @{edge->?} */ SetRingEdgeInfo(k,PEdge(PrevE(g[0]))); SetRingEdgeInfo(NextE(k),PEdge(PrevE(g[1]))); SetRingEdgeInfo(PrevE(k),PEdge(g[3])); /* Update the origins */ SetOrg(k, OrgV(g[0])); SetOrg(Clock(k), OrgV(PrevE(g[0]))); SetOrg(NextE(k), OrgV(Clock(k))); SetOrg(Clock(NextE(k)), OrgV(g[1])); SetOrg(PrevE(k), OrgV(Clock(NextE(k)))); SetOrg(Clock(PrevE(k)), OrgV(k)); /* the wall "l" will be splice with: g[0], g[2], g[3] */ SetNextF(l, Clock(PrevF(g[0]))); SetNextF(NextE(l), Clock(PrevF(PrevE(g[3])))); SetNextF(PrevE(l), NextF(NextE(g[2]))); /* Update the component @{edge->?} */ SetRingEdgeInfo(l, PEdge(g[0])); SetRingEdgeInfo(NextE(l), PEdge(PrevE(g[3]))); SetRingEdgeInfo(PrevE(l), PEdge(NextE(g[2]))); /* Update the origins */ SetOrg(l, OrgV(Clock(g[0]))); SetOrg(Clock(l), OrgV(g[0])); SetOrg(NextE(l), OrgV(Clock(l))); SetOrg(Clock(NextE(l)), OrgV(NextE(g[2]))); SetOrg(PrevE(l), OrgV(Clock(NextE(l)))); SetOrg(Clock(PrevE(l)), OrgV(l)); /* Subdivision of the Octahedron (in four tetrahedron) delimitate by the medial nodes */ ??? o1 = MakeTriangle(FALSE); ??? o2 = MakeTriangle(FALSE); ??? o3 = MakeTriangle(FALSE); ??? o4 = MakeTriangle(FALSE); { assert(i == PrevF(Clock(g[1]))); assert(i == Clock(NextF(g[1]))); assert(i == Clock(PrevE(NextF(g[2])))); assert(j == PrevF(NextE(g[0]))); assert(k == Clock(PrevF(PrevE(g[0])))); assert(k == Clock(PrevE(PrevF(g[3])))); assert(l == Clock(PrevF(g[0]))); /* triangle o1 */ SetNextF(o1, NextE(k)); SetNextF(PrevE(g[1]), o1); SetNextF(PrevE(o1), NextE(g[0])); SetNextF(j, PrevE(o1)); /* Update the component @{edge->?} */ SetRingEdgeInfo(o1, PEdge(NextE(k))); SetRingEdgeInfo(PrevE(o1), PEdge(j)); /* Update the origins */ SetOrg(o1, OrgV(PrevE(g[1]))); SetOrg(Clock(o1), OrgV(Clock(PrevE(g[1])))); SetOrg(PrevE(o1), OrgV(NextE(g[0]))); SetOrg(Clock(PrevE(o1)), OrgV(Clock(NextE(g[0])))); SetOrg(NextE(o1), OrgV(PrevE(k))); SetOrg(Clock(NextE(o1)), OrgV(j)); /* triangle o2 */ SetNextF(NextE(o2), NextE(g[3])); SetNextF(NextE(i), NextE(o2)); SetNextF(PrevE(o2), PrevE(l)); SetNextF(NextE(g[2]), PrevE(o2)); /* Update the component @{edge->?} */ SetRingEdgeInfo(NextE(o2), PEdge(NextE(g[3]))); SetRingEdgeInfo(PrevE(o2), PEdge(NextE(g[2]))); /* Update the origins */ SetOrg(o2, OrgV(j)); SetOrg(Clock(o2), OrgV(PrevE(k))); SetOrg(NextE(o2), OrgV(NextE(g[3]))); SetOrg(Clock(NextE(o2)), OrgV(Clock(NextE(g[3])))); SetOrg(PrevE(o2), OrgV(NextE(g[2]))); SetOrg(Clock(PrevE(o2)), OrgV(Clock(NextE(g[2])))); /* triangle o3 */ SetNextF(o3, l); SetNextF(Clock(g[0]), o3); SetNextF(NextE(o3), g[3]); SetNextF(Clock(PrevE(k)), NextE(o3)); /* Update the component @{edge->?} */ SetRingEdgeInfo(o3, PEdge(l)); SetRingEdgeInfo(NextE(o3), PEdge(g[3])); /* Update the origins */ SetOrg(o3, OrgV(Clock(g[0]))); SetOrg(Clock(o3), OrgV(g[0])); SetOrg(NextE(o3), OrgV(g[3])); SetOrg(Clock(NextE(o3)), OrgV(Clock(g[3]))); SetOrg(PrevE(o3), OrgV(PrevE(k))); SetOrg(Clock(PrevE(o3)), OrgV(j)); /* triangle o4 */ SetNextF(g[1], NextE(o4)); SetNextF(NextE(o4), Clock(i)); SetNextF(PrevE(j), PrevE(o4)); SetNextF(PrevE(o4), Clock(PrevE(g[2]))); SetNextF(o4, o2); SetNextF(o2, Clock(PrevE(o3))); SetNextF(Clock(PrevE(o3)), Clock(NextE(o1))); /* Update the component @{edge->?} */ SetRingEdgeInfo(NextE(o4), PEdge(Clock(i))); SetRingEdgeInfo(PrevE(o4), PEdge(PrevE(j))); SetRingEdgeInfo(o2, PEdge(o2)); /* Update the origins */ SetOrg(NextE(o4), OrgV(g[1])); SetOrg(Clock(NextE(o4)), OrgV(Clock(g[1]))); SetOrg(PrevE(o4), OrgV(Clock(PrevE(g[2])))); SetOrg(Clock(PrevE(o4)), OrgV(PrevE(g[2]))); SetOrg(o4, OrgV(o2)); SetOrg(Clock(o4), OrgV(Clock(o2))); /* making eigth tetrahedral cells */ SetPnegOfNearbyWalls(j); SetPnegOfNearbyWalls(Clock(PrevF(PrevE(NextF(k))))); SetPnegOfNearbyWalls(Clock(PrevF(NextE(i)))); SetPnegOfNearbyWalls(Clock(PrevF(NextE(NextF(l))))); SetPnegOfNearbyWalls(PrevE(NextF(o1))); SetPnegOfNearbyWalls(o2); SetPnegOfNearbyWalls(NextF(o3)); SetPnegOfNearbyWalls(o4); if (o->assert){ /* Some strong assertions */ assert(Org(Srot(j)) == PnegP(j)); assert(OrgV(Srot(Clock(PrevF(PrevE(NextF(k)))))) == PnegP(Clock(PrevF(PrevE(NextF(k)))))); assert(OrgV(Srot(Clock(PrevF(NextE(i))))) == PnegP(Clock(PrevF(NextE(i))))); assert(OrgV(Srot(Clock(PrevF(NextE(NextF(l)))))) == PnegP(Clock(PrevF(NextE(NextF(l)))))); assert(OrgV(Srot(NextF(o3))) == PnegP(NextF(o3))); assert(OrgV(Srot(PrevE(NextF(o1)))) == PnegP(PrevE(NextF(o1)))); assert(OrgV(Srot(o4)) == PnegP(o4)); assert(OrgV(Srot(o2)) == PnegP(o2)); Place_t v = Srot(j); Place_t v1 = Clock(PrevE(Srot(j))); { assert(PnegP(Tors(v)) == PnegP(Tors(v1))); } Place_t v = Srot(NextF(o3)); Place_t v1 = Clock(PrevE(Srot(NextF(o3)))); { assert(PnegP(Tors(v)) == PnegP(Tors(v1))); } Place_t v = Srot(o2); Place_t v1 = Clock(PrevE(Srot(o2))); { assert(PnegP(Tors(v)) == PnegP(Tors(v1))); } Place_t v = Srot(o4); Place_t v1 = Clock(PrevE(Srot(o4))); { assert(PnegP(Tors(v)) == PnegP(Tors(v1))); } Place_t v = Srot(PrevE(NextF(o1))); Place_t v1 = Clock(PrevE(Srot(PrevE(NextF(o1))))); { assert(PnegP(Tors(v)) == PnegP(Tors(v1))); } Place_t v = Srot(Clock(PrevF(PrevE(NextF(k))))); Place_t v1 = Clock(PrevE(Srot(Clock(PrevF(PrevE(NextF(k))))))); { assert(PnegP(Tors(v)) == PnegP(Tors(v1))); } Place_t v = Srot(Clock(PrevF(NextE(i)))); Place_t v1 = Clock(PrevE(Srot(Clock(PrevF(NextE(i)))))); { assert(PnegP(Tors(v)) == PnegP(Tors(v1))); } Place_t v = Srot(Clock(PrevF(NextE(NextF(l))))); Place_t v1 = Clock(PrevE(Srot(Clock(PrevF(NextE(NextF(l))))))); { assert(PnegP(Tors(v)) == PnegP(Tors(v1))); } } /* Now, refine the subdivision of the octahedron in more four tetrahedra, through the subdivision of the diagonal @{edge->?} of the octahedron. */ Subdivide@{Edge->?}(o2,in,o->net); } } } RefineCell; bool_t CreateNode(Place_t @p) Place_t t = a; VAR { ??? fe = NARROW(PWedge(a)-> Wedge_t); { ){ if (! fe.mark) { fe.mark = TRUE; do { ??? fe = NARROW(PWedge(t), Wedge_t); { ){ fe.mark = TRUE; } t = NextF(t) } while (t != a); return TRUE; } return FALSE; } } CreateNode; void SetAllVh(Place_t @p, Node_t v) /* Set all adjacents @places to "a" with the same medial node "v" */ void SetVh(b: Place_t) /* Set the @place "b" with the medial node "v". */ { ??? fe = NARROW(PWedge(b)-> Wedge_t); { ){ fe.vh = v; } } SetVh; Place_t t = a; VAR { do { SetVh(t); t = NextF(t); } while (t != a); } SetAllVh; void SetPnegOfNearbyWalls(Place_t @p) /* set the (12) places belonging to same 3-cell */ Place_t t = a; VAR { Cell_t p = MakeCell(0); { SetPnegOfWall(t,p); do { SetPnegOfWall(Clock(PrevF(t)),p); t = PrevE(t); } while (t != a); } } SetPnegOfNearbyWalls; uint Vhnum(Place_t @p) /* given the @place "a", this procedure return the number of its medial node "vh". */ { ??? fe = NARROW(PWedge(a)-> Wedge_t); { ){ return fe.vh->num; } } Vhnum; { /* Copy the original nodes, save icorrespondence in "vnew" array: */ for (iu = 0; iu < top->node.nel; iu++) { ??? u = top->node[iu]; Node_t v = MakeNode(u->num); { v->exists = u->exists; v.fixed = u.fixed) || (o->fixOri; v.color = u.color; v.radius = u.radius; v.label = u.label; vnew[iu] = v } } /* Create two new wedges for each original wedge "fe". The new wedge corresponding to the origin half of "fe", with same spin and orientation, will be | Half(a) == Spin^s(half[2*PWedge(a)->num + o]) where s == SpinBit(a), o == OrientationBit(a) */ REF Node_vec_t ve = Node_vec_new(top->NE); uint i = 0; { for (ie = 0; ie < top->wall.nelE; ie++) { ??? a = top->wedge[ie]; ??? oa = OrientationBit(a); ??? sa = SpinBit(a); ??? ho = half[2*ie + oa]; ??? hd = half[2*ie + 1 - oa]; ??? fe = NARROW(PWedge(a)-> Wedge_t); { double fn = fe.wall->num; double en = fe.edge->num ){ if (( CreateNode(a))) { ve[i] = MakeNode(top->node.nel + NNV); NNV++; ve[i]->exists = top->edge[en]->exists; ve[i].fixed = FALSE; ve[i].color = top->edge[en].color; ve[i].transp = top->edge[en].transp; ve[i].radius = top->edge[en].radius; ve[i].label = "VE"; SetAllVh(a, ve[i]); i++; } ho = MakeWedge(); ??? hoe = NARROW(PWedge(ho)-> Wedge_t); { ){ hoe.edge->exists = top->edge[en]->exists; hoe.edge.color = top->edge[en].color; hoe.edge.transp = top->edge[en].transp; hoe.edge.radius = top->edge[en].radius; /* set the "root" @{edge->?} */ hoe.edge->root = top->edge[en]->root; hoe.wall->exists = top->wall[fn]->exists; hoe.wall.color = top->wall[fn].color; hoe.wall.transp = top->wall[fn].transp; } if (sa == 1){ ho = Spin(ho); } hd = MakeWedge(); ??? hde = NARROW(PWedge(hd)-> Wedge_t); { ){ hde.edge->exists = top->edge[en]->exists; hde.edge.color = top->edge[en].color; hde.edge.transp = top->edge[en].transp; hde.edge.radius = top->edge[en].radius; /* set the "root" @{edge->?} */ hde.edge->root = top->edge[en]->root; hde.wall->exists = top->wall[fn]->exists; hde.wall.color = top->wall[fn].color; hde.wall.transp = top->wall[fn].transp; } if (sa == 1){ hd = Spin(hd); } SpliceEdges(ho, Clock(hd)); ??? m = Vhnum(a); { SetOrg(Clock(ho),ve[m-top->node.nel]); SetOrg(Clock(hd),ve[m-top->node.nel]); } SetOrg(ho, vnew[OrgV(a)->num]); SetOrg(hd, vnew[OrgV(Clock(a))->num]); } } } /* Connect the half-wedges as in the original triangulation */ for (ie = 0; ie < top->wall.nelE; ie++) { ??? a = top->wedge[ie]; Place_t b = NextF(a); Place_t c = NextE(a); ??? ha = Half(a); ??? hac = Half(Clock(a)); ??? hb = Half(b); ??? hc = Half(c); { if ((b!=a) && (NextF(ha)!=hb)) { SetNextF(ha, hb); /* so, NextF(ha) == hb */ SetEdgeInfo(hb, PEdge(ha)); } if ((c!=a) && (NextE(Clock(hac))!=hc)) { SetNextE(Clock(hac), hc); /* so, NextE(Clock(hac) == hc */ } } Place_t a = Clock(top->wedge[ie]); Place_t b = NextF(a); Place_t c = NextE(a); ??? ha = Half(a); ??? hac = Half(Clock(a)); ??? hb = Half(b); ??? hc = Half(c); { if ((b!=a) && (NextF(ha)!=hb)) { SetNextF(ha, hb); /* so, NextF(ha) == hb */ SetEdgeInfo(hb, PEdge(ha)); } if ((NextE(Clock(hac))!=hc)) { SetNextE(Clock(hac), hc); /* so, NextE(Clock(hac) == hc */ } } } for (j = 0; j < top->wall.nel; j++) { Place_t *d,e,f,g; { ??? fa = top->wall[j]; ??? fr = fa->root; ??? a = fa.pa; Place_t b = NextE(a); Place_t c = PrevE(a); ??? ha = Half(a); with (hac == Half(Clock(a)), double hb = Half(b), hbc == Half(Clock(b)); double hc = Half(c), hcc == Half(Clock(c)) ){ if ((((fa.node^).nel)!=3)){ BadElement(); } d = MakeWedge(); PEdge(d)->exists = Thin@{Edge->?}Exists; e = MakeWedge(); PEdge(e)->exists = Thin@{Edge->?}Exists; f = MakeWedge(); PEdge(f)->exists = Thin@{Edge->?}Exists; if ((fa->exists) && (o->net)) { with ( double co = (frgb_t){1.00,1.00,0.50}, /* color and radius for the thin */ double ra = 0.0025, /* @{edge->?} on the net */ double de = PEdge(d); double ee = PEdge(e); double fe = PEdge(f) ){ de->exists = TRUE; ee->exists = TRUE; fe->exists = TRUE; de.color = co; de.radius = ra; ee.color = co; ee.radius = ra; fe.color = co; fe.radius = ra; } } /* Make the first @{link->?} wedge */ SetNextE(hc, d); SetNextE(d, Clock(hbc)); SetNextE(Clock(hc), hbc); /* Update the origins */ SetOrg(d, OrgV(Clock(hc))); SetOrg(Clock(d), OrgV(Clock(hbc))); /* Update the component wall */ SetWallInfo(hc, PWall(d)); SetWallInfo(hbc, PWall(d)); PWall(d)->exists = top->wall[j]->exists; PWall(d).color = top->wall[j].color; PWall(d).transp = top->wall[j].transp; /* Make the second @{link->?} wedge */ SetNextE(ha, e); SetNextE(e, Clock(hcc)); SetNextE(Clock(ha), hcc); /* Update the origins */ SetOrg(e, OrgV(Clock(ha))); SetOrg(Clock(e), OrgV(Clock(hcc))); /* Update the component wall */ SetWallInfo(ha, PWall(e)); SetWallInfo(hcc, PWall(e)); PWall(e)->exists = top->wall[j]->exists; PWall(e).color = top->wall[j].color; PWall(e).transp = top->wall[j].transp; /* Make the third @{link->?} wedge */ SetNextE(hac, f); SetNextE(f, Clock(hb)); SetNextE(Clock(hac), hb); /* Update the origins */ SetOrg(f, OrgV(Clock(hac))); SetOrg(Clock(f), OrgV(Clock(hb))); /* Update the component wall */ SetWallInfo(hac, PWall(f)); SetWallInfo(hb, PWall(f)); PWall(f)->exists = top->wall[j]->exists; PWall(f).color = top->wall[j].color; PWall(f).transp = top->wall[j].transp; /* Create Triangular wall to insert */ if (top->wall[j]->exists ){ g = MakeTriangle(TRUE); } else { g = MakeTriangle(FALSE); } /* making the connections */ SetNextF(NextE(hc), Clock(g)); SetNextF(NextE(hbc), g); SetNextF(NextE(g), NextE(hcc)); SetNextF(Clock(NextE(g)), Clock(NextE(hcc))); SetNextF(Clock(NextE(hb)), PrevE(g)); SetNextF(Clock(NextE(hac)),Clock(PrevE(g))); /* Update the edge component */ SetRingEdgeInfo(g ,PEdge(d)); SetRingEdgeInfo(NextE(g) ,PEdge(e)); SetRingEdgeInfo(PrevE(g) ,PEdge(f)); /* Update the origins */ SetOrg(g, OrgV(Clock(d))); SetOrg(Clock(g), OrgV(d)); SetOrg(NextE(g), OrgV(Clock(e))); SetOrg(Clock(NextE(g)), OrgV(e)); SetOrg(PrevE(g), OrgV(f)); SetOrg(Clock(PrevE(g)), OrgV(Clock(f))); PWall(g).color = top->wall[j].color; PWall(g).transp = top->wall[j].transp; /* set the "root" wall */ PWall(d)->root = fr; PWall(e)->root = fr; PWall(f)->root = fr; PWall(g)->root = fr; /* end of the setting */ } } } for (i = 0; i < top->wall.nelE; i++) { ??? a = top->wedge[i]; ??? ha = Half(a); ??? hac = Half(Clock(a)); { SetRingEdgeInfo(ha, PEdge(ha)); SetRingEdgeInfo(hac, PEdge(hac)); } }; for (i = 0; i < top->cell.nel; i++) { ??? v = Srot(top.cell[i]); Place_t a = Tors(v); Place_t af = PrevF(a); Place_t ae = PrevE(a); Place_t aee = NextE(a); { assert(PnegP(a)->num == i); assert(OrgV(v)->num == i); if ((PposP(af)!=NULL)){ assert(PnegP(a) == PposP (af)); } if ((NUMBER(PnegP(a).node^)!=4)){ BadElement(); } gi[i,0] = Clock(NextE(PrevF(NextE(Half(a))))); gi[i,1] = Clock(NextE(NextF(NextE(Half(af))))); gi[i,2] = Clock(NextE(NextF(NextE(Half(PrevF(ae)))))); gi[i,3] = NextF(NextE(Half(PrevF(aee)))); } } for (i = 0; i < top->cell.nel; i++) { RefineCell(RowT{gi[i,0],gi[i,1],gi[i,2],gi[i,3]},i); } if (top->bdr == 0) { ntop = MakeElemTable(Half(top->wedge[1]),0); } else { ntop = MakeElemTable(Half(top->wedge[1]),1); } ??? nc = NEW(REF Coords_t; with (ntop.node.nel)^, double com = "Refined from: " & o->inFileTp & ".tp\n" & "Created by RefineTriang: " & o->outFile & ".tp" ){ assert(ntop.node.nel == top->node.nel + top->NE + top->cell.nel); assert(ntop.edge.nel == 2*top->NE + 3*top->wall.nel + 6*top->cell.nel); assert(ntop.wall.nel == 4*top->wall.nel + 16*top->cell.nel); assert(ntop.cell.nel == 12*top->cell.nel); if (top->bdr == 0){ assert(ntop.wedge.nel == 72*top->cell.nel); } for (j = 0; j < top->wall.nelE; j++) { ??? a = top->wedge[j]; Node_t ou = OrgV(a)->num; Node_t ov = OrgV(Clock(a))->num; Node_t nu = OrgV(Half(a))->num; Node_t nv = OrgV(Half(Clock(a)))->num; Node_t nx = OrgV(Clock(Half(a)))->num; Node_t ny = OrgV(Clock(Half(Clock(a))))->num; { assert(nx == ny); nc[nu] = c[ou]; nc[nv] = c[ov]; nc[nx] = r4_Scale(0.5, r4_Add(c[ou], c[ov])); } } if (o->net) { /* we compute the number of existing walls in the previus topology. */ VAR nefp: uint = 0; /* existing walls in the previus topology*/ neep: uint = 0; /* existing @{edge->?}s in the previus topology*/ neea: uint = 0; /* existing @{edge->?}s in the actual topology*/ { for (i = 0; i < top->wall.nel; i++) { ??? f = top->wall[i]; { if (f->exists){ nefp++;} } } for (i = 0; i < top->NE; i++) { ??? e = top->edge[i]; { if (e->exists){ neep++; } } } for (i = 0; i < ntop.edge.nel; i++) { ??? e = ntop.edge[i]; { if (e->exists){ neea++;} } } assert(neea == 3 * nefp + 2 * neep); } } for (i = 0; i < top->cell.nel; i++) { Node_t ou = OrgV(gi[i; with (1])->num, double ov = OrgV(PrevE(gi[i,2]))->num ){ nc[x[i]->num] = r4_Scale(0.5, r4_Add(nc[ou], nc[ov])); } } WriteTopology (o->outFile, ntop, com); WriteTable (o->outFile, ntop, com); WriteState (o->outFile, ntop, nc, com); WriteMaterials(o->outFile, ntop, com, FALSE); /* Now, unmark the attribute "mark" of Wedge_t */ for (i = 0; i < top->wall.nelE; i++) { ??? fe = NARROW(PWedge(top->wedge[i]), Wedge_t); { ){ if (fe.mark){ fe.mark = FALSE; } } } } } return 0; } PROCEDURE BadElement() { fprintf(stderr,"This topology isn't a triangulation\n"); Process.Exit¦(1); } /* END BadElement */ PROCEDURE Subdivide@{Edge->?}(an : Place_t; i: INTEGER; bool_t net) == VAR a,bn,m1,m2,m3,t0,t1,Place_vec_t t2; wn, p: REF ARRAY OF @{Node->?}; CONST double n = 4; { a = NEW(REF Place_vec_t,n); bn = NEW(REF Place_vec_t,n); m1 = NEW(REF Place_vec_t,n); m2 = NEW(REF Place_vec_t,n); m3 = NEW(REF Place_vec_t,n); t0 = NEW(REF Place_vec_t,n); t1 = NEW(REF Place_vec_t,n); t2 = NEW(REF Place_vec_t,n); /* save the @places */ a = NEW(REF Place_vec_t,n); a[0] = an; for (i = 1; i < n; i++) { a[i] = NextF(a[i-1]); } /* save the nodes */ wn = NEW(REF ARRAY OF @{Node->?},n); for (i = 0; i < n; i++) { wn[i] = OrgV(PrevE(a[i])); } /* save the cells */ p = NEW(REF ARRAY OF @{Node->?},n); for (i = 0; i < n; i++) { p[i] = PnegP(a[i]); } /* save other @places */ bn = NEW(REF Place_vec_t,n); for (i = 0; i < n; i++) { bn[i] = Clock(PrevE(NextF(PrevE(a[i])))); } /* insert wedges and edges */ for (i = 0; i < n; i++) { m1[i] = MakeWedge(); m2[i] = MakeWedge(); m3[i] = MakeWedge(); t0[i] = MakeTriangle(FALSE); t1[i] = NextE(t0[i]); t2[i] = NextE(t1[i]); /* If net==TRUE the simulates the net with cylinders ans spheres, */ if ((net) && (PEdge(bn[i])->exists)) { PEdge(t0[i])->exists = TRUE; PEdge(t0[i]).radius = 0.0025; PEdge(t0[i]).color = (frgb_t){1.0,1.0,0.5}; } ??? f = PWedge(t0[i]).wall; ??? e1 = PEdge(t1[i]); ??? e2 = PEdge(t2[i]); { f->exists = FALSE; e1->exists = FALSE; e2->exists = FALSE; } } /* Now subdivide @{edge->?} and extend the subdivision on the star of the edge */ x[i] = MakeNode(0); ??? v = NARROW(x[i], Node); { ){ v->exists = FALSE; v.label = "VE"; } for (j = 0; j < n; j++) { ??? ee = PEdge(a[j]); ??? ff = PEdge(m2[j]); Place_t b = NextE(a[j]); Edge_t be = PEdge(b); Place_t c = NextE(b); Edge_t ce = PEdge(c); ??? u = OrgV(a[j]); ??? v = OrgV(b); ??? w = OrgV(c); with ( /* save the attributes of the edge-wall component of the @place a[j] */ f = PWedge(a[j]).wall, g = PWedge(m3[j]).wall, double ge = g->exists; h = PEdge(m3[j]) ){ ee->exists = FALSE; ff->exists = FALSE; SetNextE(a[j],m1[j]); SetNextE(m1[j],c); SetNextE(m2[j],m3[j]); SetNextE(m3[j],Clock(b)); SetOrg(a[j], u); SetOrg(Clock(a[j]), x[i]); SetOrg(m2[j],v); SetOrg(Clock(m2[j]), x[i]); SetOrg(m3[j], x[i]); SetOrg(Clock(m3[j]), w); SetOrg(m1[j], x[i]); SetOrg(Clock(m1[j]), w); SetNextF(m1[j],m3[j]); /* set the attributes for the wall component */ ge = FALSE; SetRingWallInfo(a[j],f); SetRingWallInfo(m3[j],g); SetRingEdgeInfo(m3[j],h); SetRingEdgeInfo(b,be); SetRingEdgeInfo(c,ce); SetRingEdgeInfo(a[j],ee); SetRingEdgeInfo(m2[j],ff); SetRingWallInfo(bn[j],PWedge(bn[j]).wall); SetRingWallInfo(PrevF(bn[j]),PWedge(PrevF(bn[j])).wall); } } for (j = 0; j < n; j++) { SetNextF(Clock(m2[j]),Clock(m2[(j+1) % n])); } ??? ff = PEdge(m2[0]); { SetRingEdgeInfo(m2[0],ff); } for (j = 0; j < n; j++) { Place_t cn = NextF(bn[j]); ??? e0 = PEdge(t0[j]); ??? e1 = PEdge(t1[j]); ??? e2 = PEdge(t2[j]); { SetOrgAll(t0[j],wn[j]); SetOrgAll(t1[j],wn[(j+1) % n]); SetOrgAll(t2[j],x[i]); SetNextF(bn[j], t0[j]); SetNextF(t0[j],cn); SetNextF(m1[j],t2[j]); SetNextF(t2[j],m3[j]); SetNextF(Clock(m1[(j+1) % n]), t1[j]); SetNextF(t1[j],Clock(m3[(j+1) % n])); SetRingEdgeInfo(t0[j],e0); SetRingEdgeInfo(t1[j],e1); SetRingEdgeInfo(t2[j],e2); SetRingEdgeInfo(NextE(t0[j]),PEdge(NextE(t0[j]))); SetRingEdgeInfo(NextE(t1[j]),PEdge(NextE(t1[j]))); SetRingEdgeInfo(NextE(t2[j]),PEdge(NextE(t2[j]))); } } /* insert cells */ for (j = 0; j < n; j++) { ??? q = Triangulation.MakeCell(0); { Triangulation.SetPnegOfNearbyWalls(a[j],p[j]); Triangulation.SetPnegOfNearbyWalls(m2[j],q); } } } /* END Subdivide@{Edge->?} */ 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); o->fixOri = argparser_keyword_present(pp, "-fixOri"); o->net = argparser_keyword_present(pp, "-net"); argparser_finish(pp); ----------------------------------- #define _HELP \ fprintf(stderr, "Usage: RefineTriang \\\n"); fprintf(stderr, " -inFileTp \\\n"); fprintf(stderr, " -inFileSt \\\n"); fprintf(stderr, " -outFile \\\n"); fprintf(stderr, " [ -fixOri ] [-assert] [-net]\n"); END¦ } } return o; } GetOptions; Place_t MakeTriangle(bool_t exists) /* Make one triangular wall and set of the three places with the same wall component. If exists is TRUE then the triangular wall exists, FALSE otherwise. */ { Place_t a = MakeWedge(); Place_t b = MakeWedge(); Place_t c = MakeWedge(); Wall_t f = PWall(a); Node_t u = MakeNode(0); Node_t v = MakeNode(1); Node_t w = MakeNode(2); { if (exists){ f->exists = TRUE }else{ f->exists = FALSE; } PEdge(a)->exists = Thin@{Edge->?}Exists; SetOrg(a, u); SetOrg(Clock(a),v); PEdge(b)->exists = Thin@{Edge->?}Exists; SetNextE(a,b); SetWallInfo (b,f); SetOrg (b,v); SetOrg(Clock(b),w); PEdge(c)->exists = Thin@{Edge->?}Exists; SetNextE(b,c); SetWallInfo (c,f); SetOrg (c, w); SetOrg(Clock(c), OrgV(a)); return a; } } /* END MakeTriangle */ { DoIt() } RefineTriang. /* Copyright © 2000 Universidade Estadual de Campinas (UNICAMP) */