// #FILE timestamp.txt /* Last edited on 2007-01-31 14:25:24 by stolfi */ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // #FILE cur/libenergy/JUNK/OldCurvature2D.c /* See {OldCurvature2D.h} */ // #INCLUDE [!!!JUNK!!!] // #INCLUDE // #INCLUDE #include // #INCLUDE FROM Triangulation IMPORT Topology, Face, OrgV; FROM Energy IMPORT Coords, Gradient; FROM LR4 IMPORT Add, Neg, Dot; FROM Octf IMPORT Enext, Enext_1; CONST double zero = (r4_t){0.0,0.0,0.0,0.0}; double IniStackSize = 100000; double Epsilon = 0.0000000001; VAR str,stc: Stat.T; /* statistical accumulators to the number of "root" elements and the number of "children" elements inside each "root" element. */ TYPE double bool_vec_t = ARRAY OF bool_t; double StackF = REF Face_vec_t; double Number = RECORD INTEGER nre ; /* number of "root" faces. */ uint nce ; /* number of "children" faces inside each "root" face.*/ ; } REVEAL double T = Public BRANDED OBJECT double K ; /* The energy normalization factor */ top : Topology; /* The topology */ bool_vec_t vVar; /* TRUE if vertex is variable */ Number num ; /* Number of "root" and "children" faces*/ ChilFace: REF ARRAY OF StackF; /* The "children" faces */ OVERRIDES init = Init; defTop = DefTop; defVar = DefVar; eval = Eval; name = Name; ; } PROCEDURE Statistics(Topology *top) : Number == Number *num; { for (i = 0; i < top->NF; i++){ ??? f = top->face[i]; ??? fr = FLOAT(f.root, REAL), ); with( ){ Stat.Accum(str,fr); if ((fr == 0.0)){ Stat.Accum(stc,fr);}; ;} ; } num.nre = FLOOR(str.maximum)+1; num.nce = FLOOR(stc.num); return num; } /* END Statistics */ PROCEDURE CropChilFaces( *top: Triangulation.Topology; Number *num; ) : REF ARRAY OF StackF == VAR topi : REF uint_vec_t ; /* Crop the "children" faces for each "root" face. */ { /* initialize the "top" indexes for each of the "num.nre" stacks of faces. */ topi = NEW(REF uint_vec_t , num.nre); for (k = 0; k <= (num.nre-1); k++){ topi[k] = 0; } ??? t = NEW(REF ARRAY OF StackF, num.nre); with ( ){ for (k = 0; k <= (num.nre-1); k++){ t[k] = Face_vec_new(IniStackSize); ; } for (j = 0; j < top->NF; j++){ ??? f = top->face[j]; with ( double fr = f.root DO if ((fr!=-1)){ SaveF(t[fr],topi[fr],f);}; ;} ; } return t; ; } } /* END CropChilFaces */ void SaveF( StackF *Stack ; uint *top; Face *face ; ) == /* Save the face "face" on the stack "Stack" */ { Stack[top] = face; top = top +1 } /* END SaveF */ SELF_Init(SELF_T *erg) { return erg } /* END Init */ void SELF_DefTop(SELF_T *erg, Topology_t *top) { erg->K = 1.0; erg->top = top; erg->vVar = bool_vec_new(top->NV); erg->num = Statistics(top); erg->ChilFace = CropChilFaces(top,erg->num); /* Just in case the client forgets to call "defVar": */ for (i = 0; i < top->NV; i++){ erg->vVar.el[i] = FALSE; } } /* END DefTop */ void SELF_DefVar(SELF_T *erg, bool_vec_t *variable) { int NV = erg->top->NV; ??? vVar = erg->vVar^; { assert(((variable).nel) == NV); vVar = variable; ;} } /* END DefVar */ void SELF_Eval( SELF_T erg; Coords *c; double *e; bool_t grad; Gradient *eDc; ) { ??? NV = erg->top->NV; ??? K = erg->K; ??? vVar = erg->vVar^; ??? ChilFace = erg->ChilFace; ??? num = erg->num; { PROCEDURE AddTerm(*iu,iv,iw,ix: uint) == /* Adds one term of the curvature energy to "e" (and its derivative to "eDc, if "grad" is TRUE). The terms correspond to the faces "f1 == u v w" and "f2 == u w x". See the follow picture: v /|\ / | \ / | \ / / | \ \ ------ w / | \ x ----- \ \ f1 | f2 / / \ | / \ | / \ | / \|/ u */ VAR eterm: double; eDdu, eDdv, eDdw, eDdx: r4_t; { ??? u = c[iu]; ??? v = c[iv]; ??? w = c[iw]; ??? x = c[ix]; { term(u,v,w,x, eterm, eDdu, eDdv, eDdw, eDdx); e = e + eterm; if ((grad)){ if ((vVar.el[iu])){ eDc[iu] = r4_Add(eDc[iu],eDdu); ; } if ((vVar.el[iv])){ eDc[iv] = r4_Add(eDc[iv],eDdv); ; } if ((vVar.el[iw])){ eDc[iw] = r4_Add(eDc[iw],eDdw); ; } if ((vVar.el[ix])){ eDc[ix] = r4_Add(eDc[ix],eDdx); ;} ;} ;} ;} AddTerm; PROCEDURE term(u,v,w,x: r4_t; double *eterm; VAR dedu,dedv,dedw,dedx: r4_t; ) == VAR dedDv,dedDw,dedDx: r4_t; { r4_t Dv = r4_Sub(v, u); /* V */ with ( r4_t Dw = r4_Sub(w, u); /* V */ r4_t Dx = r4_Sub(x, u); ){ EangleAux(Dv,Dw, Dx, eterm, dedDv,dedDw,dedDx); dedv = dedDv; /* V */ dedw = dedDw; /* V */ dedx = dedDx; /* V */ dedu = Neg(Add(Add(dedDv,dedDw),dedDx)); ;} ;} term; PROCEDURE EangleAux(f,a,b: r4_t; double *eterm; VAR dedf,deda,dedb: r4_t; ) == /* compute the derivative of the orthogonal vectors "f" and "g" where "f== s - Proj(s,r)" and "g == r". */ VAR dedr,deds: r4_t; { ??? m = r4_Dot(f,f)+Epsilon; /* S */ with ( u == r4_Dot(f,a), /* S */ v == r4_Dot(f,b), /* S */ U == u/m, /* S */ V == v/m, /* S */ double Uf = r4_Scale(U,f), /* V */ double Vf = r4_Scale(V,f), /* V */ double R = r4_Sub(a,Uf), /* V */ double S = r4_Sub(b,Vf) /* V */ ){ Eangle(R,S,eterm,dedr,deds); ??? dedV = - Dot(deds,f); /* S */ with ( dedU == - Dot(dedr, f), /* S */ dedu == dedU/m, /* S */ dedv == dedV/m, /* S */ dedm == (-1.0/m) * (dedU*U + dedV*V), /* S */ dedm2 == r4_Scale(2.0,f), /* S */ double dedm2f = r4_Scale(dedm, dedm2), /* V */ dedua == r4_Scale(dedu, a), /* V */ dedvb == r4_Scale(dedv, b), /* V */ dedrU == r4_Neg(r4_Scale(U, dedr)), /* V */ dedsV == r4_Neg(r4_Scale(V, deds)), /* V */ t1 == r4_Add(dedm2f,dedua), /* V */ t2 == r4_Add(t1, dedvb), /* V */ t3 == r4_Add(t2, dedrU), /* V */ t4 == r4_Add(t3, dedsV), /* V */ deduf == r4_Scale(dedu, f), /* V */ dedvf == r4_Scale(dedv, f), /* V */ c1 == r4_Add(deduf, dedr), /* V */ d1 == r4_Add(dedvf, deds) /* V */ ){ dedf = t4; deda = c1; dedb = d1; ;} ;} ;} EangleAux; PROCEDURE Eangle( *R,S: r4_t; double *E; VAR EDR,EDS: r4_t; ) == /* Given two vectors "R" and "S" compute the "cos" of the angle between the vectors, the curvature term and the derivatives of energy respect to the two vectors: eeDR and eeDS. */ { ??? m = r4_Norm(R) + Epsilon; ??? n = r4_Norm(S) + Epsilon; ??? o = r4_Dot(R, S); with ( double d = m*n; double q = o/d ){ if ((d!=0.0)){ E = K * (1.0 + q); if ((grad)){ ??? eDq = 1.0 * K; ??? eDo = eDq / d; ??? eDd = - eDq * q / d; ??? eDm = eDd * n; ??? eDn = eDd * m; { EDR = r4_Mix(eDo, S, eDm/m, R); EDS = r4_Mix(eDo, R, eDn/n, S); ;} ;} ;} ;} ;} Eangle; { /* Initialize */ for (i = 0; i < NV; i++){ eDc[i] = zero; } /* Compute energy "e", and the gradient "eDr": */ e = 0.0; for (i = 0; i <= (num.nre-1); i++){ for (j = 0; j <= (num.nce-1); j++){ for (k = j+1; k <= (num.nce-1); k++){ if ((ChilFace[i,j].exists) AND AND (ChilFace[i,k].exists)){ if ((Adjacent(ChilFace[i,j], ChilFace[i,k]))){ ??? u1 = OrgV(ChilFace[i,j].pa).num; with ( double v1 = OrgV(Enext (ChilFace[i,j].pa)).num; double w1 = OrgV(Enext_1(ChilFace[i,j].pa)).num; double u2 = OrgV(ChilFace[i,k].pa).num; double v2 = OrgV(Enext (ChilFace[i,k].pa)).num; double w2 = OrgV(Enext_1(ChilFace[i,k].pa)).num ){ if (((v1 == v2)) AND AND ((u1 == w2))){ AddTerm(v1,u1,w1,u2); }else if (((v1 == w2)) AND AND ((u1 == u2))){ AddTerm(v1,u1,w1,v2); }else if (((v1 == u2)) AND AND ((u1 == v2))){ AddTerm(v1,u1,w1,w2); }else if (((u1 == v2)) AND AND ((w1 == w2))){ AddTerm(u1,w1,v1,u2); }else if (((u1 == w2)) AND AND ((w1 == u2))){ AddTerm(u1,w1,v1,v2); }else if (((u1 == u2)) AND AND ((w1 == v2))){ AddTerm(u1,w1,v1,w2); }else if (((w1 == v2)) AND AND ((v1 == w2))){ AddTerm(w1,v1,u1,u2); }else if (((w1 == w2)) AND AND ((v1 == u2))){ AddTerm(w1,v1,u1,v2); }else if (((w1 == u2)) AND AND ((v1 == v2))){ AddTerm(w1,v1,u1,w2); }else if (((v1 == w2)) AND AND ((u1 == v2))){ AddTerm(v1,u1,w1,u2); }else if (((v1 == v2)) AND AND ((u1 == u2))){ AddTerm(v1,u1,w1,w2); }else if (((v1 == u2)) AND AND ((u1 == w2))){ AddTerm(v1,u1,w1,v2); }else if (((u1 == w2)) AND AND ((w1 == v2))){ AddTerm(u1,w1,v1,u2); }else if (((u1 == v2)) AND AND ((w1 == u2))){ AddTerm(u1,w1,v1,w2); }else if (((u1 == u2)) AND AND ((w1 == w2))){ AddTerm(u1,w1,v1,v2); }else if (((w1 == w2)) AND AND ((v1 == v2))){ AddTerm(w1,v1,u1,u2); }else if (((w1 == v2)) AND AND ((v1 == u2))){ AddTerm(w1,v1,u1,w2); }else if (((w1 == u2)) AND AND ((v1 == w2))){ AddTerm(w1,v1,u1,v2); ;} ;} ;} ;} ;} ;} ;} ;} ;} } /* END Eval */ PROCEDURE Adjacent(f1,f2: Face) : bool_t == VAR count : INTEGER = 0; { ??? u1 = OrgV(f1.pa).num; ??? v1 = OrgV(Enext (f1.pa)).num; ??? w1 = OrgV(Enext_1(f1.pa)).num; ??? u2 = OrgV(f2.pa).num; ??? v2 = OrgV(Enext (f2.pa)).num; ??? w2 = OrgV(Enext_1(f2.pa)).num; { if ((u1 == u2) || (u1 == v2) || (u1 == w2)){ count = count + 1;; } if ((v1 == u2) || (v1 == v2) || (v1 == w2)){ count = count + 1;; } if ((w1 == u2) || (w1 == v2) || (w1 == w2)){ count = count + 1;; } if ((count == 2)){ return TRUE }else{ return FALSE ;} ;} } /* END Adjacent */ PROCEDURE Name(<* UNUSED); erg: T): char *== { return "Curv2D()" } /* END Name */ { ;} OldCurvature2D. // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // #FILE cur/libenergy/JUNK/OldCurvature2D.h #ifndef OldCurvature2D_H #define OldCurvature2D_H [!!! JUNK !!!] // #INCLUDE TYPE T <: Public; double Public = Energy.T OBJECT METHODS init(): T; ; } /* The "Curvature 2D" energy measures the outer dihedral angles between adjacent faces. This energy is computed only for faces that: - exist, - belong to the same original face, and - are adjacent to one common edge. The energy depends on the actual angle "A", as "1 + cos(A); which reaches its minimum when "A==Pi". */ ;} OldCurvature2D. // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // #FILE cur/liboct/JUNK/JSTri.c /* See {JSTri.h} */ // #INCLUDE /* Formerly {JSTriangulation.c}; a junk version of {Triangulation.c}? */ // #TIMESTAMP /* Last edited on 2007-01-25 09:28:10 by stolfi */ #define JSTri_C_author \ "Created 2000 by L. A. P. Lozada, based on 1996 procedures by\n" \ "R. M. Rosi and J. Stolfi.\n" \ "Revisions:\n" \ " 03-08-2000 : Optimized version of the Read and Write procedures by J. Stolfi.\n" \ " 30-08-2000 : Modified the Read and Write procedures for include the case\n" \ " when the cells are octahedra.\n" \ " 07-10-2000 : Added procedure for compute the barycenter of a tetrahedron.\n" \ " 27-01-2001 : Modified for exploding cubic cells." #include #include #include #include // #INCLUDE // #INCLUDE // #INCLUDE // #INCLUDE #include void JSTriWedgeInit(Wedge_t w) { WedgeInit(&(w.w)); w.org[0] = NULL; /* For now, node of primal: C */ w.org[1] = NULL; /* For now, node of dual: C' */ w.org[2] = NULL; /* For now, node of primal: C */ w.org[3] = NULL; /* For now, node of dual: C' */ } /* ===== ELEMENT CREATION ===== */ Place_t JSTriMakeWedge(void) { Wedge_t w = (Wedge_t)malloc(sizeof(WedgeRec_t)), "no mem" wedge_init()) JSTriWedgeInit(w); Place_t a = PlaceInWedge(e, 0); /* [!!! WAS: !!!] w->edge->pa = a; */ /* [!!! WAS: !!!] w->wall->pa = a; */ return a; } /* void SetOrgAll(Place_t @p, Node_t n) VAR Place_t t = a; Place_t tn; { do { SetOrg(t,n); tn = Clock(PrevE(t)); do { SetOrg(tn,n); tn = NextF(tn); } while (!(tn == Clock(PrevE(t)))); t = NextF(t); } while (t != a); } /* END SetOrgAll */ */ /* void SetOrgAll(Place_t @p, Node_t n) Place_t c = a; { ??? nei = Octf.RenumberEdgesOutOfNodes((Place_vec_t){a}); { for(i = 0; i < nei.nel; i++) { ??? b = nei[i]; { c = b; do { SetOrg(c,n); c = NextF(c) } while (c != b); } } } } /* END SetOrgAll */ void SetOrgAll(Place_t @p, Node_t n) { SetOrgAll(a,n); } /* END SetOrgAll */ Node_t Pneg(Place_t @p) { return Org(Srot(a)); } /* END Pneg */ void SetPneg(Place_t @p, Node_t n) { SetOrg(Srot(a), n) } /* END SetPneg */ void SetPnegOfWall(Place_t @p, Node_t n) Place_t t = a; { do { SetPneg(t, n); t = PrevE(t); } while (t != a); } /* END SetPnegOfWall */ void SetPnegOfNearbyWalls(Place_t @p, Node_t n) Place_t t = a; { SetPnegOfWall(t,n); do { SetPnegOfWall(Clock(PrevF(t)),n); t = PrevE(t); } while (t != a); } /* END SetPnegOfNearbyWalls */ Node_t Ppos(Place_t @p) { return Pneg(Clock(a)); } /* END Ppos */ void SetPpos(Place_t @p, Node_t n) { SetOrg(Tors(a), n) } /* END SetPpos */ void SetPposOfWall(Place_t @p, Node_t n) Place_t t = a; { do { SetPpos(t, n); t = PrevE(t); } while (t != a); } /* END SetPposOfWall */ void SetPposQUX(Place_t @p, Node_t n) Place_t t = Clock(a); { SetPnegOfNearbyWalls(t,n); do { SetPnegOfNearbyWalls(PrevF(t),n); t = PrevE(t); } while (!(t == Clock(a))); } /* END SetPposQUX */ Node OrgV(Place_t @p) { return (Node_t)(Org(a)); } /* END OrgV */ Node DesV(Place_t @p) { return OrgV(Clock(a)); } /* END DesV */ Cell PnegP(Place_t @p) { return (Cell_t)(Pneg(a)); } /* END PnegP */ Cell PposP(Place_t @p) { return PnegP(Clock(a)); } /* END PposP */ FourNodes_t TetraNegNodes(Place_t @p) /* Valid for the versions Spin(a), Clock(a) and SpinClock(a). */ { assert(Pneg(a)!=NULL); Node_t p = OrgV(a); Node_t q = OrgV(NextE(a)); Node_t r = OrgV(PrevE(a)); Node_t s = OrgV(PrevE(PrevF(a))); { assert(Pneg(a) == Ppos(PrevE(PrevF(a)))); assert(p != q) && (q != r) && (r != s); return (FourNodes_t){{p,q,r,s}}; } } /* END TetraNegNodes */ FourWalls_t TetraWalls(Place_t @p) { assert(Pneg(a)!=NULL); Wall_t f0 = PWall(a); Wall_t f1 = PWall(PrevF(a)); Wall_t f2 = PWall(PrevF(NextE(a))); Wall_t f3 = PWall(PrevF(PrevE(a))); { assert(f0 != f1) && (f1 != f2) && (f2 != f3); return (FourWalls_t){{f0,f1,f2,f3}}; } } /* END TetraWalls */ SixEdges_t TetraEdges(Place_t @p) { assert(Pneg(a)!=NULL); Edge_t e0 = PEdge(a); Edge_t e1 = PEdge(NextE(a)); Edge_t e2 = PEdge(PrevE(a)); Edge_t e3 = PEdge(NextE(PrevF(a))); Edge_t e4 = PEdge(PrevE(PrevF(a))); Edge_t e5 = PEdge(NextE(PrevF(NextE(a)))); { assert(e0#e1) && (e1#e2) && (e2#e3) && (e3#e4) && (e4#e5); return (SixEdges_t){{e0,e1,e2,e3,e4,e5}}; } } /* END TetraEdges */ ThreeEdges_t TriangEdges(Place_t @p) { Edge_t e0 = PEdge(a); Edge_t e1 = PEdge(NextE(a)); Edge_t e2 = PEdge(PrevE(a)); { assert(e0 != e1) && (e1 != e2); return (ThreeEdges_t){{e0,e1,e2}}; } } /* END TriangEdges */ FourNodes_t TetraPosNodes(Place_t @p) /* Valid for the versions Spin(a), Clock(a) and SpinClock(a). */ { assert(Ppos(a)!=NULL); Node_t p = OrgV(a); Node_t q = OrgV(NextE(a)); Node_t r = OrgV(PrevE(a)); Node_t s = OrgV(PrevE(NextF(a))); { assert(Ppos(a) == Pneg(PrevE(NextF(a)))); return (FourNodes_t){{p,q,r,s}}; } } /* END TetraPosNodes */ FiveNodes_t TetraNegPosNodes(Place_t @p) == /* Valid for the versions Spin(a), Clock(a) and SpinClock(a). We change the otion OrgV by Org such as, this procedure can be used in the dual space but the assertion is valid only in the primal space. */ { ??? p = Org(a); ??? q = Org(NextE(a)); ??? r = Org(PrevE(a)); ??? s = Org(PrevE(PrevF(a))); ??? t = Org(PrevE(NextF(a))); { assert(Ppos(a) == Pneg(PrevE(NextF(a))) AND Pneg(a) == Ppos(PrevE(PrevF(a))) ); return (FiveNodes_t){p,q,r,s,t}; } } /* END TetraNegPosNodes */ /* ========================== CONSTRUCTION TOOLS =============== */ EightPlaces_t MakeTetraTopo(uint nx, uint ny) VAR uint wedgeCount = 0; uint cellCount = 0; Place_t MakeTriangle(void) /* Make one triangular wall and set of the three places with the same wall component. */ { Place_t a = MakeWedge(); Place_t b = MakeWedge(); Place_t c = MakeWedge(); Wall_t f = PWall(a); Node_t u = MakeNode(); Node_t v = MakeNode(); Node_t w = MakeNode(); { PWedge(a)->num = wedgeCount; wedgeCount++; SetOrg(a, u); SetOrg(Clock(a),v); PWedge(b)->num = wedgeCount; wedgeCount++; SetNextE(a,b); SetWallInfo(b,f); SetOrg(b,v); SetOrg(Clock(b),w); PWedge(c)->num = wedgeCount; wedgeCount++; SetNextE(b,c); SetWallInfo(c,f); SetOrg(c, w); SetOrg(Clock(c), Org(a)); return a; } } /* END MakeTriangle */ Place_t AddTwoTriangles(Place_t @p) /* Build a new tetrahedral cell by insertion of two triangular walls and operations SpliceWalls. The argument "a" is the @place return by MakeTriangle(). The @place "f" is return by MakeCell(). */ Place_t *c; { c = NextE(NextF(a)); SetEdgeInfo(a, PEdge(PrevE(c))); ??? f = MakeTriangle(), g == MakeTriangle(); { SetNextF(PrevE(a), PrevE(f)); SetEdgeInfo(PrevE(f), PEdge(PrevE(a))); SetOrgAll(PrevE(a), Org(PrevE(a))); SetOrgAll(Clock(PrevE(a)), Org(Clock(PrevE(a)))); SetNextF(Clock(f), NextE(c)); SetEdgeInfo(Clock(f), PEdge(NextE(c))); SetOrgAll(NextE(c), Org(NextE(c))); SetOrgAll(Clock(NextE(c)), Org(Clock(NextE(c)))); SetNextF(NextE(g), NextE(f)); SetEdgeInfo(NextE(g), PEdge(NextE(f))); SetOrgAll(NextE(f), Org(NextE(f))); SetOrgAll(Clock(NextE(f)), Org(Clock(NextE(f)))); SetNextF(g, c); SetEdgeInfo(g, PEdge(c)); SetOrgAll(c, Org(c)); SetOrgAll(Clock(c), Org(Clock(c))); SetNextF(PrevE(g), Clock(NextE(a))); SetEdgeInfo(PrevE(g), PEdge(Clock(NextE(a)))); SetOrgAll(NextE(a), Org(NextE(a))); SetOrgAll(Clock(NextE(a)), Org(Clock(NextE(a)))); Cell_t p = MakeCell(); { p->num = cellCount; cellCount++; SetPnegOfNearbyWalls(f,p); } return f; } } /* END AddTwoTriangles */ Place_t FinishCellRow(Place_t @p) /* Adds one new row of tetradedral cells and return the @place wedge belong to the cell more rigth. */ { for (col = 0; col < nx; col++) { a = AddTwoTriangles(a); } return a; } FinishCellRow; VAR Place_t t,b; EightPlaces_t ca; { /* ======================= create bottom row of triangles ================= */ for (col = 0; col < nx; col++) { Place_t c = MakeTriangle(); { if (col == 0) { t = c; } else { SpliceWalls(NextE(b), Clock(PrevE(c))); SetEdgeInfo(NextE(b), PEdge(Clock(PrevE(c)))); SetOrgAll(NextE(b), Org(NextE(b))); SetOrgAll(Clock(NextE(b)), Org(Clock(NextE(b)))); } b = c; } } ca.p[2] = t; ca.p[7] = Spin(Clock(b)); for (row = 0; row < ny; row++) { Place_t a = MakeTriangle(); { if (row == 0){ ca.p[1] = Clock(PrevE(a)); } if (row == ny-1){ ca.p[4] = Spin(PrevE(a)); } SpliceWalls(Clock(a), Clock(PrevE(t))); SetOrgAll(PrevE(t), Org(PrevE(t))); SetOrgAll(Clock(PrevE(t)), Org(Clock(PrevE(t)))); Place_t aa = FinishCellRow(a); { t = PrevF(t); SetOrg(PrevE(t), Org(PrevE(a))); b = PrevF(b); SetOrg(Clock(NextE(b)), Org(PrevE(aa))); if (row == 0){ ca.p[0] = Clock(Spin(PrevE(aa))); } if (row == ny-1){ ca.p[5] = PrevE(aa); } } } } ca.p[3] = Spin(t); ca.p[6] = Clock(b); return ca; } /* END MakeTetraTopo */ void EmphasizeTetrahedron(Place_t @p, Place_t b, uint n) void HideNode(Node_t v) /* Hide the node "v". */ { v->exists = FALSE; } HideNode; void HideEdge(Place_t @p) /* Hide the edge "PEdge(a)".*/ { Edge_t e = PEdge(a); { e->exists = FALSE; } } HideEdge; void HideWall(Place_t @p) /* Hide the wall "PWall(a)".*/ { Wall_t f = PWall(a) { f->exists = FALSE; } } HideWall; void HideWallsOfEdge(Place_t @p) /* Hide the ring wall "PWall(a)". */ Place_t an = PrevF(a); VAR { for (j = 1; j < n; j++) { HideWall(an); an = PrevF(an); } } HideWallsOfEdge; VAR ta,Place_t tb[100+1]; { ta[0] = a; tb[0] = b; for (i = 1; i < n; i++) { ta[i] = Clock(PrevE(NextF(NextE(ta[i-1])))); tb[i] = Clock(PrevE(NextF(NextE(tb[i-1])))); } for (i = 0; i < n; i++) { HideWallsOfEdge(ta[i]); HideWallsOfEdge(tb[i]); } for (i = 1; i < n; i++) { HideNode(OrgV(tb[i])); HideNode(OrgV(ta[i])); Place_t dn = tb[i]; VAR { for (j = 0; j <= (n); j++) { HideEdge(PrevE(dn)); dn = PrevF(dn); } } } for (i = 0; i <= (n-2); i++) { HideEdge(NextE(ta[i])); HideEdge(NextE(NextF(ta[i]))); } } /* END EmphasizeTetrahedron */ Place_t Glue( Place_t @p,b; uint n; bool_t setorg = TRUE; ) /* The @place "a" and "b" have the same Orientation and Spin bits, such as, after of the glue procedure performs: Pneg(a) == Pneg(b). */ VAR ta,Place_t tb[100+1]; { /* sanity check */ assert(n >= 1); ta[0] = a; tb[0] = b; if (n > 1) { for (i = 1; i < n; i++) { ta[i] = Clock(PrevE(PrevF(NextE(ta[i-1])))); tb[i] = Clock(PrevE(NextF(NextE(tb[i-1])))); assert(ta[i]!=a); assert(tb[i]!=b); } } Meld(b, a); /* Update the edge slots for i==0 */ SetRingEdgeInfo(a, PEdge(a)); SetRingEdgeInfo(NextE(a), PEdge(NextE(a))); SetRingEdgeInfo(PrevE(a), PEdge(PrevE(a))); if (setorg) { /* Update the node slots for i==0 */ SetOrgAll(a, Org(a)); SetOrgAll(Clock(a), Org(Clock(a))); SetOrgAll(NextE(a), Org(NextE(a))); SetOrgAll(Clock(NextE(a)), Org(Clock(NextE(a)))); SetOrgAll(PrevE(a), Org(PrevE(a))); SetOrgAll(Clock(PrevE(a)), Org(Clock(PrevE(a)))); } /* Update the cell slots for i==0 */ SetPneg(a, Pneg(b)); SetPneg(PrevE(a), Pneg(PrevE(b))); SetPneg(NextE(a), Pneg(NextE(b))); for (i = 1; i < n; i++) { Meld(tb[i],ta[i]); /* Update the edge slots */ SetRingEdgeInfo(ta[i], PEdge(ta[i])); SetRingEdgeInfo(NextE(ta[i]), PEdge(NextE(ta[i]))); SetRingEdgeInfo(PrevE(ta[i]), PEdge(PrevE(ta[i]))); if (setorg) { /* Update the node slots */ SetOrgAll(ta[i], Org(ta[i])); SetOrgAll(Clock(ta[i]), Org(Clock(ta[i]))); SetOrgAll(NextE(ta[i]), Org(NextE(ta[i]))); SetOrgAll(Clock(NextE(ta[i])), Org(Clock(NextE(ta[i])))); SetOrgAll(PrevE(ta[i]), Org(PrevE(ta[i]))); SetOrgAll(Clock(PrevE(ta[i])), Org(Clock(PrevE(ta[i])))); } /* Update the cell slots */ ??? f = Pneg(tb[i]); ??? g = Pneg(PrevE(tb[i])); ??? h = Pneg(NextE(tb[i])); { SetPneg(ta[i],f); SetPneg(PrevE(ta[i]),g); SetPneg(NextE(ta[i]),h); } } if (setorg) { SetOrgAll(Clock(PrevE(PrevF(PrevE(ta[n-1])))),Org(PrevE(a))); } return ta[n-1]; } /* END Glue */ /* ===== GLOBAL PROCEDURES ===== */ VAR Node_vec_t seenNode = Node_vec_new(INIT_STACK_SIZE); uint_vec_t nodeOnum = uint_vec_new(INIT_STACK_SIZE); uint seenNodeCount = 0; void MarkNode(Node_t n) { Node_vec_expand(&seenNode,seenNodeCount); uint_vec_expand(&nodeOnum,seenNodeCount); seenNode[k] = seenNodeCount; nodeOnum[k] = seenNodeCount->num; n->num = seenNodeCount; seenNodeCount++ } void NodeIsMarked(Node_t n) : bool_t { return (n->num < seenNodeCount)) && ((seenNode[n->num] == n); } /* END @{Node->?}IsMarked */ void UnmarkMarkedNodes() { while (seenNodeCount > 0){ ??? k = seenNodeCount-1, n == seenNode[k]; { n->num = nodeOnum[k]; seenNodeCount--; } } } /* END UnmarkMarked@{Node->?}s */ void EnumNodes(Place_t @p, visit: VisitProc) CONST INIT_STACK_SIZE == 1000; VAR Place_vec_t pstack = Place_vec_new(INIT_STACK_SIZE); uint top; /* top of pstack */ void Stack(Place_t c) Place_t cn = c; VAR Place_t dn; { do { Place_vec_expand(&pstack,top); pstack.el[top] = cn; top++; dn = Clock(PrevE(cn)); do { Place_vec_expand(&pstack,top); pstack.el[top] = dn; top++; dn = NextF(dn); } while (!(( dn == Clock(PrevE(cn))))); cn = NextF(cn); } while (cn != c); } Stack; void VisitAndMark(Place_t c) /* If org(c) is diferent that the origins @places stacked, stack the @place "c". */ Place_t *cn; { ??? n = Org(c); { if ((n!=NULL) && (NOT NodeIsMarked(n))) { visit(c); MarkNode(n); Stack(c); cn = c; do { Stack(cn); cn = Onext(cn); } while (!( (cn == c))); } } } VisitAndMark; uint *seen; { top = 0; seen = 0; assert(seenNodeCount == 0); VisitAndMark(a); while (seen < top) { Place_t b = (Place_t){festack.el[seen], bstack.el[seen]}; { VisitAndMark(NextF(b)); VisitAndMark(NextF(NextE(b))); VisitAndMark(NextE(b)); VisitAndMark(PrevF(PrevE(b))); } seen++; } UnmarkMarkedNodes(); } /* END EnumNodes */ uint NumberNodes(Place_t @p) uint n = 0; void Visit(Place_t c) void Visit_(Place_t c) Place_t *cn,dn; { ??? v = Org(c); { cn = c; do { ??? vn = Org(cn); { assert(vn == v); vn->num = n; dn = Clock(PrevE(cn)); do { ??? wn = Org(dn); { assert(wn == vn); wn->num = n; dn = NextF(dn); } } while (!(( dn == Clock(PrevE(cn))))); cn = NextF(cn); } } while (!( cn == c)); } } Visit_; Place_t cn = c; VAR { do { Visit_(cn); cn = Onext(cn); } while (cn != c); n++; } Visit; void VisitDual(Place_t c) void VisitDual_(Place_t c) VAR Place_t cn; p,pn: Node_t; { p = Org(c); if (p != NULL) { cn = c; do { pn = Org(cn); if (pn != NULL) { assert(pn == p); pn->num = n; } cn = NextF(cn); } while (!( cn == c)); } } VisitDual_; Place_t cn = c; VAR { VisitDual_(cn); VisitDual_(Clock(PrevE(cn))); n++; } VisitDual; { if ((Octf.DualBit(a) == 0)) { EnumNodes(a, Visit); } else { EnumNodes(a, VisitDual); } return n; } /* END NumberNodes */ Topology MakeTopology(Place_t @p, uint bdr) Topology top ; INTEGER euler; { top->NV = NumberNodes(a); top->node = Node_vec_new(top->NV); top->wedge = Octf.RenumberWedges((Place_vec_t){a}); top->NFE = ((top->wedge^).nel); top->edge = Octf.RenumberEdges((Place_vec_t){top->wedge[0]}); top->wall = Octf.RenumberWalls((Place_vec_t){top->wedge[0]}); top->out = Place_vec_new(top->NV); top->NE = ((top->edge^).nel); top->NF = ((top->wall^).nel); top->NP = NumberNodes(Srot(a)); top->cell = NEW(REF ARRAY OF Cell, top->NP); top->region = Place_vec_new(top->NP); top->der = Octf.DegreeOfEdge(top->wall[0].pa); top->bdr = bdr; for (i = 0; i < top->NFE; i++) { VAR Place_t c = top->wedge[i]; Node_t v,p; uint vi,qi; { for (k = 0; k < 2; k++) { v = Org(c); vi = v->num; top->node[vi] = v; top->out[vi] = c; p = Pneg(c); if (p != NULL) { qi = p->num; top->cell[qi] = p; top->region[qi] = Srot(c); assert(Pneg(c) == Org(top->region[qi])); } c = Clock(c); } } } ??? m = NumDigits(top->NFE)+1; { fprintf(stderr, "\n"); fprintf(stderr, "nv = " & Fmt.Pad(Fmt.Int(top->NV),m) & "\n"); fprintf(stderr, "ne = " & Fmt.Pad(Fmt.Int(top->NE),m) & "\n"); fprintf(stderr, "nf = " & Fmt.Pad(Fmt.Int(top->NF),m) & "\n"); fprintf(stderr, "np = " & Fmt.Pad(Fmt.Int(top->NP),m) & "\n"); fprintf(stderr, "nfe = " & Fmt.Pad(Fmt.Int(top->NFE),m) & "\n"); } /* Euler's Number */ euler = top->NV-top->NE+top->NF-top->NP; if (euler == 0) { fprintf(stderr,"\nThe Map performs the Euler's Number == 0\n"); } else { fprintf(stderr,"\nThe Euler's Number: " & Fmt.Int(euler) & "\n"); } return top; } /* END MakeTopology */ uint DegreeOfNode(Place_t @p) { ??? @{edge->?} = Octf.RenumberEdgesOutOfNodes((Place_vec_t){a}); ??? degree = ((@{edge->?}^).nel); { return degree; } } /* END DegreeOfNode */ AdjacencyMatrix *MakeAdjacencyMatrix(Topology *top)== m = NEW(REF ARRAY OF ARRAY OF bool_t, top->NV, top->NV); { for (i = 0; i < top->NV; i++){ ??? mi = m[i]; { for (j = 0; j < top->NV; j++) { mi[j] = FALSE } } } for (ei = 0; ei < top->NFE; ei++) { ??? a = top->wedge[ei]; { ??? i = Org(a)->num, j == Org(Clock(a))->num; { m[i,j] = TRUE; m[j,i] = TRUE; } } } return m; } /* END MakeAdjacencyMatrix */ bool_t TriviallyIsomorphic(Topology_t *ta, Topology_t *tb) { if ((ta.NV!=tb.NV ) || (ta.NF!=tb.NF ) || (ta.NE!=tb.NE ) || (ta.NP!=tb.NP ) || (ta.NFE!=tb.NFE )){ return FALSE; } int NFE = ta.NFE; { for (i = 0; i < NFE; i++) { VAR sPlace_t @p = PWedge(ta)->i]; sb: Place_t = PWedge(tb)->i]; { for (r = 0; r < 4; r++) { ??? na = Org(sa)->num; ??? nb = Org(sb)->num; { if (na != nb){ return FALSE;} } ??? za = GetPlaceNum_t(NextF(sa)); ??? zb = GetPlaceNum_t(NextF(sb)); { if (za != zb){ return FALSE;} } sa = Srot(sa); sb = Srot(sb); } } } } return TRUE; } /* END TriviallyIsomorphic */ void CheckOutAndRegion(Topology *top) { for (i = 0; i < top->NV; i++){ ??? v = top->node[i]; ??? e = top->out[i]; { assert(v->num == i); assert(Org(e) == v); } } for (i = 0; i < top->NP; i++) { ??? p = top->cell[i]; ??? r = top->region[i]; { assert(p->num == i); assert(Org(r) == p); } } } /* END CheckOutAndRegion */ void GetVariableNodes(Topology_t *top, VAR vr: ARRAY OF bool_t) { for (i = 0; i < top->NV; i++){ ??? v = top->node[i]; { vr[i] = NOT v.fixed; } } } /* END GetVariableNodes */ /* ===== GEOMETRIC TOOLS ===== */ void InitCoords( rand_t *coins; Coords_t *c; float radius = 1.0) { ??? r = ((double)radius); { for (i = 0 TO (c.nel-1)) { c[i] = (r4_t) { coins.longreal(-r, r), coins.longreal(-r, r), coins.longreal(-r, r), coins.longreal(-r, r) } } } } /* END InitCoords */ Coords_t *GenCoords(Topology_t *t) { ??? coins = NEW(Random.Default).init(TRUE); ??? r = r4_vec_new(t.NV); ??? c = r^; { for (i = 0 TO (c.nel-1)) { c[i] = (r4_t) { coins.longreal(-1.0, +1.0), coins.longreal(-1.0, +1.0), coins.longreal(-1.0, +1.0), coins.longreal(-1.0, +1.0) } } return r; } } /* END GenCoords */ r4_t Barycenter(Topology_t *top, Coords_t *c) B: r4_t = (r4_t){0.0, ..}; uint n = 0; { for (i = 0 TO (c.nel-1)){ ??? v = top->node[i]; { if (v->exists){ B = r4_Add(B, c[i]); N++; } } } return r4_Scale(1.0/((double)N), B); } /* END Barycenter */ double MeanNodeDistance(Topology_t *top, *c: Coords_t) S: double = 0.0; uint n = 0; { for (i = 0 TO (c.nel-1)){ ??? v = top->node[i]; { if (v->exists) { S = S + r4_NormSqr(c[i]); N++ } } } return Math.sqrt(S/FLOAT(N,double)); } /* END MeanNodeDistance */ double MeanEdgeLength(Topology_t *top, *c: Coords_t) S: double = 0.0; uint n = 0; { for (i = 0; i < top->NE; i++){ ??? e = top->edge[i]; { if (e->exists) { ??? o = e.node[0]->num, d == e.node[1]->num; { S = S + r4_DistSqr(c[o], c[d]) } N++; } } } return Math.sqrt(S/FLOAT(N,double)); } /* END MeanEdgeLength */ void Displace(Topology_t *top, r4_t *d, Coords_t *c) { for (i = 0 TO (c.nel-1)){ if (top->node[i]->exists) { ??? vc = c[i]; { vc = r4_Add(vc, d);} } } } /* END Displace */ void Scale(Topology_t *top, double s, Coords_t *c) { for (i = 0 TO (c.nel-1)){ if (top->node[i]->exists) { ??? vc = c[i]; { vc = r4_Scale(s, vc); } } } } /* END Scale */ void NormalizeNodeDistance(Topology_t *top, Coords_t *c) { ??? b = Barycenter(top, c); { Displace(top, r4_Neg(b), c); } ??? s = MeanNodeDistance(top, c); { Scale(top, 1.0/s, c); } } /* END NormalizeNodeDistance */ void NormalizeEdgeLengths(Topology_t *top, Coords_t *c) { ??? b = Barycenter(top, c); { Displace(top, r4_Neg(b), c); } ??? s = MeanEdgeLength(top, c); { Scale(top, 1.0/s, c); } } /* END NormalizeEdgeLengths */ r4_t WallCross(Place_t @p, Coords_t *c) { if ((! PWall(a)->exists)){ return (r4_t){1.0, 0.0, 0.0, 0.0}; } else { Node_t ov = OrgV(a); Node_t dv = OrgV(Clock(a)); Node_t pv = OrgV(PrevE(a)); Node_t qv = OrgV(PrevE(NextF(a))); Node_t rv = OrgV(PrevE(PrevF(a))); { ??? o = c[ov->num]; ??? d = c[dv->num]; ??? p = c[pv->num]; ??? q = c[qv->num]; ??? r = c[rv->num]; ??? n1 = LR4Extras.Cross(r4_Sub(p,o), r4_Sub(r,o), r4_Sub(d,o)); ??? n2 = LR4Extras.Cross(r4_Sub(p,o), r4_Sub(d,o), r4_Sub(q,o)); { return r4_Add(n1,n2); } } } } WallCross; r4_t WallNormal(Place_t @p, Coords_t *c) { ??? n = WallCross(a,c); ??? s = r4_Norm(n); { if (s == 0.0) { return (r4_t){1.0, 0.0, 0.0, 0.0}; } else { return r4_Scale(1.0/s, n); } } } /* END WallNormal */ r4_t PolyCross(Place_t @p, Coords_t *c) { if ((! PnegP(a)->exists)){ return (r4_t){1.0, 0.0, 0.0, 0.0}; } else { Node_t ov = OrgV(a); Node_t dv = OrgV(Clock(a)); Node_t pv = OrgV(PrevE(a)); Node_t rv = OrgV(PrevE(PrevF(a))); { ??? o = c[ov->num]; ??? d = c[dv->num]; ??? p = c[pv->num]; ??? r = c[rv->num]; { return LR4Extras.Cross(r4_Sub(p,o), r4_Sub(r,o), r4_Sub(d,o)); } } } } PolyCross; r4_t PolyNormal(Place_t @p, Coords_t *c) { ??? n = PolyCross(a, c); ??? s = r4_Norm(n); { if (s == 0.0) { return (r4_t){1.0, 0.0, 0.0, 0.0}; } else { return r4_Scale(1.0/s, n); } } } /* END PolyNormal */ r4_t WallBarycenter(Place_t @p, Coords_t *c) VAR Place_t @po; n : uint = 0; sum = (r4_t){0.0, ..}; { ao = a; do { ??? aoc = c[OrgV(ao)->num]; { sum = r4_Add(sum, aoc); n++; ao = NextE(ao); } } while (ao != a); if (n == 0) { return sum; } else { return r4_Scale(1.0/FLOAT(n,double), sum); } } /* END WallBarycenter */ r4_t TetraBarycenter(Place_t @p, Coords_t *c) uint n = 0; sum = (r4_t){0.0, ..}; { ??? tetra = TetraNegNodes(a); { for (i = 0; i < 4; i++) { ??? aoc = c[tetra[i]->num]; { sum = r4_Add(sum, aoc); n++; } } } return r4_Scale(1.0/FLOAT(n,double), sum); } /* END TetraBarycenter */ r4_t EdgeCross(Place_t @p, Coords_t *c) sum: r4_t = (r4_t){0.0,..}; Place_t @po; { Node_t uv = OrgV(a); ??? u = c[uv->num]; { if ((! PEdge(a)->exists)) { return (r4_t){1.0, 0.0, 0.0, 0.0}; } else { ao = a; do { Place_t an = PrevF(ao); Node_t dv = OrgV(Clock(ao)); Node_t pv = OrgV(PrevE(ao)); Node_t rv = OrgV(PrevE(an)); { ??? d = c[dv->num]; ??? p = c[pv->num]; ??? r = c[rv->num]; ??? n = LR4Extras.Cross(r4_Sub(p,u), r4_Sub(d,u), r4_Sub(r,u)); { if (ao == a){ sum = n; }else{ sum = r4_Add(sum,n); } } ao = an; } } while (ao != a); return sum; } } } /* END EdgeCross */ r4_t EdgeNormal(Place_t @p, Coords_t *c) { ??? n = EdgeCross(a, c); ??? s = r4_Norm(n); { if (s == 0.0) { return (r4_t){1.0, 0.0, 0.0, 0.0}; } else { return r4_Scale(1.0/s, n); } } } /* END EdgeNormal */ r4_t NodeCross(Place_t @p, Coords_t *c, *top : Topology) sum,n: r4_t = (r4_t){0.0, ..}; { Node_t uv = OrgV(a); ??? poly = StarOfNode(a,top); { if (! uv->exists) { return (r4_t){1.0, 0.0, 0.0, 0.0}; } else { for(i = 0; i < poly.nel; i++) { if ((poly[i][0]->exists) && (poly[i][1]->exists AND poly[i][2]->exists) && (poly[i][3]->exists )){ ??? u = c[poly[i][0]->num]; ??? d = c[poly[i][1]->num]; ??? p = c[poly[i][2]->num]; ??? r = c[poly[i][3]->num]; { n= LR4Extras.Cross(r4_Sub(p,u), r4_Sub(d,u), r4_Sub(r,u)); if (i == 0){ sum = n; }else{ sum = r4_Add(sum,n); } } } } } return sum; } } /* END NodeCross */ r4_t NodeNormal(Place_t @p, Coords_t *c, Topology_t *top) { ??? n = NodeCross(a, c, top); ??? s = r4_Norm(n); { if (s == 0.0) { return (r4_t){1.0, 0.0, 0.0, 0.0}; } else { return r4_Scale(1.0/s, n); } } } /* END NodeNormal */ r4_t NeighborBarycenter(Node_vec_t *neigh; *c: Coords_t) uint n = 0; b = (r4_t){0.0, ..}; { ??? r = ((neigh^).nel); { for (i = 0; i < r; i++) { ??? v = neigh^[i]; { if (v->exists){ b = r4_Add(b, c[v->num]); n++;} } } } if (n == 0) { return b; } else { return r4_Scale(1.0/((double)n), b) } } /* END NeighborBarycenter */ Node_vec_t *NeighborNode(Place_t @p, Topology_t *top) == REF ARRAY OF Triv tri ; Node_vec_t vstack = Node_vec_new(INIT_STACK_SIZE); uint nstack = 0; uint c; bool_t Present(Node_t v) /* Return TRUE if "v" is on the stack, FALSE c.c */ uint nstack1 = nstack; VAR { while (nstack1 > 0) { nstack1 = nstack1 - 1; if (vstack.el[nstack1] == v){ return TRUE; } } return FALSE; } Present; { ??? poly = StarOfNode(a,top); ??? np = NumberPolyOfStar(poly); { tri = NEW(REF ARRAY OF Triv, np); for(i = 0; i < poly.nel; i++) { c = 0; for (j = 0; j < 4; j++) { if (( poly[i][j]!=OrgV(a) )){ tri[i,c] = poly[i][j]; c++ } } assert(c == 3); } for(j = 0; j < poly.nel; j++) { for (l = 0; l <= (2); l++) { if ((! Present(tri[j,l]) )){ vstack.el[nstack] = tri[j,l]; nstack++; } } } Node_vec_trim(&vstack,nstack); } } /* END NeighborNode */ Node_vec_t *Neighbors(Place_t @p, Topology_t *top) == VAR Node_vec_t vstack = Node_vec_new(INIT_STACK_SIZE); uint nstack = 0; uint v; { v = OrgV(a)->num; for (i = 0; i < top->NE; i++) { ??? ei = top->edge[i]; ??? ei0 = ei.node[0]->num; ??? ei1 = ei.node[1]->num; { if ((v == ei0) || ((v == ei1) )){ if (v == ei0 ){ vstack.el[nstack] = top->node[ei1]; nstack++; } else { vstack.el[nstack] = top->node[ei0]; nstack++; } } } } Node_vec_trim(&vstack,nstack); } /* END Neighbors */ uint NumberNeighborNode(Node_vec_t *neigh) { ??? n = ((neigh^).nel); { return n; }; } /* END NumberNeighborNode */ Quadv_vec_t *StarOfNode(Place_t @p, Topology_t *top)== PROCEDURE NodePoly(*top : Topology): REF ARRAY OF Quadp == /* Return four places such as the origins corresponding to nodes extremes of each tetrahedron of Triangulation. */ REF ARRAY OF Quadp poly = NEW(REF ARRAY OF Quadp, top->NP); { for (i = 0; i < top->NP; i++) { ??? da = top->region[i]; Place_t a0 = Tors(da); Place_t db = Clock(PrevE(da)); Place_t b0 = Tors(db); { /*assert(DegreeOfNode(a0) == 4); */ /* consuming time */ Place_t a1 = NextE(a0); Place_t a2 = NextE(a1); Place_t a3 = PrevE(b0); { if ((NextE(a2)!=a0)) { fprintf(stderr, "\nTriangulation: This topology isn't" \ " a triangulation\n"); assert(NextE(a2) == a0); } assert(Pneg(a0)->num == i); if ((Pneg(b0)!=NULL)) { assert(Pneg(a0) == Pneg(b0)); } assert(a0!=b0); poly[i] = Quadp{a0,a1,a2,a3}; } } } return poly; } NodePoly; poly1 : REF Quadv_vec_t = Quadv_vec_new(top->NP); n : uint = 0; { ??? poly = NodePoly(top), v == OrgV(a); { for(i = 0; i < poly.nel; i++) { if (((OrgV(poly[i][0]) == v)) || ((OrgV(poly[i][1]) == v) OR (OrgV(poly[i][2]) == v)) || ((OrgV(poly[i][3]) == v) )){ for (j = 0; j < 4; j++) { if ((OrgV(poly[i][j]) == v )){ Node_t uv = OrgV(poly[i][j]); Node_t dv = OrgV(NextE(poly[i][j])); Node_t pv = OrgV(PrevE(poly[i][j])); Node_t rv = OrgV(PrevE(PrevF(poly[i][j]))); { poly1[n] = Quadv{uv,dv,pv,rv}; n++; } } } } } Quadv_vec_trim(&poly1,n); } } /* END StarOfNode */ uint NumberPolyOfStar(Quadv_vec_t *qv) { ??? n = (qv->nel); { return n; }; } /* END NumberPolyOfStar */ r4_vec_t *ComputeAllNodeNormals( *Topology_t *top; Coords_t *c; )== { ??? rvn = r4_vec_new(top->NV); ??? vn = rvn^; { for (i = 0; i < top->NV; i++) { vn[i] = NodeNormal(top->out[i], c, top) } return rvn } } /* END ComputeAllNodeNormals */ r4_vec_t *ComputeAllEdgeNormals(Topology_t *top, Coords_t *c; )== { ??? rvn = r4_vec_new(top->NE); ??? vn = rvn^; { for (i = 0; i < top->NE; i++) { vn[i] = EdgeNormal(top->edge[i].pa, c); } return rvn } } /* END ComputeAllEdgeNormals */ r4_vec_t *ComputeAllWallNormals(Topology_t *top, Coords_t *c; )== { ??? rvn = r4_vec_new(top->NF); ??? vn = rvn^; { for (i = 0; i < top->NF; i++) { vn[i] = WallNormal(top->wall[i].pa, c); } return rvn } } /* END ComputeAllWallNormals */ r4_vec_t *ComputeAllCellNormals(Topology_t *top, Coords_t *c; )== { ??? rvn = r4_vec_new(top->NP); ??? vn = rvn^; { for (i = 0; i < top->NP; i++) { vn[i] = PolyNormal(Tors(top->region[i]), c); } return rvn } } /* END ComputeAllCellNormals */ CellTopology MakeCellTopology(Place_t @p) uint *ne, nv; b: Place_t; ptop: CellTopology; { ??? star = Octf.RenumberEdgesOutOfNodes((Place_vec_t){a})^; { /* Gather walls: */ ptop.NF = (star.nel); ptop.fRef = Place_vec_new(ptop.NF); ne = 0; for (i = 0; i < ptop.NF; i++) { ??? side = Octf.Dual(star[i]); { ptop.fRef[i] = side; b = side; do { ne++; Edge_t e = PEdge(b); { e.xmark = FALSE; OrgV(b).xmark = FALSE; } b = NextE(b) } while (b != side); } } /* Gather @{edge->?}s: */ assert(ne % 2 == 0); ptop.NE = ne DIV 2; ptop.eRef = Place_vec_new(ptop.NE); ne = 0; for (i = 0; i < ptop.NF; i++) { ??? side = Octf.Dual(star[i]); { b = side; do { Edge_t e = PEdge(b); { if (! e.xmark) { ptop.eRef[ne] = b; ne++; e.xmark = TRUE; } } b = NextE(b) } while (!( b == side)); } } /* Gather nodes: */ ptop.NV = 2 + ptop.NE - ptop.NF; ptop.vRef = Place_vec_new(ptop.NV); nv = 0; for (i = 0; i < ptop.NE; i++) { ??? eu = ptop.eRef[i]; Node_t u = OrgV(eu); Place_t ev = Clock(PrevF(eu)); Node_t v = OrgV(ev); { if (! u.xmark){ ptop.vRef[nv] = eu; nv++; u.xmark = TRUE; } if (! v.xmark){ ptop.vRef[nv] = ev; nv++; v.xmark = TRUE; } PEdge(eu).xmark = FALSE; }; } for (i = 0; i < ptop.NV; i++){ Node_t u = OrgV(ptop.vRef[i]); { u.xmark = FALSE; } } } return ptop; } /* END MakeCellTopology */ /* =============== INPUT/OUTPUT =========================== */ CONST BoolChars == Mis.BoolChars; AlphaChars == Mis.AlphaChars; TopCom_t ReadTopology(FILE *rd) VAR Topology_t *top; char *cmt; CHAR n; { /* Topology */ ReadHeader(rd,"topology","99-08-25"); cmt = Mis.ReadCommentsJS(rd, '|'); /* Element counts: */ Lex.Skip(rd, cs = AlphaChars); top->NV = fget_int(rd); Lex.Skip(rd, cs = AlphaChars); top->NE = fget_int(rd); Lex.Skip(rd, cs = AlphaChars); top->NF = fget_int(rd); Lex.Skip(rd, cs = AlphaChars); top->NP = fget_int(rd); Lex.Skip(rd, cs = AlphaChars); top->NFE = fget_int(rd); Lex.Skip(rd, cs = AlphaChars); top->der = fget_int(rd); Lex.Skip(rd, cs = AlphaChars); top->bdr = fget_int(rd); Lex.Skip(rd); ??? map = NEW(REF ARRAY OF Octf.wedge, top->NFE)^; { /* Create node records: */ top->node = Node_vec_new(top->NV); top->out = Place_vec_new(top->NV); for (i = 0; i < top->NV; i++) { top->node[i] = MakeNode(); top->node[i]->num = i; } /* Create @{edge->?}s records */ top->edge = NEW(REF ARRAY OF Octf.edge, top->NE); for (i = 0; i < top->NE; i++) { top->edge[i] = Octf.Make@{Edge->?}(); top->edge[i]->num = i; } /* Create wall records: */ top->wall = Wall_vec_new(top->NF); for (i = 0; i < top->NF; i++) { top->wall[i] = Octf.MakeWall(); top->wall[i]->num = i; } /* Create cells records: */ top->cell = NEW(REF ARRAY OF Cell, top->NP); top->region = Place_vec_new(top->NP); for (i = 0; i < top->NP; i++) { top->cell[i] = MakeCell(); top->cell[i]->num = i; } /* Create wedge records: */ top->wedge = NEW(REF ARRAY OF Octf_Place_t, top->NFE); for (i = 0; i < top->NFE; i++) { top->wedge[i] = MakeWedge(); PWedge(top->wedge[i])->num = i; map[i] = PWedge(top->wedge[i]); } /* Read @{edge->?} data: */ EVAL Mis.ReadCommentsJS(rd, '|'); for (j = 0; j < top->NE; j++) { Lex.Skip(rd); ??? ne = fget_int(rd); /* index to @{edge->?} */ ??? e == top->edge[ne]; ??? pa == Octf.ReadPlace(rd,map); { assert(ne == j); e->pa = pa } } Lex.Skip(rd); /* Read wall data: */ EVAL Mis.ReadCommentsJS(rd, '|'); for (j = 0; j < top->NF; j++) { Lex.Skip(rd); ??? nf = fget_int(rd); /* index to wall */ ??? f = top->wall[nf]; ??? pa = Octf.ReadPlace(rd,map); { assert(nf == j); f->pa = pa; } } Lex.Skip(rd); /* Read wedge data: */ EVAL Mis.ReadCommentsJS(rd, '|'); for (j = 0; j < top->NFE; j++) { Lex.Skip(rd); ??? nfe = fget_int(rd); /* index to wedge */ ??? fe = NARROW(PWedge(top->wedge[nfe]), wedge); { assert(nfe == j); assert(top->wedge[nfe].bits == 0); fe->num = nfe; Octf.ReadWedge(rd, fe, map); for (k = 0; k < 4; k++) { Lex.Skip(rd); n = Rd.GetChar(rd); if (n!='-') { Rd.UnGetChar(rd); ??? m = fget_int(rd); ??? pa = PlaceInWedge(fe, 2*k); ??? vf = fe.org[k]; { if ((Rd.GetChar(rd) == 'v')) { vf = top->node[m]; top->out[m] = pa; } else { vf = top->cell[m]; top->region[m] = pa; } } } } Lex.Skip(rd); ??? nf = fget_int(rd); { fe.wall = top->wall[nf]; } assert(Rd.GetChar(rd) == 'f'); Lex.Skip(rd); ??? ne = fget_int(rd); { fe.edge = top->edge[ne] } assert(Rd.GetChar(rd) == 'e'); } } ReadFooter(rd,"topology"); fprintf(stderr, "rebuilding node tables...\n"); RebuildNodeTables(top); return (TopCom_t){top, cmt}; } } /* END ReadTopology */ void RebuildNodeTables(Topology *top) Place_t *b; { for (i = 0; i < top->NE; i++){ ??? pa = top->edge[i].pa; { assert(PEdge(pa) == top->edge[i]); assert(top->edge[i]->num == i); top->edge[i].node[0] = OrgV(pa); top->edge[i].node[1] = OrgV(Clock(pa)) } } for (i = 0; i < top->NF; i++) { ??? a = top->wall[i].pa; ??? deg = Octf.DegreeOfWall(a); ??? vr = NEW(REF ARRAY OF Node_t, deg); { assert(PWall(a) == top->wall[i]); assert(top->wall[i]->num == i); top->wall[i].node = vr; b = a; for (j = 0; j < deg; j++) { vr[j] = Org(b); b = NextE(b); } } } for (i = 0; i < top->NP; i++) { ??? a = top->region[i]; ??? ptop = MakeCellTopology(a); ??? vr = NEW(REF ARRAY OF Node_t, ptop.NV); { assert(Org(a) == top->cell[i]); assert(top->cell[i]->num == i); top->cell[i].node = vr; for (j = 0; j < ptop.NV; j++) { vr[j] = Org(ptop.vRef[j]) } } } } /* END RebuildNodeTables */ PROCEDURE ReadMaterials( FILE *rd; *Topology_t *top; bool_t ro_te = FALSE; ) { /* Materials */ ReadHeader(rd,"materials","99-08-25"); /* Read node data Materials: */ EVAL Mis.ReadCommentsJS(rd, '|'); for (j = 0; j < top->NV; j++) { Lex.Skip(rd); ??? nv = fget_int(rd); /* index to node */ ??? v = top->node[nv]; { assert(nv == j); v->num = nv; Lex.Skip(rd); v->exists = Mis.ReadBool(rd); Lex.Skip(rd); v.fixed = Mis.ReadBool(rd); Lex.Skip(rd); ??? cc = v.color; { cc[0] = Lex.Real(rd); Lex.Skip(rd); cc[1] = Lex.Real(rd); Lex.Skip(rd); cc[2] = Lex.Real(rd); } Lex.Skip(rd); ??? tt = v.transp; { tt[0] = Lex.Real(rd); Lex.Skip(rd); tt[1] = Lex.Real(rd); Lex.Skip(rd); tt[2] = Lex.Real(rd); } Lex.Skip(rd); v.radius = Lex.Real(rd); Lex.Skip(rd); v.label = Rd.GetText(rd,2); } } Lex.Skip(rd); /* Read @{edge->?} data materials: */ EVAL Mis.ReadCommentsJS(rd, '|'); for (j = 0; j < top->NE; j++) { Lex.Skip(rd); ??? ne = fget_int(rd); /* index to @{edge->?} */ ??? e = top->edge[ne]; { assert(ne == j); e->num = ne; Lex.Skip(rd); e->exists = Mis.ReadBool(rd); Lex.Skip(rd); ??? cc = e.color; { cc[0] = Lex.Real(rd); Lex.Skip(rd); cc[1] = Lex.Real(rd); Lex.Skip(rd); cc[2] = Lex.Real(rd); } Lex.Skip(rd); ??? tt = e.transp; { tt[0] = Lex.Real(rd); Lex.Skip(rd); tt[1] = Lex.Real(rd); Lex.Skip(rd); tt[2] = Lex.Real(rd); } Lex.Skip(rd); e.radius = Lex.Real(rd); Lex.Skip(rd); e.degenerate = Mis.ReadBool(rd); Lex.Skip(rd); ??? n = Rd.GetChar(rd); { if (n!='-') { Rd.UnGetChar(rd); e->root = fget_int(rd); } else { e->root = -1; } } } } Lex.Skip(rd); /* Read wall data materials: */ EVAL Mis.ReadCommentsJS(rd, '|'); for (j = 0; j < top->NF; j++) { Lex.Skip(rd); ??? nf = fget_int(rd); /* index to wall */ ??? f = top->wall[nf]; { assert(nf == j); f->num = nf; Lex.Skip(rd); f->exists = Mis.ReadBool(rd); Lex.Skip(rd); ??? cc = f.color; { cc[0] = Lex.Real(rd); Lex.Skip(rd); cc[1] = Lex.Real(rd); Lex.Skip(rd); cc[2] = Lex.Real(rd); } Lex.Skip(rd); ??? tt = f.transp; { tt[0] = Lex.Real(rd); Lex.Skip(rd); tt[1] = Lex.Real(rd); Lex.Skip(rd); tt[2] = Lex.Real(rd); } Lex.Skip(rd); f.degenerate = Mis.ReadBool(rd); Lex.Skip(rd); ??? n = Rd.GetChar(rd); { if (n!='-') { Rd.UnGetChar(rd); f->root = fget_int(rd); } else { f->root = -1; } } } } Lex.Skip(rd); /* Read cell data materials: */ if (top->NP!=0) { EVAL Mis.ReadCommentsJS(rd, '|'); } for (j = 0; j < top->NP; j++) { Lex.Skip(rd); ??? np = fget_int(rd); /* index to cell */ ??? p = top->cell[np]; { assert(np == j); p->num = np; Lex.Skip(rd); p->exists = Mis.ReadBool(rd); Lex.Skip(rd); ??? cc = p.color; { cc[0] = Lex.Real(rd); Lex.Skip(rd); cc[1] = Lex.Real(rd); Lex.Skip(rd); cc[2] = Lex.Real(rd); } Lex.Skip(rd); ??? tt = p.transp; { tt[0] = Lex.Real(rd); Lex.Skip(rd); tt[1] = Lex.Real(rd); Lex.Skip(rd); tt[2] = Lex.Real(rd); } Lex.Skip(rd); p.degenerate = Mis.ReadBool(rd); if (ro_te) { Lex.Skip(rd); ??? n = Rd.GetChar(rd); { if (n!='-') { Rd.UnGetChar(rd); p->root = fget_int(rd); } else { p->root = -1; } } } } } ReadFooter(rd,"materials"); Rd.Close(rd); CheckOutAndRegion(top); } /* END ReadMaterials */ TopCom_t ReadToTaMa(char *name, bool_t ro_te /* DF FALSE */) /* Where ro_te meaning "root tetrahedron". */ VAR TopCom_t tc; { ??? ntp = name & ".tp"; ??? rtp = FileRd.Open(ntp); { fprintf(stderr, "reading " & ntp & "\n"); tc = ReadTopology(rtp); fprintf(stderr, "OK leio tp\n"); Rd.Close(rtp); } ??? nma = name & ".ma"; ??? rma = FileRd.Open(nma); { fprintf(stderr, "reading " & nma & "\n"); ReadMaterials(rma, tc.top, ro_te); Rd.Close(rma); } fprintf(stderr, "OK leio ma\n"); return tc; } /* END ReadToTaMa */ Coords_t *ReadState(char *name)== <* FATAL Rd.Failure, Thread.Alerted,FloatMode.Trap, Lex.Error, OSError.E); Coords_t c; char *cmt; uint nv; { ??? rs = FileRd.Open(name & ".st"); { /* Read Headers File Formats*/ ReadHeader(rs,"state","99-08-25"); /* Element counts: */ Lex.Skip(rs, cs = AlphaChars); nv = fget_int(rs); Lex.Skip(rs); cmt = Mis.ReadCommentsJS(rs, '|'); c = NEW(REF Coords_t, nv); /* Read node data state: */ for (j = 0; j < nv; j++) { Lex.Skip(rs); ??? nv = fget_int(rs); { ??? cv = c[nv]; { cv[0] = Lex.LongReal(rs); Lex.Skip(rs); cv[1] = Lex.LongReal(rs); Lex.Skip(rs); cv[2] = Lex.LongReal(rs); Lex.Skip(rs); cv[3] = Lex.LongReal(rs); } } } ReadFooter(rs,"state"); Rd.Close(rs); return c; } } /* END ReadState */ void WriteState(char *name, Topology_t *top, Coords_t *c, char *cmt /* DF "" */) <* FATAL OSError.E); { ??? st = FileWr.Open(name & ".st"); ??? vWidth = NumDigits(top->NV - 1); { void WriteCoord(x: double) { Wr.PutText(st, Fmt.Pad(Fmt.LongReal(x, Fmt.Style.Sci, prec = 3), 7)); } WriteCoord; void WritePoint(*c: r4_t) { WriteCoord(c[0]); Wr.PutText(st, " "); WriteCoord(c[1]); Wr.PutText(st, " "); WriteCoord(c[2]); Wr.PutText(st, " "); WriteCoord(c[3]); } WritePoint; { WriteHeader(st,"state","99-08-25"); Wr.PutText(st, "nodes "); Wr.PutText(st, Fmt.Int(top->NV) & "\n"); if ((! Text.Empty(cmt))) { Mis.WriteCommentsJS(st, cmt & "\n", '|') } Mis.WriteCommentsJS(st, "\nNode data:\n", '|'); for (i = 0; i < top->NV; i++) { ??? v = top->node[i]; { /* state */ Wr.PutText(st, Fmt.Pad(Fmt.Int(v->num), vWidth)); Wr.PutText(st, " "); WritePoint(c[v->num]); Wr.PutText(st, "\n"); } } } WriteFooter(st, "state"); fclose(st); } } /* END WriteState */ void WriteMaterials(char *name, Topology_t *top, char *cmt /* DF "" */, bool_t ro_te = FALSE; /* root tetrahedron */ ) <* FATAL Wr.Failure, Thread.Alerted, OSError.E); INTEGER *pWidth; { if (top->NP == 0){ pWidth = 2 }else{ pWidth = NumDigits(top->NP-1); } ??? ma = FileWr.Open(name & ".ma"); ??? vWidth = NumDigits(top->NV - 1); ??? fWidth = NumDigits(top->NF - 1); ??? eWidth = NumDigits(top->NE - 1); { void WriteIntensity(r: REAL) { Wr.PutText(ma, Fmt.Real(r, Fmt.Style.Fix, prec = 2)); } WriteIntensity; void WriteColor(*c: frgb_t) { WriteIntensity(c[0]); Wr.PutText(ma, " "); WriteIntensity(c[1]); Wr.PutText(ma, " "); WriteIntensity(c[2]); } WriteColor; void WriteRadius(r: REAL) { if (r == 0.02 ){ Wr.PutText(ma, "0.020"); } else { Wr.PutText(ma,Fmt.Real(r, prec = 4)); } } WriteRadius; void WriteLabel(char *label) { Wr.PutText(ma, label); } WriteLabel; { WriteHeader(ma,"materials","99-08-25"); if ((! Text.Empty(cmt))) { WriteCommentsJS(ma, cmt & "\n", '|') } ??? m = NumDigits(top->NP); { WriteCommentsJS(ma,"nodes " & Fmt.Pad(Fmt.Int(top->NV),m), '|'); WriteCommentsJS(ma,"@{edge->?}s " & Fmt.Pad(Fmt.Int(top->NE),m), '|'); WriteCommentsJS(ma,"walls " & Fmt.Pad(Fmt.Int(top->NF),m), '|'); WriteCommentsJS(ma,"cells " & Fmt.Int(top->NP), '|'); } WriteCommentsJS(ma, "\nNode data:\n", '|'); for (i = 0; i < top->NV; i++) { ??? v = top->node[i]; { /* materials */ Wr.PutText(ma, Fmt.Pad(Fmt.Int(v->num), vWidth)); Wr.PutText(ma, " "); Wr.PutText(ma, Fmt.Char(BoolChars[v->exists])); Wr.PutText(ma, " "); Wr.PutText(ma, Fmt.Char(BoolChars[v.fixed])); Wr.PutText(ma, " "); WriteColor(v.color); Wr.PutText(ma, " "); WriteColor(v.transp); Wr.PutText(ma, " "); WriteRadius(v.radius); Wr.PutText(ma, " "); WriteLabel(v.label); Wr.PutText(ma, "\n"); } } WriteCommentsJS(ma, "\n@{Edge->?} data:\n", '|'); for (i = 0; i < top->NE; i++) { ??? e = top->edge[i]; { /* materials */ Wr.PutText(ma, Fmt.Pad(Fmt.Int(e->num), eWidth)); Wr.PutText(ma, " "); Wr.PutText(ma, Fmt.Char(BoolChars[e->exists])); Wr.PutText(ma, " "); WriteColor(e.color); Wr.PutText(ma, " "); WriteColor(e.transp); Wr.PutText(ma, " "); WriteRadius(e.radius); Wr.PutText(ma, " "); Wr.PutText(ma, Fmt.Char(BoolChars[e.degenerate])); Wr.PutText(ma, " "); if (e->root == -1) { Wr.PutText(ma, " - "); } else { Wr.PutText(ma, Fmt.Pad(Fmt.Int(e->root), eWidth)); } Wr.PutText(ma, "\n"); } } WriteCommentsJS(ma, "\nWall data:\n", '|'); for (i = 0; i < top->NF; i++) { ??? f = top->wall[i]; { Wr.PutText(ma, Fmt.Pad(Fmt.Int(f->num), fWidth)); Wr.PutText(ma, " "); Wr.PutText(ma, Fmt.Char(BoolChars[f->exists])); Wr.PutText(ma, " "); WriteColor(f.color); Wr.PutText(ma, " "); WriteColor(f.transp); Wr.PutText(ma, " "); Wr.PutText(ma, Fmt.Char(BoolChars[f.degenerate])); Wr.PutText(ma, " "); if (f->root == -1) { Wr.PutText(ma, " - "); } else { Wr.PutText(ma, Fmt.Pad(Fmt.Int(f->root), fWidth)); } Wr.PutText(ma, "\n"); } } if (top->NP!=0) { WriteCommentsJS(ma, "\nCell data:\n", '|'); } for (i = 0; i < top->NP; i++) { ??? p = top->cell[i]; { Wr.PutText(ma, Fmt.Pad(Fmt.Int(p->num), pWidth)); Wr.PutText(ma, " "); Wr.PutText(ma, Fmt.Char(BoolChars[p->exists])); Wr.PutText(ma, " "); WriteColor(p.color); Wr.PutText(ma, " "); WriteColor(p.transp); Wr.PutText(ma, " "); Wr.PutText(ma, Fmt.Char(BoolChars[p.degenerate])); if (ro_te) { Wr.PutText(ma, " "); if (p->root == -1) { Wr.PutText(ma, " - "); } else { Wr.PutText(ma, Fmt.Pad(Fmt.Int(p->root), pWidth)); } } Wr.PutText(ma, "\n"); } } } WriteFooter(ma, "materials"); fclose(ma); } } /* END WriteMaterials */ void FindDegeneracies(Topology *top) VAR bool_t allTriangles, allTetrahedra; vei, vej: REF ARRAY OF INTEGER; vfi, vfj: REF ARRAY OF INTEGER; vpi, vpj: REF ARRAY OF INTEGER; { /* Clean the marks for the attribute degenerate */ for (i = 0; i < top->NE; i++) { ??? ei = top->edge[i]; { ei.degenerate = FALSE } } /* Now, find @{edge->?} degeneracies */ vei = NEW(REF ARRAY OF INTEGER, 2); vej = NEW(REF ARRAY OF INTEGER, 2); for (i = 0; i < top->NE; i++) { for (j = i+1; j < top->NE; j++) { ??? ei = NARROW(top->edge[i], @{Edge->?}); ??? ej = NARROW(top->edge[j], @{Edge->?}); { for (k = 0; k < 2; k++) { vei[k] = ei.node[k]->num; vej[k] = ej.node[k]->num } Mis.InsertionSort(1,vei); Mis.InsertionSort(1,vej); if ((vei[0] == vej[0]) && (vei[1] == vej[1])) { ei.degenerate = TRUE; ej.degenerate = TRUE; } } } } /* Clean the marks for the attribute degenerate */ allTriangles = TRUE; for (i = 0; i < top->NF; i++) { ??? fi = top->wall[i]; ??? di = ((fi.node^).nel); { fi.degenerate = (di!=3); allTriangles = allTriangles) && ((di == 3) } } if (! allTriangles){ return; } /* Now, find degenerate @places of triangles */ vfi = NEW(REF ARRAY OF INTEGER, 3); vfj = NEW(REF ARRAY OF INTEGER, 3); for (i = 0; i < top->NF; i++) { for (j = i+1; j < top->NF; j++) { ??? fi = top->wall[i]; ??? fj = top->wall[j]; { for (k = 0; k <= (2); k++) { vfi[k] = fi.node^[k]->num; vfj[k] = fj.node^[k]->num } Mis.InsertionSort(2,vfi); Mis.InsertionSort(2,vfj); if ((vfi[0] == vfj[0]) && (vfi[1] == vfj[1]) && (vfi[2] == vfj[2])) { fi.degenerate = TRUE; fj.degenerate = TRUE } } } } /* Clean the marks for the attribute degenerate */ allTetrahedra = TRUE; for (i = 0; i < top->NP; i++) { ??? pi = top->cell[i]; ??? di = ((pi.node^).nel); { pi.degenerate = (di!=4); allTetrahedra = allTetrahedra) && ((di == 4) } } if (! allTetrahedra){ return; } /* Now, find cell degeneracies */ vpi = NEW(REF ARRAY OF INTEGER, 4); vpj = NEW(REF ARRAY OF INTEGER, 4); for (i = 0; i < top->NP; i++) { for (j = i+1; j < top->NP; j++) { ??? pi = top->cell[i]; ??? pj = top->cell[j]; { for (k = 0; k < 4; k++) { vpi[k] = pi.node^[k]->num; vpj[k] = pj.node^[k]->num } Mis.InsertionSort(3,vpi); Mis.InsertionSort(3,vpj); if ((vpi[0]==vpj[0]) && (vpi[1]==vpj[1]) && (vpi[2]==vpj[2]) && (vpi[3]==vpj[3])) { pi.degenerate = TRUE; pj.degenerate = TRUE; } } } } } /* END FindDegeneracies */ void WriteTopology(char *name, Topology_t *top, char *cmt /* DF " " */) <* FATAL Wr.Failure, Thread.Alerted, OSError.E); INTEGER *pWidth; { if (top->NP == 0){ pWidth= 2 }else{ pWidth = NumDigits(top->NP-1); } ??? tp = FileWr.Open(name & ".tp"); ??? vWidth = NumDigits(top->NV - 1); ??? eWidth = NumDigits(top->NE - 1); ??? fWidth = NumDigits(top->NF - 1); ??? feWidth = NumDigits(top->NFE -1); { WriteHeader(tp,"topology","99-08-25"); if ((! Text.Empty(cmt))) { WriteCommentsJS(tp, cmt & "\n", '|') } ??? m = NumDigits(top->NFE); { Wr.PutText(tp, "nodes "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NV),m) & "\n"); Wr.PutText(tp, "@{edge->?}s "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NE),m) & "\n"); Wr.PutText(tp, "walls "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NF),m) & "\n"); Wr.PutText(tp, "cells "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NP),m) & "\n"); Wr.PutText(tp, "wedges "); Wr.PutText(tp, Fmt.Int(top->NFE) & "\n"); Wr.PutText(tp, "der "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->der),m) & "\n"); Wr.PutText(tp, "bdr "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->bdr),m) & "\n"); } WriteCommentsJS(tp, "\n@{Edge->?} data:\n", '|'); for (i = 0; i < top->NE; i++) { ??? e = top->edge[i]; { Wr.PutText(tp, Fmt.Pad(Fmt.Int(e->num), eWidth)); Wr.PutText(tp, " "); Octf.PrintPlace(tp, e.pa, eWidth+1, TRUE); } } WriteCommentsJS(tp, "\nWall data:\n", '|'); for (i = 0; i < top->NF; i++) { ??? f = top->wall[i]; { Wr.PutText(tp, Fmt.Pad(Fmt.Int(f->num), fWidth)); Wr.PutText(tp, " "); Octf.PrintPlace(tp, f.pa, fWidth+1, TRUE); } } WriteCommentsJS(tp, "\nwedge data:\n", '|'); for (i = 0; i < top->NFE; i++) { ??? fe = NARROW(PWedge(top->wedge[i]), wedge); { Wr.PutText(tp, Fmt.Pad(Fmt.Int(fe->num), feWidth)); Wr.PutText(tp, " "); Octf.PrintWedge(tp, fe, feWidth); Wr.PutText(tp, " "); for (j = 0; j < 4; j++) { ??? n = fe.org[j]; { TYPECASE n OF | NULL ==> for (i = 0; i <= (pWidth-2); i++){ Wr.PutText(tp," "); } Wr.PutText(tp, " - "); | Node(v) ==> Wr.PutText(tp, Fmt.Pad(Fmt.Int(v->num), vWidth) & "v "); | Cell(p) ==> Wr.PutText(tp, Fmt.Pad(Fmt.Int(p->num), pWidth) & "p "); } else { /* nothing */ } } } Wr.PutText(tp, " "); ??? f = fe.wall; { Wr.PutText(tp, Fmt.Pad(Fmt.Int(f->num), fWidth)); } Wr.PutText(tp, "f "); ??? e = fe.edge; { Wr.PutText(tp, Fmt.Pad(Fmt.Int(e->num), eWidth)); } Wr.PutText(tp, "e\n"); } } WriteFooter(tp, "topology"); fclose(tp); } } /* END WriteTopology */ void WriteDualTopology(char *name, Topology_t *top, char *cmt /* DF " " */) <* FATAL Wr.Failure, Thread.Alerted, OSError.E); void PrintDualPlace_t(FILE *wr, Place_t @p; feWidth: uint) { Wr.PutText(wr,Fmt.Pad(Fmt.Int(PWedge(a)->num), feWidth) & ":" \ Fmt.Int((Octf.SrotBits(a)+ 3) % 4) & ":" & Fmt.Int(Octf.SpinBit(a))); } PrintDualPlace_t; void PrintDualwedge(FILE *wr, Wedge_t w, feWidth: uint) Place_t *b; { b = Srot(PlaceInWedge(w, 0)); for (i = 0; i < 4; i++) { PrintDualPlace_t(wr, NextF(b), feWidth); Wr.PutText(wr, " "); b = Srot(b) } } PrintDualwedge; { ??? tp = FileWr.Open(name & ".tp"); ??? vWidth = NumDigits(MAX(1, top->NV - 1)); ??? eWidth = NumDigits(MAX(1,top->NE - 1)); ??? fWidth = NumDigits(MAX(1,top->NF - 1)); ??? pWidth = NumDigits(MAX(1,top->NP - 1)); ??? feWidth = NumDigits(MAX(1,top->NFE -1)); { WriteHeader(tp,"topology","99-08-25"); if ((! Text.Empty(cmt))) { WriteCommentsJS(tp, cmt & "\n", '|') } ??? m = NumDigits(top->NFE); { Wr.PutText(tp, "nodes "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NP),m) & "\n"); Wr.PutText(tp, "@{edge->?}s "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NF),m) & "\n"); Wr.PutText(tp, "walls "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NE),m) & "\n"); Wr.PutText(tp, "cells "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NV),m) & "\n"); Wr.PutText(tp, "wedges "); Wr.PutText(tp, Fmt.Int(top->NFE) & "\n"); Wr.PutText(tp, "der "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(1),m) & "\n"); Wr.PutText(tp, "bdr "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(1),m) & "\n"); } WriteCommentsJS(tp, "\n@{Edge->?} data:\n", '|'); for (i = 0; i < top->NF; i++) { ??? f = top->wall[i]; { Wr.PutText(tp, Fmt.Pad(Fmt.Int(f->num), fWidth)); Wr.PutText(tp, " "); Octf.PrintPlace(tp, f.pa, feWidth); /* [sic] */ Wr.PutText(tp, "\n") } } WriteCommentsJS(tp, "\nWall data:\n", '|'); for (i = 0; i < top->NE; i++) { ??? e = top->edge[i]; { Wr.PutText(tp, Fmt.Pad(Fmt.Int(e->num), eWidth)); Wr.PutText(tp, " "); Octf.PrintPlace(tp, e.pa, feWidth+1); /* [sic] */ Wr.PutText(tp, "\n") } } WriteCommentsJS(tp, "\nwedge data:\n", '|'); for (i = 0; i < top->NFE; i++) { ??? fe = PWedge(top->wedge[i]); { Wr.PutText(tp, Fmt.Pad(Fmt.Int(fe->num), feWidth)); Wr.PutText(tp, " "); PrintDualwedge(tp, fe, feWidth); Wr.PutText(tp, " "); for (j = 0; j < 4; j++) { ??? n = fe.org[(j+1) % 4]; { TYPECASE n OF | NULL ==> for (i = 0; i <= (pWidth-2); i++){ Wr.PutText(tp," "); } Wr.PutText(tp, " - "); | Node(v) ==> Wr.PutText(tp, Fmt.Pad(Fmt.Int(v->num), vWidth) & "p "); | Cell(p) ==> Wr.PutText(tp, Fmt.Pad(Fmt.Int(p->num), pWidth) & "v "); } else { /* nothing */ } } } Wr.PutText(tp, " "); ??? e = fe.edge; { Wr.PutText(tp, Fmt.Pad(Fmt.Int(e->num), eWidth)); } Wr.PutText(tp, "f "); ??? f = fe.wall; { Wr.PutText(tp, Fmt.Pad(Fmt.Int(f->num), fWidth)); } Wr.PutText(tp, "e"); Wr.PutText(tp, "\n"); } } WriteFooter(tp, "topology"); fclose(tp); } } /* END WriteDualTopology */ void WriteTable(char *name, Topology_t *top, char *cmt /* DF "" */, <*UNUSED); bool_t debug = FALSE; ) == <* FATAL Wr.Failure, Thread.Alerted, OSError.E); { RebuildNodeTables(top); ??? tab = FileWr.Open(name & ".tb"); ??? vWidth = NumDigits(MAX(1, top->NV - 1)); ??? eWidth = NumDigits(MAX(1, top->NE - 1)); ??? fWidth = NumDigits(MAX(1, top->NF - 1)); ??? pWidth = NumDigits(MAX(1, top->NP - 1)); { WriteHeader(tab,"table","99-08-25"); if ((! Text.Empty(cmt))) { Mis.WriteCommentsJS(tab, cmt & "\n", '|') } /* Now, writes @{edge->?} nodes */ Mis.WriteCommentsJS(tab,"\n@{Edge->?} data:\n", '|'); for (i = 0; i < top->NE; i++) { ??? e = top->edge[i]; { Wr.PutText(tab, Fmt.Pad(Fmt.Int(e->num), eWidth)); Wr.PutText(tab, " "); Wr.PutText(tab, Fmt.Pad(Fmt.Int(e.node[0]->num),eWidth) & "v "); Wr.PutText(tab, Fmt.Pad(Fmt.Int(e.node[1]->num),eWidth) & "v\n"); } } /* Now, writes wall nodes */ Mis.WriteCommentsJS(tab, "\nWall data:\n", '|'); for (i = 0; i < top->NF; i++) { ??? f = top->wall[i]; ??? fv = f.node^; { Wr.PutText(tab, Fmt.Pad(Fmt.Int(f->num), fWidth)); Wr.PutText(tab, " "); for (j = 0 TO (fv.nel-1)) { Wr.PutText(tab,Fmt.Pad(Fmt.Int(fv[j]->num), vWidth) & "v "); } Wr.PutText(tab, "\n"); } } /* Now, writes cell nodes */ if (top->NP!=0) { Mis.WriteCommentsJS(tab, "\nCell data:\n", '|'); } for (i = 0; i < top->NP; i++) { ??? p = top->cell[i]; ??? pv = p.node^; { Wr.PutText(tab, Fmt.Pad(Fmt.Int(p->num), pWidth)); Wr.PutText(tab, " "); for (j = 0 TO (pv.nel-1)) { Wr.PutText(tab,Fmt.Pad(Fmt.Int(pv[j]->num), vWidth) & "v "); } Wr.PutText(tab, "\n"); } } WriteFooter(tab, "table"); fclose(tab); } } /* END WriteTable */ void WriteStDe( FILE *wr, Coords_t *c, Coords_t *Dc; uint prec /* DF 4 */, char *cmt /* DF "" */) void WriteCoord(x: double) { Wr.PutText(wr, Fmt.LongReal(x, Fmt.Style.Sci, prec = prec)) } WriteCoord; void WritePoint(*p: r4_t) { WriteCoord(p[1]); Wr.PutText(wr, " "); WriteCoord(p[0]); Wr.PutText(wr, " "); WriteCoord(p[2]); Wr.PutText(wr, " "); WriteCoord(p[3]); Wr.PutText(wr, " "); } WritePoint; { int NV = (c.nel); ??? d = NumDigits(NV-1); { WriteHeader(wr,"state-derivatives","99-08-25"); Mis.WriteCommentsJS(wr, "\n" & cmt & "\n",'|'); Wr.PutText(wr, "nodes == " & Fmt.Int(NV) & "\n"); for (i = 0; i < NV; i++) { Wr.PutText(wr, Fmt.Pad(Fmt.Int(i), d) & ": "); WritePoint(c[i]); Wr.PutText(wr, " "); /* if the derivatives of nodes are zero them writes "0 0 0 0" else writes velocites with the format for write points defined here. */ if ((Dc[i][0] == 0.0) && (Dc[i][1] == 0.0) && (Dc[i][2] == 0.0 AND Dc[i][3] == 0.0 )){ Wr.PutText(wr, "0 0 0 0") } else { WritePoint(Dc[i]) } Wr.PutText(wr, "\n"); } WriteFooter(wr, "state-derivatives"); Wr.PutText(wr, "\n"); fflush(wr); } } /* END WriteStDe */ { } Triangulation. /* Copyright © 2000 Universidade Estadual de Campinas (UNICAMP) */ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // #FILE cur/liboct/JUNK/JSTri.h #ifndef JSTri_H #define JSTri_H // #TIMESTAMP /* Last edited on 2007-01-25 03:14:38 by stolfi */ /* Formerly {JSTriangulation.h}; a junk version of {Triangulation.h}? */ // #INCLUDE // #INCLUDE // #INCLUDE // #INCLUDE #include #include Place_t JSTriMakeWedge(void); void JSTriWedgeInit(Wedge_t w); /* Creates a new wedge in a triangulation, intializes it with {JSTriWedgeInit}, returns it base place. */ Place_t JSTriMakeCell(void); void JSTriCellInit(Wedge_t w); /* Creates a new cell record, uniattached. */ Topology_t *JSTriMakeTopology(Place_t a, uint bdr); /* This procedure compute the structure topological of tridimensional complexes without/with border. It assumes that the procedures "NumberVertices", "Octf->numberEdges","Octf->numberFacets","Octf->numberfacetedges" of library "libm3triang" was called. */ void JSTriWriteTable ( char *name, Topology_t *top, char *cmt /* DF " " */, bool_t debug /* DF FALSE */ ); /* This procedure create tables of vertex/edges/faces/polyhedra for one triangulation. */ Node_vec_t *JSTriNeighborNode(Place_t a, Topology_t *top); /* [!!! Marked OBSOLETE in M3 Sources!!!] Neighbor vertices of "OrgV(a)" in a triangulation. */ uint JSTriNumberNeighborNode(Node_vec_t *neigh); /* Return the number neighbor vertex of OrgV(a). */ void JSTriFindDegeneracies(Topology_t *top); /* Finds geometric degeneracies. Update the attribute "degenerate" of the elements: edge, face and polyhedron. */ TopCom_t JSTriReadToTaMa(char *name, bool_t ro_te /* DF FALSE */); /* Reads three disk files created by "WriteTopology", "MakeTopologyTable" and "WriteMaterials". The files must have the given "name" with ".tp", ".tb" and ".ma". */ #endif // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // #FILE cur/liboct/JUNK/JSTriRest.c /* See {JStri.h}. */ // #INCLUDE // #TIMESTAMP /* Last edited on 2007-01-25 03:22:19 by stolfi */ /* Rest of implementation of {Triangulation.c} [JS 2007-01-24-222100 ] */ /* GLOBAL PROCEDURES */ static Node_vec_t seenNode = Node_vec_new(INIT_STACK_SIZE); /* A stack of {Node_t}s used by {NumberNodes} below. */ static uint seenNodeCount = 0; /* During {NumberNodes}, {seenNode[0..seenNodeCount-1]} are the node records renumbered so far. */ static uint_vec_t nodeOnum = uint_vec_new(INIT_STACK_SIZE); /* During {NumberNodes}, {nodeOnum[i]} is the saved node number of {seenNode[i]}, for {i} in {0..seenNodeCount-1}. */ void MarkNode(Node_t n) { /* Stack the node {n}: */ Node_vec_expand(&seenNode,seenNodeCount); seenNode[k] = seenNodeCount; /* Save its old number, and renumber it with position in {seenNode}: */ uint_vec_expand(&nodeOnum,seenNodeCount); nodeOnum[k] = seenNodeCount->num; n->num = seenNodeCount; /* Now we have one more node in {seenNode}: */ seenNodeCount++ } void NodeIsMarked(Node_t n) : bool_t { return ((n->num < seenNodeCount) && (seenNode[n->num] == n)); } void UnmarkMarkedNodes() { while (seenNodeCount > 0){ { int k = seenNodeCount-1; Node_t n = seenNode.el[k]; n->num = nodeOnum.el[k]; seenNodeCount--; } } void EnumNodes(Place_t @p, visit: VisitProc) CONST INIT_STACK_SIZE == 1000; VAR Place_vec_t pstack = Place_vec_new(INIT_STACK_SIZE); uint top; /* top of pstack */ void Stack(Place_t c) Place_t cn = c; VAR Place_t dn; { do { Place_vec_expand(&pstack,top); pstack.el[top] = cn; top++; dn = Clock(PrevE(cn)); do { Place_vec_expand(&pstack,top); pstack.el[top] = dn; top++; dn = NextF(dn); } while (!(( dn == Clock(PrevE(cn))))); cn = NextF(cn); } while (cn != c); } Stack; void VisitAndMark(Place_t c) /* If org(c) is diferent that the origins @places stacked, stack the @place "c". */ Place_t *cn; { ??? n = Org(c); { if ((n!=NULL) && (NOT NodeIsMarked(n))) { visit(c); MarkNode(n); Stack(c); cn = c; do { Stack(cn); cn = Onext(cn); } while (!( (cn == c))); } } } VisitAndMark; uint *seen; { top = 0; seen = 0; assert(seenNodeCount == 0); VisitAndMark(p); while (seen < top) { Place_t b = (Place_t){festack.el[seen], bstack.el[seen]}; { VisitAndMark(NextF(b)); VisitAndMark(NextF(NextE(b))); VisitAndMark(NextE(b)); VisitAndMark(PrevF(PrevE(b))); } seen++; } UnmarkMarkedNodes(); } /* END EnumNodes */ uint NumberNodes(Place_t @p) uint n = 0; void Visit(Place_t c) void Visit_(Place_t c) Place_t *cn,dn; { ??? v = Org(c); { cn = c; do { ??? vn = Org(cn); { assert(vn == v); vn->num = n; dn = Clock(PrevE(cn)); do { ??? wn = Org(dn); { assert(wn == vn); wn->num = n; dn = NextF(dn); } } while (!(( dn == Clock(PrevE(cn))))); cn = NextF(cn); } } while (!( cn == c)); } } Visit_; Place_t cn = c; VAR { do { Visit_(cn); cn = Onext(cn); } while (cn != c); n++; } Visit; void VisitDual(Place_t c) void VisitDual_(Place_t c) VAR Place_t cn; p,pn: Node_t; { p = Org(c); if (p != NULL) { cn = c; do { pn = Org(cn); if (pn != NULL) { assert(pn == p); pn->num = n; } cn = NextF(cn); } while (!( cn == c)); } } VisitDual_; Place_t cn = c; VAR { VisitDual_(cn); VisitDual_(Clock(PrevE(cn))); n++; } VisitDual; { if ((Octf.DualBit(p) == 0)) { EnumNodes(p, Visit); } else { EnumNodes(p, VisitDual); } return n; } /* END NumberNodes */ Topology MakeTopology(Place_t @p, uint bdr) Topology top ; INTEGER euler; { top->NV = NumberNodes(p); top->node = Node_vec_new(top->NV); top->wedge = Octf.RenumberWedges((Place_vec_t){p}); top->NFE = ((top->wedge^).nel); top->edge = Octf.RenumberEdges((Place_vec_t){top->wedge[0]}); top->wall = Octf.RenumberWalls((Place_vec_t){top->wedge[0]}); top->out = Place_vec_new(top->NV); top->NE = ((top->edge^).nel); top->NF = ((top->wall^).nel); top->NP = NumberNodes(Srot(p)); top->cell = NEW(REF ARRAY OF Cell, top->NP); top->region = Place_vec_new(top->NP); top->der = Octf.DegreeOfEdge(top->wall[0].pa); top->bdr = bdr; for (i = 0; i < top->NFE; i++) { VAR Place_t c = top->wedge[i]; Node_t v,p; uint vi,qi; { for (k = 0; k < 2; k++) { v = Org(c); vi = v->num; top->node[vi] = v; top->out[vi] = c; p = Pneg(c); if (p != NULL) { qi = p->num; top->cell[qi] = p; top->region[qi] = Srot(c); assert(Pneg(c) == Org(top->region[qi])); } c = Clock(c); } } } ??? m = NumDigits(top->NFE)+1; { fprintf(stderr, "\n"); fprintf(stderr, "nv = " & Fmt.Pad(Fmt.Int(top->NV),m) & "\n"); fprintf(stderr, "ne = " & Fmt.Pad(Fmt.Int(top->NE),m) & "\n"); fprintf(stderr, "nf = " & Fmt.Pad(Fmt.Int(top->NF),m) & "\n"); fprintf(stderr, "np = " & Fmt.Pad(Fmt.Int(top->NP),m) & "\n"); fprintf(stderr, "nfe = " & Fmt.Pad(Fmt.Int(top->NFE),m) & "\n"); } /* Euler's Number */ euler = top->NV-top->NE+top->NF-top->NP; if (euler == 0) { fprintf(stderr,"\nThe Map performs the Euler's Number == 0\n"); } else { fprintf(stderr,"\nThe Euler's Number: " & Fmt.Int(euler) & "\n"); } return top; } /* END MakeTopology */ uint DegreeOfNode(Place_t @p) { ??? @{edge->?} = Octf.RenumberEdgesOutOfNodes((Place_vec_t){p}); ??? degree = ((@{edge->?}^).nel); { return degree; } } /* END DegreeOfNode */ AdjacencyMatrix *MakeAdjacencyMatrix(Topology *top)== m = NEW(REF ARRAY OF ARRAY OF bool_t, top->NV, top->NV); { for (i = 0; i < top->NV; i++){ ??? mi = m[i]; { for (j = 0; j < top->NV; j++) { mi[j] = FALSE } } } for (ei = 0; ei < top->NFE; ei++) { ??? p = top->wedge[ei]; { ??? i = Org(p)->num, j == Org(Clock(p))->num; { m[i,j] = TRUE; m[j,i] = TRUE; } } } return m; } /* END MakeAdjacencyMatrix */ bool_t TriviallyIsomorphic(Topology_t *ta, Topology_t *tb) { if ((ta.NV!=tb.NV ) || (ta.NF!=tb.NF ) || (ta.NE!=tb.NE ) || (ta.NP!=tb.NP ) || (ta.NFE!=tb.NFE )){ return FALSE; } int NFE = ta.NFE; { for (i = 0; i < NFE; i++) { VAR sPlace_t @p = PWedge(ta)->i]; sb: Place_t = PWedge(tb)->i]; { for (r = 0; r < 4; r++) { ??? na = Org(sa)->num; ??? nb = Org(sb)->num; { if (na != nb){ return FALSE;} } ??? za = GetPlaceNum_t(NextF(sa)); ??? zb = GetPlaceNum_t(NextF(sb)); { if (za != zb){ return FALSE;} } sa = Srot(sa); sb = Srot(sb); } } } } return TRUE; } /* END TriviallyIsomorphic */ void CheckOutAndRegion(Topology *top) { for (i = 0; i < top->NV; i++){ ??? v = top->node[i]; ??? e = top->out[i]; { assert(v->num == i); assert(Org(e) == v); } } for (i = 0; i < top->NP; i++) { ??? p = top->cell[i]; ??? r = top->region[i]; { assert(p->num == i); assert(Org(r) == p); } } } /* END CheckOutAndRegion */ void GetVariableNodes(Topology_t *top, VAR vr: ARRAY OF bool_t) { for (i = 0; i < top->NV; i++){ ??? v = top->node[i]; { vr[i] = NOT v.fixed; } } } /* END GetVariableNodes */ /* ===== GEOMETRIC TOOLS ===== */ void InitCoords( rand_t *coins; Coords_t *c; float radius = 1.0) { ??? r = ((double)radius); { for (i = 0 TO (c.nel-1)) { c[i] = (r4_t) { coins.longreal(-r, r), coins.longreal(-r, r), coins.longreal(-r, r), coins.longreal(-r, r) } } } } /* END InitCoords */ Coords_t *GenCoords(Topology_t *t) { ??? coins = NEW(Random.Default).init(TRUE); ??? r = r4_vec_new(t.NV); ??? c = r^; { for (i = 0 TO (c.nel-1)) { c[i] = (r4_t) { coins.longreal(-1.0, +1.0), coins.longreal(-1.0, +1.0), coins.longreal(-1.0, +1.0), coins.longreal(-1.0, +1.0) } } return r; } } /* END GenCoords */ r4_t Barycenter(Topology_t *top, Coords_t *c) B: r4_t = (r4_t){0.0, ..}; uint n = 0; { for (i = 0 TO (c.nel-1)){ ??? v = top->node[i]; { if (v->exists){ B = r4_Add(B, c[i]); N++; } } } return r4_Scale(1.0/((double)N), B); } /* END Barycenter */ double MeanNodeDistance(Topology_t *top, *c: Coords_t) S: double = 0.0; uint n = 0; { for (i = 0 TO (c.nel-1)){ ??? v = top->node[i]; { if (v->exists) { S = S + r4_NormSqr(c[i]); N++ } } } return Math.sqrt(S/FLOAT(N,double)); } /* END MeanNodeDistance */ double MeanEdgeLength(Topology_t *top, *c: Coords_t) S: double = 0.0; uint n = 0; { for (i = 0; i < top->NE; i++){ ??? e = top->edge[i]; { if (e->exists) { ??? o = e.node[0]->num, d == e.node[1]->num; { S = S + r4_DistSqr(c[o], c[d]) } N++; } } } return Math.sqrt(S/FLOAT(N,double)); } /* END MeanEdgeLength */ void Displace(Topology_t *top, r4_t *d, Coords_t *c) { for (i = 0 TO (c.nel-1)){ if (top->node[i]->exists) { ??? vc = c[i]; { vc = r4_Add(vc, d);} } } } /* END Displace */ void Scale(Topology_t *top, double s, Coords_t *c) { for (i = 0 TO (c.nel-1)){ if (top->node[i]->exists) { ??? vc = c[i]; { vc = r4_Scale(s, vc); } } } } /* END Scale */ void NormalizeNodeDistance(Topology_t *top, Coords_t *c) { ??? b = Barycenter(top, c); { Displace(top, r4_Neg(b), c); } ??? s = MeanNodeDistance(top, c); { Scale(top, 1.0/s, c); } } /* END NormalizeNodeDistance */ void NormalizeEdgeLengths(Topology_t *top, Coords_t *c) { ??? b = Barycenter(top, c); { Displace(top, r4_Neg(b), c); } ??? s = MeanEdgeLength(top, c); { Scale(top, 1.0/s, c); } } /* END NormalizeEdgeLengths */ r4_t WallCross(Place_t @p, Coords_t *c) { if ((! PWall(p)->exists)){ return (r4_t){1.0, 0.0, 0.0, 0.0}; } else { Node_t ov = OrgV(p); Node_t dv = OrgV(Clock(p)); Node_t pv = OrgV(PrevE(p)); Node_t qv = OrgV(PrevE(NextF(p))); Node_t rv = OrgV(PrevE(PrevF(p))); { ??? o = c[ov->num]; ??? d = c[dv->num]; ??? p = c[pv->num]; ??? q = c[qv->num]; ??? r = c[rv->num]; ??? n1 = LR4Extras.Cross(r4_Sub(p,o), r4_Sub(r,o), r4_Sub(d,o)); ??? n2 = LR4Extras.Cross(r4_Sub(p,o), r4_Sub(d,o), r4_Sub(q,o)); { return r4_Add(n1,n2); } } } } WallCross; r4_t WallNormal(Place_t @p, Coords_t *c) { ??? n = WallCross(p,c); ??? s = r4_Norm(n); { if (s == 0.0) { return (r4_t){1.0, 0.0, 0.0, 0.0}; } else { return r4_Scale(1.0/s, n); } } } /* END WallNormal */ r4_t PolyCross(Place_t @p, Coords_t *c) { if ((! PnegP(p)->exists)){ return (r4_t){1.0, 0.0, 0.0, 0.0}; } else { Node_t ov = OrgV(p); Node_t dv = OrgV(Clock(p)); Node_t pv = OrgV(PrevE(p)); Node_t rv = OrgV(PrevE(PrevF(p))); { ??? o = c[ov->num]; ??? d = c[dv->num]; ??? p = c[pv->num]; ??? r = c[rv->num]; { return LR4Extras.Cross(r4_Sub(p,o), r4_Sub(r,o), r4_Sub(d,o)); } } } } PolyCross; r4_t PolyNormal(Place_t @p, Coords_t *c) { ??? n = PolyCross(p, c); ??? s = r4_Norm(n); { if (s == 0.0) { return (r4_t){1.0, 0.0, 0.0, 0.0}; } else { return r4_Scale(1.0/s, n); } } } /* END PolyNormal */ r4_t WallBarycenter(Place_t @p, Coords_t *c) VAR Place_t @po; n : uint = 0; sum = (r4_t){0.0, ..}; { ao = p; do { ??? aoc = c[OrgV(ao)->num]; { sum = r4_Add(sum, aoc); n++; ao = NextE(ao); } } while (ao != p); if (n == 0) { return sum; } else { return r4_Scale(1.0/FLOAT(n,double), sum); } } /* END WallBarycenter */ r4_t TetraBarycenter(Place_t @p, Coords_t *c) uint n = 0; sum = (r4_t){0.0, ..}; { ??? tetra = TetraNegNodes(p); { for (i = 0; i < 4; i++) { ??? aoc = c[tetra[i]->num]; { sum = r4_Add(sum, aoc); n++; } } } return r4_Scale(1.0/FLOAT(n,double), sum); } /* END TetraBarycenter */ r4_t EdgeCross(Place_t @p, Coords_t *c) sum: r4_t = (r4_t){0.0,..}; Place_t @po; { Node_t uv = OrgV(p); ??? u = c[uv->num]; { if ((! PEdge(p)->exists)) { return (r4_t){1.0, 0.0, 0.0, 0.0}; } else { ao = p; do { Place_t an = PrevF(ao); Node_t dv = OrgV(Clock(ao)); Node_t pv = OrgV(PrevE(ao)); Node_t rv = OrgV(PrevE(an)); { ??? d = c[dv->num]; ??? p = c[pv->num]; ??? r = c[rv->num]; ??? n = LR4Extras.Cross(r4_Sub(p,u), r4_Sub(d,u), r4_Sub(r,u)); { if (ao == p){ sum = n; }else{ sum = r4_Add(sum,n); } } ao = an; } } while (ao != p); return sum; } } } /* END EdgeCross */ r4_t EdgeNormal(Place_t @p, Coords_t *c) { ??? n = EdgeCross(p, c); ??? s = r4_Norm(n); { if (s == 0.0) { return (r4_t){1.0, 0.0, 0.0, 0.0}; } else { return r4_Scale(1.0/s, n); } } } /* END EdgeNormal */ r4_t NodeCross(Place_t @p, Coords_t *c, *top : Topology) sum,n: r4_t = (r4_t){0.0, ..}; { Node_t uv = OrgV(p); ??? poly = StarOfNode(p,top); { if (! uv->exists) { return (r4_t){1.0, 0.0, 0.0, 0.0}; } else { for(i = 0; i < poly.nel; i++) { if ((poly[i][0]->exists) && (poly[i][1]->exists AND poly[i][2]->exists) && (poly[i][3]->exists )){ ??? u = c[poly[i][0]->num]; ??? d = c[poly[i][1]->num]; ??? p = c[poly[i][2]->num]; ??? r = c[poly[i][3]->num]; { n= LR4Extras.Cross(r4_Sub(p,u), r4_Sub(d,u), r4_Sub(r,u)); if (i == 0){ sum = n; }else{ sum = r4_Add(sum,n); } } } } } return sum; } } /* END NodeCross */ r4_t NodeNormal(Place_t @p, Coords_t *c, Topology_t *top) { ??? n = NodeCross(p, c, top); ??? s = r4_Norm(n); { if (s == 0.0) { return (r4_t){1.0, 0.0, 0.0, 0.0}; } else { return r4_Scale(1.0/s, n); } } } /* END NodeNormal */ r4_t NeighborBarycenter(Node_vec_t *neigh; *c: Coords_t) uint n = 0; b = (r4_t){0.0, ..}; { ??? r = ((neigh^).nel); { for (i = 0; i < r; i++) { ??? v = neigh^[i]; { if (v->exists){ b = r4_Add(b, c[v->num]); n++;} } } } if (n == 0) { return b; } else { return r4_Scale(1.0/((double)n), b) } } /* END NeighborBarycenter */ Node_vec_t *NeighborNode(Place_t @p, Topology_t *top) == REF ARRAY OF Triv tri ; Node_vec_t vstack = Node_vec_new(INIT_STACK_SIZE); uint nstack = 0; uint c; bool_t Present(Node_t v) /* Return TRUE if "v" is on the stack, FALSE c.c */ uint nstack1 = nstack; VAR { while (nstack1 > 0) { nstack1 = nstack1 - 1; if (vstack.el[nstack1] == v){ return TRUE; } } return FALSE; } Present; { ??? poly = StarOfNode(p,top); ??? np = NumberPolyOfStar(poly); { tri = NEW(REF ARRAY OF Triv, np); for(i = 0; i < poly.nel; i++) { c = 0; for (j = 0; j < 4; j++) { if (( poly[i][j]!=OrgV(p) )){ tri[i,c] = poly[i][j]; c++ } } assert(c == 3); } for(j = 0; j < poly.nel; j++) { for (l = 0; l <= (2); l++) { if ((! Present(tri[j,l]) )){ vstack.el[nstack] = tri[j,l]; nstack++; } } } Node_vec_trim(&vstack,nstack); } } /* END NeighborNode */ Node_vec_t *Neighbors(Place_t @p, Topology_t *top) == VAR Node_vec_t vstack = Node_vec_new(INIT_STACK_SIZE); uint nstack = 0; uint v; { v = OrgV(a)->num; for (i = 0; i < top->NE; i++) { ??? ei = top->edge[i]; ??? ei0 = ei.node[0]->num; ??? ei1 = ei.node[1]->num; { if ((v == ei0) || ((v == ei1) )){ if (v == ei0 ){ vstack.el[nstack] = top->node[ei1]; nstack++; } else { vstack.el[nstack] = top->node[ei0]; nstack++; } } } } Node_vec_trim(&vstack,nstack); } /* END Neighbors */ uint NumberNeighborNode(Node_vec_t *neigh) { ??? n = ((neigh^).nel); { return n; }; } /* END NumberNeighborNode */ Quadv_vec_t *StarOfNode(Place_t @p, Topology_t *top)== PROCEDURE NodePoly(*top : Topology): REF ARRAY OF Quadp == /* Return four places such as the origins corresponding to nodes extremes of each tetrahedron of Triangulation. */ REF ARRAY OF Quadp poly = NEW(REF ARRAY OF Quadp, top->NP); { for (i = 0; i < top->NP; i++) { ??? da = top->region[i]; Place_t a0 = Tors(da); Place_t db = Clock(PrevE(da)); Place_t b0 = Tors(db); { /*assert(DegreeOfNode(a0) == 4); */ /* consuming time */ Place_t a1 = NextE(a0); Place_t a2 = NextE(a1); Place_t a3 = PrevE(b0); { if ((NextE(a2)!=a0)) { fprintf(stderr, "\nTriangulation: This topology isn't" \ " a triangulation\n"); assert(NextE(a2) == a0); } assert(Pneg(a0)->num == i); if ((Pneg(b0)!=NULL)) { assert(Pneg(a0) == Pneg(b0)); } assert(a0!=b0); poly[i] = Quadp{a0,a1,a2,a3}; } } } return poly; } NodePoly; poly1 : REF Quadv_vec_t = Quadv_vec_new(top->NP); n : uint = 0; { ??? poly = NodePoly(top), v == OrgV(a); { for(i = 0; i < poly.nel; i++) { if (((OrgV(poly[i][0]) == v)) || ((OrgV(poly[i][1]) == v) OR (OrgV(poly[i][2]) == v)) || ((OrgV(poly[i][3]) == v) )){ for (j = 0; j < 4; j++) { if ((OrgV(poly[i][j]) == v )){ Node_t uv = OrgV(poly[i][j]); Node_t dv = OrgV(NextE(poly[i][j])); Node_t pv = OrgV(PrevE(poly[i][j])); Node_t rv = OrgV(PrevE(PrevF(poly[i][j]))); { poly1[n] = Quadv{uv,dv,pv,rv}; n++; } } } } } Quadv_vec_trim(&poly1,n); } } /* END StarOfNode */ uint NumberPolyOfStar(Quadv_vec_t *qv) { ??? n = (qv->nel); { return n; }; } /* END NumberPolyOfStar */ r4_vec_t *ComputeAllNodeNormals( *Topology_t *top; Coords_t *c; )== { ??? rvn = r4_vec_new(top->NV); ??? vn = rvn^; { for (i = 0; i < top->NV; i++) { vn[i] = NodeNormal(top->out[i], c, top) } return rvn } } /* END ComputeAllNodeNormals */ r4_vec_t *ComputeAllEdgeNormals(Topology_t *top, Coords_t *c; )== { ??? rvn = r4_vec_new(top->NE); ??? vn = rvn^; { for (i = 0; i < top->NE; i++) { vn[i] = EdgeNormal(top->edge[i].pa, c); } return rvn } } /* END ComputeAllEdgeNormals */ r4_vec_t *ComputeAllWallNormals(Topology_t *top, Coords_t *c; )== { ??? rvn = r4_vec_new(top->NF); ??? vn = rvn^; { for (i = 0; i < top->NF; i++) { vn[i] = WallNormal(top->wall[i].pa, c); } return rvn } } /* END ComputeAllWallNormals */ r4_vec_t *ComputeAllCellNormals(Topology_t *top, Coords_t *c; )== { ??? rvn = r4_vec_new(top->NP); ??? vn = rvn^; { for (i = 0; i < top->NP; i++) { vn[i] = PolyNormal(Tors(top->region[i]), c); } return rvn } } /* END ComputeAllCellNormals */ CellTopology MakeCellTopology(Place_t @p) uint *ne, nv; b: Place_t; ptop: CellTopology; { ??? star = Octf.RenumberEdgesOutOfNodes((Place_vec_t){a})^; { /* Gather walls: */ ptop.NF = (star.nel); ptop.fRef = Place_vec_new(ptop.NF); ne = 0; for (i = 0; i < ptop.NF; i++) { ??? side = Octf.Dual(star[i]); { ptop.fRef[i] = side; b = side; do { ne++; Edge_t e = PEdge(b); { e.xmark = FALSE; OrgV(b).xmark = FALSE; } b = NextE(b) } while (b != side); } } /* Gather @{edge->?}s: */ assert(ne % 2 == 0); ptop.NE = ne DIV 2; ptop.eRef = Place_vec_new(ptop.NE); ne = 0; for (i = 0; i < ptop.NF; i++) { ??? side = Octf.Dual(star[i]); { b = side; do { Edge_t e = PEdge(b); { if (! e.xmark) { ptop.eRef[ne] = b; ne++; e.xmark = TRUE; } } b = NextE(b) } while (!( b == side)); } } /* Gather nodes: */ ptop.NV = 2 + ptop.NE - ptop.NF; ptop.vRef = Place_vec_new(ptop.NV); nv = 0; for (i = 0; i < ptop.NE; i++) { ??? eu = ptop.eRef[i]; Node_t u = OrgV(eu); Place_t ev = Clock(PrevF(eu)); Node_t v = OrgV(ev); { if (! u.xmark){ ptop.vRef[nv] = eu; nv++; u.xmark = TRUE; } if (! v.xmark){ ptop.vRef[nv] = ev; nv++; v.xmark = TRUE; } PEdge(eu).xmark = FALSE; }; } for (i = 0; i < ptop.NV; i++){ Node_t u = OrgV(ptop.vRef[i]); { u.xmark = FALSE; } } } return ptop; } /* END MakeCellTopology */ /* =============== INPUT/OUTPUT =========================== */ CONST BoolChars == Mis.BoolChars; AlphaChars == Mis.AlphaChars; TopCom_t ReadTopology(FILE *rd) VAR Topology_t *top; char *cmt; CHAR n; { /* Topology */ ReadHeader(rd,"topology","99-08-25"); cmt = Mis.ReadCommentsJS(rd, '|'); /* Element counts: */ Lex.Skip(rd, cs = AlphaChars); top->NV = fget_int(rd); Lex.Skip(rd, cs = AlphaChars); top->NE = fget_int(rd); Lex.Skip(rd, cs = AlphaChars); top->NF = fget_int(rd); Lex.Skip(rd, cs = AlphaChars); top->NP = fget_int(rd); Lex.Skip(rd, cs = AlphaChars); top->NFE = fget_int(rd); Lex.Skip(rd, cs = AlphaChars); top->der = fget_int(rd); Lex.Skip(rd, cs = AlphaChars); top->bdr = fget_int(rd); Lex.Skip(rd); ??? map = NEW(REF ARRAY OF Octf.wedge, top->NFE)^; { /* Create node records: */ top->node = Node_vec_new(top->NV); top->out = Place_vec_new(top->NV); for (i = 0; i < top->NV; i++) { top->node[i] = MakeNode(); top->node[i]->num = i; } /* Create @{edge->?}s records */ top->edge = NEW(REF ARRAY OF Octf.edge, top->NE); for (i = 0; i < top->NE; i++) { top->edge[i] = Octf.Make@{Edge->?}(); top->edge[i]->num = i; } /* Create wall records: */ top->wall = Wall_vec_new(top->NF); for (i = 0; i < top->NF; i++) { top->wall[i] = Octf.MakeWall(); top->wall[i]->num = i; } /* Create cells records: */ top->cell = NEW(REF ARRAY OF Cell, top->NP); top->region = Place_vec_new(top->NP); for (i = 0; i < top->NP; i++) { top->cell[i] = MakeCell(); top->cell[i]->num = i; } /* Create wedge records: */ top->wedge = NEW(REF ARRAY OF Octf_Place_t, top->NFE); for (i = 0; i < top->NFE; i++) { top->wedge[i] = MakeWedge(); PWedge(top->wedge[i])->num = i; map[i] = PWedge(top->wedge[i]); } /* Read @{edge->?} data: */ EVAL Mis.ReadCommentsJS(rd, '|'); for (j = 0; j < top->NE; j++) { Lex.Skip(rd); ??? ne = fget_int(rd); /* index to @{edge->?} */ ??? e == top->edge[ne]; ??? pa == Octf.ReadPlace(rd,map); { assert(ne == j); e->pa = pa } } Lex.Skip(rd); /* Read wall data: */ EVAL Mis.ReadCommentsJS(rd, '|'); for (j = 0; j < top->NF; j++) { Lex.Skip(rd); ??? nf = fget_int(rd); /* index to wall */ ??? f = top->wall[nf]; ??? pa = Octf.ReadPlace(rd,map); { assert(nf == j); f->pa = pa; } } Lex.Skip(rd); /* Read wedge data: */ EVAL Mis.ReadCommentsJS(rd, '|'); for (j = 0; j < top->NFE; j++) { Lex.Skip(rd); ??? nfe = fget_int(rd); /* index to wedge */ ??? fe = NARROW(PWedge(top->wedge[nfe]), wedge); { assert(nfe == j); assert(top->wedge[nfe].bits == 0); fe->num = nfe; Octf.ReadWedge(rd, fe, map); for (k = 0; k < 4; k++) { Lex.Skip(rd); n = Rd.GetChar(rd); if (n!='-') { Rd.UnGetChar(rd); ??? m = fget_int(rd); ??? pa = PlaceInWedge(fe, 2*k); ??? vf = fe.org[k]; { if ((Rd.GetChar(rd) == 'v')) { vf = top->node[m]; top->out[m] = pa; } else { vf = top->cell[m]; top->region[m] = pa; } } } } Lex.Skip(rd); ??? nf = fget_int(rd); { fe.wall = top->wall[nf]; } assert(Rd.GetChar(rd) == 'f'); Lex.Skip(rd); ??? ne = fget_int(rd); { fe.edge = top->edge[ne] } assert(Rd.GetChar(rd) == 'e'); } } ReadFooter(rd,"topology"); fprintf(stderr, "rebuilding node tables...\n"); RebuildNodeTables(top); return (TopCom_t){top, cmt}; } } /* END ReadTopology */ void RebuildNodeTables(Topology *top) Place_t *b; { for (i = 0; i < top->NE; i++){ ??? pa = top->edge[i].pa; { assert(PEdge(pa) == top->edge[i]); assert(top->edge[i]->num == i); top->edge[i].node[0] = OrgV(pa); top->edge[i].node[1] = OrgV(Clock(pa)) } } for (i = 0; i < top->NF; i++) { ??? a = top->wall[i].pa; ??? deg = Octf.DegreeOfWall(a); ??? vr = NEW(REF ARRAY OF Node_t, deg); { assert(PWall(a) == top->wall[i]); assert(top->wall[i]->num == i); top->wall[i].node = vr; b = a; for (j = 0; j < deg; j++) { vr[j] = Org(b); b = NextE(b); } } } for (i = 0; i < top->NP; i++) { ??? a = top->region[i]; ??? ptop = MakeCellTopology(a); ??? vr = NEW(REF ARRAY OF Node_t, ptop.NV); { assert(Org(a) == top->cell[i]); assert(top->cell[i]->num == i); top->cell[i].node = vr; for (j = 0; j < ptop.NV; j++) { vr[j] = Org(ptop.vRef[j]) } } } } /* END RebuildNodeTables */ PROCEDURE ReadMaterials( FILE *rd; *Topology_t *top; bool_t ro_te = FALSE; ) { /* Materials */ ReadHeader(rd,"materials","99-08-25"); /* Read node data Materials: */ EVAL Mis.ReadCommentsJS(rd, '|'); for (j = 0; j < top->NV; j++) { Lex.Skip(rd); ??? nv = fget_int(rd); /* index to node */ ??? v = top->node[nv]; { assert(nv == j); v->num = nv; Lex.Skip(rd); v->exists = Mis.ReadBool(rd); Lex.Skip(rd); v.fixed = Mis.ReadBool(rd); Lex.Skip(rd); ??? cc = v.color; { cc[0] = Lex.Real(rd); Lex.Skip(rd); cc[1] = Lex.Real(rd); Lex.Skip(rd); cc[2] = Lex.Real(rd); } Lex.Skip(rd); ??? tt = v.transp; { tt[0] = Lex.Real(rd); Lex.Skip(rd); tt[1] = Lex.Real(rd); Lex.Skip(rd); tt[2] = Lex.Real(rd); } Lex.Skip(rd); v.radius = Lex.Real(rd); Lex.Skip(rd); v.label = Rd.GetText(rd,2); } } Lex.Skip(rd); /* Read @{edge->?} data materials: */ EVAL Mis.ReadCommentsJS(rd, '|'); for (j = 0; j < top->NE; j++) { Lex.Skip(rd); ??? ne = fget_int(rd); /* index to @{edge->?} */ ??? e = top->edge[ne]; { assert(ne == j); e->num = ne; Lex.Skip(rd); e->exists = Mis.ReadBool(rd); Lex.Skip(rd); ??? cc = e.color; { cc[0] = Lex.Real(rd); Lex.Skip(rd); cc[1] = Lex.Real(rd); Lex.Skip(rd); cc[2] = Lex.Real(rd); } Lex.Skip(rd); ??? tt = e.transp; { tt[0] = Lex.Real(rd); Lex.Skip(rd); tt[1] = Lex.Real(rd); Lex.Skip(rd); tt[2] = Lex.Real(rd); } Lex.Skip(rd); e.radius = Lex.Real(rd); Lex.Skip(rd); e.degenerate = Mis.ReadBool(rd); Lex.Skip(rd); ??? n = Rd.GetChar(rd); { if (n!='-') { Rd.UnGetChar(rd); e->root = fget_int(rd); } else { e->root = -1; } } } } Lex.Skip(rd); /* Read wall data materials: */ EVAL Mis.ReadCommentsJS(rd, '|'); for (j = 0; j < top->NF; j++) { Lex.Skip(rd); ??? nf = fget_int(rd); /* index to wall */ ??? f = top->wall[nf]; { assert(nf == j); f->num = nf; Lex.Skip(rd); f->exists = Mis.ReadBool(rd); Lex.Skip(rd); ??? cc = f.color; { cc[0] = Lex.Real(rd); Lex.Skip(rd); cc[1] = Lex.Real(rd); Lex.Skip(rd); cc[2] = Lex.Real(rd); } Lex.Skip(rd); ??? tt = f.transp; { tt[0] = Lex.Real(rd); Lex.Skip(rd); tt[1] = Lex.Real(rd); Lex.Skip(rd); tt[2] = Lex.Real(rd); } Lex.Skip(rd); f.degenerate = Mis.ReadBool(rd); Lex.Skip(rd); ??? n = Rd.GetChar(rd); { if (n!='-') { Rd.UnGetChar(rd); f->root = fget_int(rd); } else { f->root = -1; } } } } Lex.Skip(rd); /* Read cell data materials: */ if (top->NP!=0) { EVAL Mis.ReadCommentsJS(rd, '|'); } for (j = 0; j < top->NP; j++) { Lex.Skip(rd); ??? np = fget_int(rd); /* index to cell */ ??? p = top->cell[np]; { assert(np == j); p->num = np; Lex.Skip(rd); p->exists = Mis.ReadBool(rd); Lex.Skip(rd); ??? cc = p.color; { cc[0] = Lex.Real(rd); Lex.Skip(rd); cc[1] = Lex.Real(rd); Lex.Skip(rd); cc[2] = Lex.Real(rd); } Lex.Skip(rd); ??? tt = p.transp; { tt[0] = Lex.Real(rd); Lex.Skip(rd); tt[1] = Lex.Real(rd); Lex.Skip(rd); tt[2] = Lex.Real(rd); } Lex.Skip(rd); p.degenerate = Mis.ReadBool(rd); if (ro_te) { Lex.Skip(rd); ??? n = Rd.GetChar(rd); { if (n!='-') { Rd.UnGetChar(rd); p->root = fget_int(rd); } else { p->root = -1; } } } } } ReadFooter(rd,"materials"); Rd.Close(rd); CheckOutAndRegion(top); } /* END ReadMaterials */ TopCom_t ReadToTaMa(char *name, bool_t ro_te /* DF FALSE */) /* Where ro_te meaning "root tetrahedron". */ VAR TopCom_t tc; { ??? ntp = name & ".tp"; ??? rtp = FileRd.Open(ntp); { fprintf(stderr, "reading " & ntp & "\n"); tc = ReadTopology(rtp); fprintf(stderr, "OK leio tp\n"); Rd.Close(rtp); } ??? nma = name & ".ma"; ??? rma = FileRd.Open(nma); { fprintf(stderr, "reading " & nma & "\n"); ReadMaterials(rma, tc.top, ro_te); Rd.Close(rma); } fprintf(stderr, "OK leio ma\n"); return tc; } /* END ReadToTaMa */ Coords_t *ReadState(char *name)== <* FATAL Rd.Failure, Thread.Alerted,FloatMode.Trap, Lex.Error, OSError.E); Coords_t c; char *cmt; uint nv; { ??? rs = FileRd.Open(name & ".st"); { /* Read Headers File Formats*/ ReadHeader(rs,"state","99-08-25"); /* Element counts: */ Lex.Skip(rs, cs = AlphaChars); nv = fget_int(rs); Lex.Skip(rs); cmt = Mis.ReadCommentsJS(rs, '|'); c = NEW(REF Coords_t, nv); /* Read node data state: */ for (j = 0; j < nv; j++) { Lex.Skip(rs); ??? nv = fget_int(rs); { ??? cv = c[nv]; { cv[0] = Lex.LongReal(rs); Lex.Skip(rs); cv[1] = Lex.LongReal(rs); Lex.Skip(rs); cv[2] = Lex.LongReal(rs); Lex.Skip(rs); cv[3] = Lex.LongReal(rs); } } } ReadFooter(rs,"state"); Rd.Close(rs); return c; } } /* END ReadState */ void WriteState(char *name, Topology_t *top, Coords_t *c, char *cmt /* DF "" */) <* FATAL OSError.E); { ??? st = FileWr.Open(name & ".st"); ??? vWidth = NumDigits(top->NV - 1); { void WriteCoord(x: double) { Wr.PutText(st, Fmt.Pad(Fmt.LongReal(x, Fmt.Style.Sci, prec = 3), 7)); } WriteCoord; void WritePoint(*c: r4_t) { WriteCoord(c[0]); Wr.PutText(st, " "); WriteCoord(c[1]); Wr.PutText(st, " "); WriteCoord(c[2]); Wr.PutText(st, " "); WriteCoord(c[3]); } WritePoint; { WriteHeader(st,"state","99-08-25"); Wr.PutText(st, "nodes "); Wr.PutText(st, Fmt.Int(top->NV) & "\n"); if ((! Text.Empty(cmt))) { Mis.WriteCommentsJS(st, cmt & "\n", '|') } Mis.WriteCommentsJS(st, "\nNode data:\n", '|'); for (i = 0; i < top->NV; i++) { ??? v = top->node[i]; { /* state */ Wr.PutText(st, Fmt.Pad(Fmt.Int(v->num), vWidth)); Wr.PutText(st, " "); WritePoint(c[v->num]); Wr.PutText(st, "\n"); } } } WriteFooter(st, "state"); fclose(st); } } /* END WriteState */ void WriteMaterials(char *name, Topology_t *top, char *cmt /* DF "" */, bool_t ro_te = FALSE; /* root tetrahedron */ ) <* FATAL Wr.Failure, Thread.Alerted, OSError.E); INTEGER *pWidth; { if (top->NP == 0){ pWidth = 2 }else{ pWidth = NumDigits(top->NP-1); } ??? ma = FileWr.Open(name & ".ma"); ??? vWidth = NumDigits(top->NV - 1); ??? fWidth = NumDigits(top->NF - 1); ??? eWidth = NumDigits(top->NE - 1); { void WriteIntensity(r: REAL) { Wr.PutText(ma, Fmt.Real(r, Fmt.Style.Fix, prec = 2)); } WriteIntensity; void WriteColor(*c: frgb_t) { WriteIntensity(c[0]); Wr.PutText(ma, " "); WriteIntensity(c[1]); Wr.PutText(ma, " "); WriteIntensity(c[2]); } WriteColor; void WriteRadius(r: REAL) { if (r == 0.02 ){ Wr.PutText(ma, "0.020"); } else { Wr.PutText(ma,Fmt.Real(r, prec = 4)); } } WriteRadius; void WriteLabel(char *label) { Wr.PutText(ma, label); } WriteLabel; { WriteHeader(ma,"materials","99-08-25"); if ((! Text.Empty(cmt))) { WriteCommentsJS(ma, cmt & "\n", '|') } ??? m = NumDigits(top->NP); { WriteCommentsJS(ma,"nodes " & Fmt.Pad(Fmt.Int(top->NV),m), '|'); WriteCommentsJS(ma,"@{edge->?}s " & Fmt.Pad(Fmt.Int(top->NE),m), '|'); WriteCommentsJS(ma,"walls " & Fmt.Pad(Fmt.Int(top->NF),m), '|'); WriteCommentsJS(ma,"cells " & Fmt.Int(top->NP), '|'); } WriteCommentsJS(ma, "\nNode data:\n", '|'); for (i = 0; i < top->NV; i++) { ??? v = top->node[i]; { /* materials */ Wr.PutText(ma, Fmt.Pad(Fmt.Int(v->num), vWidth)); Wr.PutText(ma, " "); Wr.PutText(ma, Fmt.Char(BoolChars[v->exists])); Wr.PutText(ma, " "); Wr.PutText(ma, Fmt.Char(BoolChars[v.fixed])); Wr.PutText(ma, " "); WriteColor(v.color); Wr.PutText(ma, " "); WriteColor(v.transp); Wr.PutText(ma, " "); WriteRadius(v.radius); Wr.PutText(ma, " "); WriteLabel(v.label); Wr.PutText(ma, "\n"); } } WriteCommentsJS(ma, "\n@{Edge->?} data:\n", '|'); for (i = 0; i < top->NE; i++) { ??? e = top->edge[i]; { /* materials */ Wr.PutText(ma, Fmt.Pad(Fmt.Int(e->num), eWidth)); Wr.PutText(ma, " "); Wr.PutText(ma, Fmt.Char(BoolChars[e->exists])); Wr.PutText(ma, " "); WriteColor(e.color); Wr.PutText(ma, " "); WriteColor(e.transp); Wr.PutText(ma, " "); WriteRadius(e.radius); Wr.PutText(ma, " "); Wr.PutText(ma, Fmt.Char(BoolChars[e.degenerate])); Wr.PutText(ma, " "); if (e->root == -1) { Wr.PutText(ma, " - "); } else { Wr.PutText(ma, Fmt.Pad(Fmt.Int(e->root), eWidth)); } Wr.PutText(ma, "\n"); } } WriteCommentsJS(ma, "\nWall data:\n", '|'); for (i = 0; i < top->NF; i++) { ??? f = top->wall[i]; { Wr.PutText(ma, Fmt.Pad(Fmt.Int(f->num), fWidth)); Wr.PutText(ma, " "); Wr.PutText(ma, Fmt.Char(BoolChars[f->exists])); Wr.PutText(ma, " "); WriteColor(f.color); Wr.PutText(ma, " "); WriteColor(f.transp); Wr.PutText(ma, " "); Wr.PutText(ma, Fmt.Char(BoolChars[f.degenerate])); Wr.PutText(ma, " "); if (f->root == -1) { Wr.PutText(ma, " - "); } else { Wr.PutText(ma, Fmt.Pad(Fmt.Int(f->root), fWidth)); } Wr.PutText(ma, "\n"); } } if (top->NP!=0) { WriteCommentsJS(ma, "\nCell data:\n", '|'); } for (i = 0; i < top->NP; i++) { ??? p = top->cell[i]; { Wr.PutText(ma, Fmt.Pad(Fmt.Int(p->num), pWidth)); Wr.PutText(ma, " "); Wr.PutText(ma, Fmt.Char(BoolChars[p->exists])); Wr.PutText(ma, " "); WriteColor(p.color); Wr.PutText(ma, " "); WriteColor(p.transp); Wr.PutText(ma, " "); Wr.PutText(ma, Fmt.Char(BoolChars[p.degenerate])); if (ro_te) { Wr.PutText(ma, " "); if (p->root == -1) { Wr.PutText(ma, " - "); } else { Wr.PutText(ma, Fmt.Pad(Fmt.Int(p->root), pWidth)); } } Wr.PutText(ma, "\n"); } } } WriteFooter(ma, "materials"); fclose(ma); } } /* END WriteMaterials */ void FindDegeneracies(Topology *top) VAR bool_t allTriangles, allTetrahedra; vei, vej: REF ARRAY OF INTEGER; vfi, vfj: REF ARRAY OF INTEGER; vpi, vpj: REF ARRAY OF INTEGER; { /* Clean the marks for the attribute degenerate */ for (i = 0; i < top->NE; i++) { ??? ei = top->edge[i]; { ei.degenerate = FALSE } } /* Now, find @{edge->?} degeneracies */ vei = NEW(REF ARRAY OF INTEGER, 2); vej = NEW(REF ARRAY OF INTEGER, 2); for (i = 0; i < top->NE; i++) { for (j = i+1; j < top->NE; j++) { ??? ei = NARROW(top->edge[i], @{Edge->?}); ??? ej = NARROW(top->edge[j], @{Edge->?}); { for (k = 0; k < 2; k++) { vei[k] = ei.node[k]->num; vej[k] = ej.node[k]->num } Mis.InsertionSort(1,vei); Mis.InsertionSort(1,vej); if ((vei[0] == vej[0]) && (vei[1] == vej[1])) { ei.degenerate = TRUE; ej.degenerate = TRUE; } } } } /* Clean the marks for the attribute degenerate */ allTriangles = TRUE; for (i = 0; i < top->NF; i++) { ??? fi = top->wall[i]; ??? di = ((fi.node^).nel); { fi.degenerate = (di!=3); allTriangles = allTriangles) && ((di == 3) } } if (! allTriangles){ return; } /* Now, find degenerate @places of triangles */ vfi = NEW(REF ARRAY OF INTEGER, 3); vfj = NEW(REF ARRAY OF INTEGER, 3); for (i = 0; i < top->NF; i++) { for (j = i+1; j < top->NF; j++) { ??? fi = top->wall[i]; ??? fj = top->wall[j]; { for (k = 0; k <= (2); k++) { vfi[k] = fi.node^[k]->num; vfj[k] = fj.node^[k]->num } Mis.InsertionSort(2,vfi); Mis.InsertionSort(2,vfj); if ((vfi[0] == vfj[0]) && (vfi[1] == vfj[1]) && (vfi[2] == vfj[2])) { fi.degenerate = TRUE; fj.degenerate = TRUE } } } } /* Clean the marks for the attribute degenerate */ allTetrahedra = TRUE; for (i = 0; i < top->NP; i++) { ??? pi = top->cell[i]; ??? di = ((pi.node^).nel); { pi.degenerate = (di!=4); allTetrahedra = allTetrahedra) && ((di == 4) } } if (! allTetrahedra){ return; } /* Now, find cell degeneracies */ vpi = NEW(REF ARRAY OF INTEGER, 4); vpj = NEW(REF ARRAY OF INTEGER, 4); for (i = 0; i < top->NP; i++) { for (j = i+1; j < top->NP; j++) { ??? pi = top->cell[i]; ??? pj = top->cell[j]; { for (k = 0; k < 4; k++) { vpi[k] = pi.node^[k]->num; vpj[k] = pj.node^[k]->num } Mis.InsertionSort(3,vpi); Mis.InsertionSort(3,vpj); if ((vpi[0]==vpj[0]) && (vpi[1]==vpj[1]) && (vpi[2]==vpj[2]) && (vpi[3]==vpj[3])) { pi.degenerate = TRUE; pj.degenerate = TRUE; } } } } } /* END FindDegeneracies */ void WriteTopology(char *name, Topology_t *top, char *cmt /* DF " " */) <* FATAL Wr.Failure, Thread.Alerted, OSError.E); INTEGER *pWidth; { if (top->NP == 0){ pWidth= 2 }else{ pWidth = NumDigits(top->NP-1); } ??? tp = FileWr.Open(name & ".tp"); ??? vWidth = NumDigits(top->NV - 1); ??? eWidth = NumDigits(top->NE - 1); ??? fWidth = NumDigits(top->NF - 1); ??? feWidth = NumDigits(top->NFE -1); { WriteHeader(tp,"topology","99-08-25"); if ((! Text.Empty(cmt))) { WriteCommentsJS(tp, cmt & "\n", '|') } ??? m = NumDigits(top->NFE); { Wr.PutText(tp, "nodes "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NV),m) & "\n"); Wr.PutText(tp, "@{edge->?}s "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NE),m) & "\n"); Wr.PutText(tp, "walls "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NF),m) & "\n"); Wr.PutText(tp, "cells "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NP),m) & "\n"); Wr.PutText(tp, "wedges "); Wr.PutText(tp, Fmt.Int(top->NFE) & "\n"); Wr.PutText(tp, "der "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->der),m) & "\n"); Wr.PutText(tp, "bdr "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->bdr),m) & "\n"); } WriteCommentsJS(tp, "\n@{Edge->?} data:\n", '|'); for (i = 0; i < top->NE; i++) { ??? e = top->edge[i]; { Wr.PutText(tp, Fmt.Pad(Fmt.Int(e->num), eWidth)); Wr.PutText(tp, " "); Octf.PrintPlace(tp, e.pa, eWidth+1, TRUE); } } WriteCommentsJS(tp, "\nWall data:\n", '|'); for (i = 0; i < top->NF; i++) { ??? f = top->wall[i]; { Wr.PutText(tp, Fmt.Pad(Fmt.Int(f->num), fWidth)); Wr.PutText(tp, " "); Octf.PrintPlace(tp, f.pa, fWidth+1, TRUE); } } WriteCommentsJS(tp, "\nwedge data:\n", '|'); for (i = 0; i < top->NFE; i++) { ??? fe = NARROW(PWedge(top->wedge[i]), wedge); { Wr.PutText(tp, Fmt.Pad(Fmt.Int(fe->num), feWidth)); Wr.PutText(tp, " "); Octf.PrintWedge(tp, fe, feWidth); Wr.PutText(tp, " "); for (j = 0; j < 4; j++) { ??? n = fe.org[j]; { TYPECASE n OF | NULL ==> for (i = 0; i <= (pWidth-2); i++){ Wr.PutText(tp," "); } Wr.PutText(tp, " - "); | Node(v) ==> Wr.PutText(tp, Fmt.Pad(Fmt.Int(v->num), vWidth) & "v "); | Cell(p) ==> Wr.PutText(tp, Fmt.Pad(Fmt.Int(p->num), pWidth) & "p "); } else { /* nothing */ } } } Wr.PutText(tp, " "); ??? f = fe.wall; { Wr.PutText(tp, Fmt.Pad(Fmt.Int(f->num), fWidth)); } Wr.PutText(tp, "f "); ??? e = fe.edge; { Wr.PutText(tp, Fmt.Pad(Fmt.Int(e->num), eWidth)); } Wr.PutText(tp, "e\n"); } } WriteFooter(tp, "topology"); fclose(tp); } } /* END WriteTopology */ void WriteDualTopology(char *name, Topology_t *top, char *cmt /* DF " " */) <* FATAL Wr.Failure, Thread.Alerted, OSError.E); void PrintDualPlace_t(FILE *wr, Place_t @p; feWidth: uint) { Wr.PutText(wr,Fmt.Pad(Fmt.Int(PWedge(a)->num), feWidth) & ":" \ Fmt.Int((Octf.SrotBits(a)+ 3) % 4) & ":" & Fmt.Int(Octf.SpinBit(a))); } PrintDualPlace_t; void PrintDualwedge(FILE *wr, Wedge_t w, feWidth: uint) Place_t *b; { b = Srot(PlaceInWedge(w, 0)); for (i = 0; i < 4; i++) { PrintDualPlace_t(wr, NextF(b), feWidth); Wr.PutText(wr, " "); b = Srot(b) } } PrintDualwedge; { ??? tp = FileWr.Open(name & ".tp"); ??? vWidth = NumDigits(MAX(1, top->NV - 1)); ??? eWidth = NumDigits(MAX(1,top->NE - 1)); ??? fWidth = NumDigits(MAX(1,top->NF - 1)); ??? pWidth = NumDigits(MAX(1,top->NP - 1)); ??? feWidth = NumDigits(MAX(1,top->NFE -1)); { WriteHeader(tp,"topology","99-08-25"); if ((! Text.Empty(cmt))) { WriteCommentsJS(tp, cmt & "\n", '|') } ??? m = NumDigits(top->NFE); { Wr.PutText(tp, "nodes "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NP),m) & "\n"); Wr.PutText(tp, "@{edge->?}s "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NF),m) & "\n"); Wr.PutText(tp, "walls "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NE),m) & "\n"); Wr.PutText(tp, "cells "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(top->NV),m) & "\n"); Wr.PutText(tp, "wedges "); Wr.PutText(tp, Fmt.Int(top->NFE) & "\n"); Wr.PutText(tp, "der "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(1),m) & "\n"); Wr.PutText(tp, "bdr "); Wr.PutText(tp, Fmt.Pad(Fmt.Int(1),m) & "\n"); } WriteCommentsJS(tp, "\n@{Edge->?} data:\n", '|'); for (i = 0; i < top->NF; i++) { ??? f = top->wall[i]; { Wr.PutText(tp, Fmt.Pad(Fmt.Int(f->num), fWidth)); Wr.PutText(tp, " "); Octf.PrintPlace(tp, f.pa, feWidth); /* [sic] */ Wr.PutText(tp, "\n") } } WriteCommentsJS(tp, "\nWall data:\n", '|'); for (i = 0; i < top->NE; i++) { ??? e = top->edge[i]; { Wr.PutText(tp, Fmt.Pad(Fmt.Int(e->num), eWidth)); Wr.PutText(tp, " "); Octf.PrintPlace(tp, e.pa, feWidth+1); /* [sic] */ Wr.PutText(tp, "\n") } } WriteCommentsJS(tp, "\nwedge data:\n", '|'); for (i = 0; i < top->NFE; i++) { ??? fe = PWedge(top->wedge[i]); { Wr.PutText(tp, Fmt.Pad(Fmt.Int(fe->num), feWidth)); Wr.PutText(tp, " "); PrintDualwedge(tp, fe, feWidth); Wr.PutText(tp, " "); for (j = 0; j < 4; j++) { ??? n = fe.org[(j+1) % 4]; { TYPECASE n OF | NULL ==> for (i = 0; i <= (pWidth-2); i++){ Wr.PutText(tp," "); } Wr.PutText(tp, " - "); | Node(v) ==> Wr.PutText(tp, Fmt.Pad(Fmt.Int(v->num), vWidth) & "p "); | Cell(p) ==> Wr.PutText(tp, Fmt.Pad(Fmt.Int(p->num), pWidth) & "v "); } else { /* nothing */ } } } Wr.PutText(tp, " "); ??? e = fe.edge; { Wr.PutText(tp, Fmt.Pad(Fmt.Int(e->num), eWidth)); } Wr.PutText(tp, "f "); ??? f = fe.wall; { Wr.PutText(tp, Fmt.Pad(Fmt.Int(f->num), fWidth)); } Wr.PutText(tp, "e"); Wr.PutText(tp, "\n"); } } WriteFooter(tp, "topology"); fclose(tp); } } /* END WriteDualTopology */ void WriteTable(char *name, Topology_t *top, char *cmt /* DF "" */, <*UNUSED); bool_t debug = FALSE; ) == <* FATAL Wr.Failure, Thread.Alerted, OSError.E); { RebuildNodeTables(top); ??? tab = FileWr.Open(name & ".tb"); ??? vWidth = NumDigits(MAX(1, top->NV - 1)); ??? eWidth = NumDigits(MAX(1, top->NE - 1)); ??? fWidth = NumDigits(MAX(1, top->NF - 1)); ??? pWidth = NumDigits(MAX(1, top->NP - 1)); { WriteHeader(tab,"table","99-08-25"); if ((! Text.Empty(cmt))) { Mis.WriteCommentsJS(tab, cmt & "\n", '|') } /* Now, writes @{edge->?} nodes */ Mis.WriteCommentsJS(tab,"\n@{Edge->?} data:\n", '|'); for (i = 0; i < top->NE; i++) { ??? e = top->edge[i]; { Wr.PutText(tab, Fmt.Pad(Fmt.Int(e->num), eWidth)); Wr.PutText(tab, " "); Wr.PutText(tab, Fmt.Pad(Fmt.Int(e.node[0]->num),eWidth) & "v "); Wr.PutText(tab, Fmt.Pad(Fmt.Int(e.node[1]->num),eWidth) & "v\n"); } } /* Now, writes wall nodes */ Mis.WriteCommentsJS(tab, "\nWall data:\n", '|'); for (i = 0; i < top->NF; i++) { ??? f = top->wall[i]; ??? fv = f.node^; { Wr.PutText(tab, Fmt.Pad(Fmt.Int(f->num), fWidth)); Wr.PutText(tab, " "); for (j = 0 TO (fv.nel-1)) { Wr.PutText(tab,Fmt.Pad(Fmt.Int(fv[j]->num), vWidth) & "v "); } Wr.PutText(tab, "\n"); } } /* Now, writes cell nodes */ if (top->NP!=0) { Mis.WriteCommentsJS(tab, "\nCell data:\n", '|'); } for (i = 0; i < top->NP; i++) { ??? p = top->cell[i]; ??? pv = p.node^; { Wr.PutText(tab, Fmt.Pad(Fmt.Int(p->num), pWidth)); Wr.PutText(tab, " "); for (j = 0 TO (pv.nel-1)) { Wr.PutText(tab,Fmt.Pad(Fmt.Int(pv[j]->num), vWidth) & "v "); } Wr.PutText(tab, "\n"); } } WriteFooter(tab, "table"); fclose(tab); } } /* END WriteTable */ void WriteStDe( FILE *wr, Coords_t *c, Coords_t *Dc; uint prec /* DF 4 */, char *cmt /* DF "" */) void WriteCoord(x: double) { Wr.PutText(wr, Fmt.LongReal(x, Fmt.Style.Sci, prec = prec)) } WriteCoord; void WritePoint(*p: r4_t) { WriteCoord(p[1]); Wr.PutText(wr, " "); WriteCoord(p[0]); Wr.PutText(wr, " "); WriteCoord(p[2]); Wr.PutText(wr, " "); WriteCoord(p[3]); Wr.PutText(wr, " "); } WritePoint; { int NV = (c.nel); ??? d = NumDigits(NV-1); { WriteHeader(wr,"state-derivatives","99-08-25"); Mis.WriteCommentsJS(wr, "\n" & cmt & "\n",'|'); Wr.PutText(wr, "nodes == " & Fmt.Int(NV) & "\n"); for (i = 0; i < NV; i++) { Wr.PutText(wr, Fmt.Pad(Fmt.Int(i), d) & ": "); WritePoint(c[i]); Wr.PutText(wr, " "); /* if the derivatives of nodes are zero them writes "0 0 0 0" else writes velocites with the format for write points defined here. */ if ((Dc[i][0] == 0.0) && (Dc[i][1] == 0.0) && (Dc[i][2] == 0.0 AND Dc[i][3] == 0.0 )){ Wr.PutText(wr, "0 0 0 0") } else { WritePoint(Dc[i]) } Wr.PutText(wr, "\n"); } WriteFooter(wr, "state-derivatives"); Wr.PutText(wr, "\n"); fflush(wr); } } /* END WriteStDe */ { } Triangulation. /* Copyright © 2000 Universidade Estadual de Campinas (UNICAMP) */ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // #FILE cur/liboct/JUNK/Mapi.c /* [!!! JUNK -- Superseded by Bary.c !!!] */ /* This module contain essentially procedures created by J. Stolfi and R. Marcone (see the copyright and authorship futher down), modified L. P. Lozada for the case tridimensional. */ // #INCLUDE // #INCLUDE // #INCLUDE FROM Octf IMPORT Spin; REVEAL double FacetEdge = PublicFacetEdge BRANDED OBJECT OVERRIDES init = MapInit; ; } FacetEdge MapInit(s: FacetEdge; uint order1, uint order2) { EVAL NARROW(s, Octf.FacetEdge).init(); s.ca = Triangulation.MakeTetraTopo(order1,order2); s.order1 = order1; s.order2 = order2; return s } /* END MapInit */ Pair MakeFacetEdge(uint order1, uint order2) { FacetEdge e = FacetEdge_new(); FacetEdge_init(order1,order2)){ return (Pair){/*facetedge*/ e, /*bits*/ 0}; ; } } /* END MakeFacetEdge */ Triangulation.Pair Corner(Pair a) { SRBits r = a.bits; { return a.facetedge.ca[r]; ; } } /* END Corner */ PROCEDURE SetCorner(Pair a, Triangulation.Pair c) { SRBits r = a.bits; { a.facetedge.ca[r] = c; ; } } /* END SetCorner */ Triangulation.Pair CCorner(Pair a) { SRBits r = (a.bits ^ 1u); { return Spin(a.facetedge.ca[r]); ; } } /* END CCorner */ PROCEDURE SetCCorner(Pair a, Triangulation.Pair c) { SRBits r = (a.bits ^ 1u); { a.facetedge.ca[r] = Spin(c); ; } } /* END SetCCorner */ { ;} Mapi. /* ***************** START OF COPYRIGHT) AND AND (AUTHORSHIP NOTICE ********** All files in this directory tree are Copyright 1996 by Jorge Stolfi, Rober Marcone Rosi, and Universidade Estadual de Campinas, Brazil--- unless stated otherwise in the files themselves. THESE FILES ARE DISTRIBUTED with (NO GUARANTEE OF ANY KIND. Neither the authors nor their employers may be held responsible for any losses or damages attributed to their use. These files may be freely copied, distributed, modified, and used for any purpose; provided that any subtantial excerpt of these files that is redistributed or incorporated in other software packages is accompanied by this copyright and authorship notice, and is made freely available under these same terms. ***************** ;} OF COPYRIGHT) AND AND (AUTHORSHIP NOTICE ************ */ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // #FILE cur/liboct/JUNK/Mapi.h #ifndef Mapi_H #define Mapi_H /* [!!! JUNK -- Superseded by Bary.h !!!] */ /* This interface contain essentially procedures created by J. Stolfi and R. Marcone (see the copyright and authorship futher down), modified by L. P. Lozada for the case tridimensional. */ // #INCLUDE // #INCLUDE TYPE double Pair = Octf_Pair; FacetEdge <: PublicFacetEdge; double PublicFacetEdge = Octf.FacetEdge OBJECT Pair ca[8]; order1, order2 : uint ; METHODS init(uint order1, uint order2): FacetEdge; /* Initializes "self" hangs on the edge component from it a new unglued topological tetrahedron of the given order: "order1xorder2". */ ; } Pair MakeFacetEdge(uint order1, uint order2); /* Returns Pair{NEW(FacetEdge).init(order1,order2), bits = 0} */ Octf_Pair Corner(Pair a); /* The pair (of triangulation) "c" belong to topological tetrahedron associated to pair (of mapi) "a", which has Org(c) == Org(a) and lies onto boundary of tetrahedron such as in one topological tetrahedron to perform: co[0] == Corner(a) co[1] == Corner(aSpin) co[2] == Corner(aSrot) co[3] == Corner(aSrotSpin) co[4] == Corner(aClock) co[5] == Corner(aClockSpin) co[6] == Corner(aTors) co[7] == Corner(aTorsSpin) */ Octf_Pair CCorner(Pair a); /* The pair (of triangulation) "c" belong to topological tetrahedron associated to pair (of mapi) "a", which has Org(c) == Org(a) and lies onto boundary of tetrahedron such as in one topological tetrahedron to perform: Spin(co[1]) == CCorner(a) Spin(co[0]) == CCorner(aSpin) Spin(co[3]) == CCorner(aSrot) Spin(co[2]) == CCorner(aSrotSpin) Spin(co[5]) == CCorner(aClock) Spin(co[4]) == CCorner(aClockSpin) Spin(co[7]) == CCorner(aTors) Spin(co[6]) == CCorner(aTorsSpin) */ void SetCorner(Pair a, Pair c); void SetCCorner(Pair a, Pair c); #endif ***************** START OF COPYRIGHT) AND AND (AUTHORSHIP NOTICE ********** All files in this directory tree are Copyright 1996 by Jorge Stolfi, Rober Marcone Rosi, and Universidade Estadual de Campinas, Brazil--- unless stated otherwise in the files themselves. THESE FILES ARE DISTRIBUTED with (NO GUARANTEE OF ANY KIND. Neither the authors nor their employers may be held responsible for any losses or damages attributed to their use. These files may be freely copied, distributed, modified, and used for any purpose; provided that any subtantial excerpt of these files that is redistributed or incorporated in other software packages is accompanied by this copyright and authorship notice, and is made freely available under these same terms. ***************** ;} OF COPYRIGHT) AND AND (AUTHORSHIP NOTICE ************ */ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // #FILE cur/liboct/JUNK/QuadEdge.c /* See {QuadEdge.h} */ // #INCLUDE /* [!!! JUNK !!!] */ // #INCLUDE #include // #INCLUDE // #INCLUDE FROM Octf IMPORT Fnext, Spin, Clock, Enext_1, SpliceEdges, SpliceFacets, MakeFacetEdge; /* ========================== Function's QuadEdge ============== */ Arc Flip(s : Arc) { if ((s.d == 0)){ return Arc{pair = Spin(Fnext(s.pair)), d = 0}; }else{ return Arc{pair = Clock(Spin(s.pair)), d = 1}; ; } } /* END Flip */ Arc Sym(s : Arc) { if ((s.d == 0)){ return Arc{pair = Clock(Fnext(s.pair)), d = 0}; }else{ return Arc{pair = Clock(Fnext(s.pair)), d = 1}; ; } } /* END Sym */ Arc Sym_(s : Arc) { if ((s.d == 0)){ return Arc{pair = Fnext(s.pair), d = 0}; }else{ return Arc{pair = Fnext(s.pair), d = 1}; ; } } /* END Sym_ */ Arc Onext(s : Arc) { if ((s.d == 0)){ return Arc{pair = Clock(Fnext(Enext_1(s.pair))), d = 0}; }else{ return Arc{pair = Enext_1(s.pair), d = 1}; ; } } /* END Onext */ Arc Dual(s : Arc) { return Arc{pair = s.pair, d = 1-s.d}; } /* END Dual */ Arc Rot(s : Arc) { return Dual(Flip(s)); } /* END Rot */ Arc Oprev(s : Arc) { return Rot(Onext(Rot(s))); } /* END Oprev */ Arc Lnext(s : Arc) { return Oprev(Sym(s)); } /* END Lnext */ Arc Rprev(s : Arc) { return Onext(Sym(s)); } /* END Rprev */ /* ============================= Operator's QuadEdge =============== */ VAR FacetEdgeCount : uint = 1; PROCEDURE MakeEdge() : Arc == <* FATAL Wr.Failure, Thread.Alerted); Pair *a,b; { a = MakeFacetEdge(); a.facetedge.num = FacetEdgeCount; INC(FacetEdgeCount); /*Octf.PrintPair(Stdio.stdout, a); Wr.PutText(Stdio.stdout, "\n"); */ b = MakeFacetEdge(); /*b.facetedge.num = FacetEdgeCount; INC(FacetEdgeCount); Octf.PrintPair(Stdio.stdout, b); Wr.PutText(Stdio.stdout, "\n"); */ SpliceFacets(a, b); SpliceEdges(a, Clock(b)); return Arc{pair = a, d = 0}; } /* END MakeEdge */ PROCEDURE MakeLoop() : Arc == <* FATAL Wr.Failure, Thread.Alerted); Pair *a,b; { a = MakeFacetEdge(); a.facetedge.num = FacetEdgeCount; INC(FacetEdgeCount); Octf.PrintPair(Stdio.stdout, a); Wr.PutText(Stdio.stdout, "\n"); b = MakeFacetEdge(); b.facetedge.num = FacetEdgeCount; INC(FacetEdgeCount); Octf.PrintPair(Stdio.stdout, b); Wr.PutText(Stdio.stdout, "\n"); SpliceFacets(a, b); return Arc{pair = a, d = 0}; } /* END MakeLoop */ void Splice(s1,s2 : Arc) == <* FATAL Wr.Failure, Thread.Alerted); { if (((s1.d == 0) AND AND (s2.d == 0) )){ Octf.PrintPair(Stdio.stdout, Spin(Clock(s1.pair)),1); Wr.PutText(Stdio.stdout, "\n"); Octf.PrintPair(Stdio.stdout, Spin(Clock(s2.pair)),1); Wr.PutText(Stdio.stdout, "\n"); SpliceEdges(Spin(Clock(s1.pair)), Spin(Clock(s2.pair))); }else{ SpliceEdges(Enext_1(s1.pair), Enext_1(s2.pair)); ; } } /* END Splice */ void Splice_(s1,s2 : Arc) == <* FATAL Wr.Failure, Thread.Alerted); { if (((s1.d == 0) AND AND (s2.d == 0) )){ Octf.PrintPair(Stdio.stdout, Spin(Clock(s1.pair)),10); Wr.PutText(Stdio.stdout, "\n"); Octf.PrintPair(Stdio.stdout, Spin(Clock(s2.pair)),10); Wr.PutText(Stdio.stdout, "\n"); SpliceEdges(Spin(Clock(s1.pair)), Spin(Clock(s2.pair))); }else{ SpliceEdges(Enext_1(s1.pair), Enext_1(s2.pair)); ; } } /* END Splice_ */ { ;} QuadEdge. // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // #FILE cur/liboct/JUNK/QuadEdge.h #ifndef QuadEdge_H #define QuadEdge_H /* [!!! JUNK !!!] */ // #INCLUDE TYPE double Pair = Octf_Pair; double Arc = RECORD pair : Pair; d: Octf.DBit; } /* Definition the scheme of representation the Edges in function the FacetEdge structure */ /* ========================== Function's QuadEdge ============== */ Arc PROCEDURE Flip(s : Arc); Arc PROCEDURE Sym(s : Arc); Arc PROCEDURE Sym_(s : Arc); Arc PROCEDURE Onext(s : Arc); Arc PROCEDURE Oprev(s : Arc); Arc PROCEDURE Dual(s : Arc); Arc PROCEDURE Rot(s : Arc); Arc PROCEDURE Lnext(s : Arc); Arc PROCEDURE Rprev(s : Arc); /* ============================= Operator's QuadEdge =============== */ Arc PROCEDURE MakeEdge(); Arc PROCEDURE MakeLoop(); void Splice(a,b : Arc); void Splice_(a,b : Arc); ;} QuadEdge. // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // #FILE cur/liboct/JUNK/Trunc.c q = r4_Cos(a,b); if ((! (-1.0 <= q) && (q <= 1.0))) { q = FLOAT(TRUNC(q), double); }; z = Math.acos(q); /* Aproximation */ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // #FILE cur/liboct/JUNK/VStar.c /* See {VStar.h} */ // #INCLUDE [!!! JUNK??? May have been superseded by Squared.c !!! ??? ] /* This module contain procedures to build several structures such as: n-gons, triangles and squares, and solid polyhedra: cube, ball, bigcube (3D array of cube) with procedures for the gluing two such bigcubes. */ #include // #INCLUDE // #INCLUDE FROM Octf IMPORT Enext_1, Clock, Fnext, Enext, SetFnext, SetEdgeAll, Fnext_1, Spin, SetFaceAll; FROM Triangulation IMPORT Pair, Org, MakeVertex, SetAllOrgs, OrgV, SetAllPneg, MakePolyhedron, Node, Vertex, Glue, Ppos, MakeTetraTopo; FROM Squared IMPORT MakeTetrahedron, MakeTriangle; /* Procedures to build V-stars */ PROCEDURE MakeTetrahedronTriang() : ARRAY [0..3] OF Pair { ??? t = MakeTetrahedron(); { SubdivideTetrahedron(t[0]); return t; ; } } /* END MakeTetrahedronTriang */ void SubdivideTetrahedron(Pair a) == /* Subdivide the tetrahedron Pneg(a) into four new tetrahedra through the insertion of a new medial vertex of type "VP", six new faces, and four edges. */ Node *y; { /* Create the medial vertex "VP" */ y = MakeVertex(); ??? v = NARROW(y, Vertex); { ){ v.exists= TRUE; v.label = "VP"; v.num = 0; ; } ??? b = Fnext_1(a); ??? c = Enext_1(b); ??? d = Fnext(c); ??? e = Enext_1(a); ??? f_ = Fnext_1(e); ??? g = Enext_1(f_); ??? h = Fnext(g); ??? fa = a.facetedge.face; ??? fb = b.facetedge.face; ??? fg = g.facetedge.face; ??? fh = h.facetedge.face; ??? ea = a.facetedge.edge; ??? ec = c.facetedge.edge; ??? ee = e.facetedge.edge; ??? eh = h.facetedge.edge; ??? eeb = Enext(b).facetedge.edge; ??? eea = Enext(a).facetedge.edge; ??? i = Enext(a); ??? j = Fnext_1(i); ??? k = Enext(b); ??? l = Fnext(k); ??? u = Org(a); ??? v = Org(Clock(a)); ??? w = Org(c); ??? x = Org(e); ??? f1 = MakeTriangle(); ??? f2 = MakeTriangle(); ??? f3 = MakeTriangle(); ??? f4 = MakeTriangle(); ??? f5 = MakeTriangle(); ??? f6 = MakeTriangle(); ??? q1 = MakePolyhedron(); ??? q2 = MakePolyhedron(); ??? q3 = MakePolyhedron(); ??? q4 = MakePolyhedron(); { /* save attributes for the original faces: "fa", "fb", "fg", "fh". */ fa.exists = TRUE; fa.color = frgb_t{1.00, 1.00, 0.20}; fa.transp = frgb_t{0.70, 0.70, 0.70}; fb.exists = TRUE; fb.color = frgb_t{1.00, 1.00, 0.20}; fb.transp = frgb_t{0.70, 0.70, 0.70}; fg.exists = TRUE; fg.color = frgb_t{1.00, 1.00, 0.20}; fg.transp = frgb_t{0.70, 0.70, 0.70}; fh.exists = TRUE; fh.color = frgb_t{1.00, 1.00, 0.20}; fh.transp = frgb_t{0.70, 0.70, 0.70}; ea.exists = TRUE; ea.color = frgb_t{0.0, 0.0, 0.0}; ea.transp = frgb_t{0.0, 0.0, 0.0}; ea.radius = 0.004; ec.exists = TRUE; ec.color = frgb_t{0.0, 0.0, 0.0}; ec.transp = frgb_t{0.0, 0.0, 0.0}; ec.radius = 0.004; ee.exists = TRUE; ee.color = frgb_t{0.0, 0.0, 0.0}; ee.transp = frgb_t{0.0, 0.0, 0.0}; ee.radius = 0.004; eh.exists = TRUE; eh.color = frgb_t{0.0, 0.0, 0.0}; eh.transp = frgb_t{0.0, 0.0, 0.0}; eh.radius = 0.004; eea.exists = TRUE; eea.color = frgb_t{0.0, 0.0, 0.0}; eea.transp = frgb_t{0.0, 0.0, 0.0}; eea.radius = 0.004; eeb.exists = TRUE; eeb.color = frgb_t{0.0, 0.0, 0.0}; eeb.transp = frgb_t{0.0, 0.0, 0.0}; eeb.radius = 0.004; /* insert f1 */ SetFnext(b,f1); SetFnext(f1,a); /* insert f2 */ SetFnext(c,f2); SetFnext(f2,d); /* insert f3 */ SetFnext(f_,f3); SetFnext(f3,e); /* set the relations among f1,f2 and f3 */ SetFnext(Clock(Enext(f2)),Enext_1(f1)); SetFnext(Enext_1(f1),Clock(Enext(f3))); SetFnext(Clock(Enext(f3)),Clock(Enext(f2))); SetEdgeAll(Enext_1(f1), Enext_1(f1).facetedge.edge); /* insert f4 */ SetFnext(j,f4); SetFnext(f4,i); /* insert f5 */ SetFnext(k,f5); SetFnext(f5,l); /* insert f6 */ SetFnext(g,f6); SetFnext(f6,h); /* set the internal relations along edge "yv" */ SetFnext(Enext_1(f5),Enext_1(f4)); SetFnext(Enext_1(f4),Clock(Enext(f1))); SetFnext(Clock(Enext(f1)), Enext_1(f5)); SetEdgeAll(Clock(Enext(f1)),Clock(Enext(f1)).facetedge.edge); /* set the internal relations along edge "wy" */ SetFnext(Enext(f5),Clock(Enext_1(f6))); SetFnext(Clock(Enext_1(f6)),Clock(Enext_1(f2))); SetFnext(Clock(Enext_1(f2)),Enext(f5)); SetEdgeAll(Enext(f5),Enext(f5).facetedge.edge); /* set the internal relations along edge "xy" */ SetFnext(Enext(f6), Enext(f4)); SetFnext(Enext(f4), Clock(Enext_1(f3))); SetFnext(Clock(Enext_1(f3)), Enext(f6)); SetEdgeAll(Enext(f4), Enext(f4).facetedge.edge); /* set the overall edge component */ SetEdgeAll(a, ea); SetEdgeAll(c, ec); SetEdgeAll(e, ee); SetEdgeAll(i, eea); SetEdgeAll(k, eeb); SetEdgeAll(h, eh); SetFaceAll(a, fa); SetFaceAll(b, fb); SetFaceAll(g, fg); SetFaceAll(h, fh); /* set the origins */ SetAllOrgs(a,u); SetAllOrgs(Clock(a),v); SetAllOrgs(c,w); SetAllOrgs(e,x); SetAllOrgs(Enext_1(f1),y); /* set the polyhedrons */ SetAllPneg(a,q1); SetAllPneg(Spin(b),q2); SetAllPneg(Spin(g),q3); SetAllPneg(Spin(f6),q4); ; } } /* END SubdivideTetrahedron */ PROCEDURE MakeOctahedronTriang() : ARRAY[0..7] OF Pair == /* Builds a triangulated octahedron, with eight tetrahedra, and return one facet-edge pair by original face. */ VAR a : ARRAY[0..7] OF ARRAY[0..7] OF Pair; b : ARRAY[0..7] OF Pair; { for (i = 0; i <= (7); i++){ a[i] = Triangulation.MakeTetraTopo(1,1); ; } /* first level gluing tetrahedra */ EVAL Glue(Spin(a[1][1]),a[0][0],1); EVAL Glue(Spin(a[2][1]),a[1][0],1); EVAL Glue(Spin(a[3][1]),a[2][0],1); EVAL Glue(Spin(a[0][1]),a[3][0],1); /* second level gluing tetrahedra */ EVAL Glue(Spin(a[0][3]),a[4][2],1); EVAL Glue(Spin(a[1][3]),a[5][2],1); EVAL Glue(Spin(a[2][3]),a[6][2],1); EVAL Glue(Spin(a[3][3]),a[7][2],1); /* gluing between levels */ EVAL Glue(Spin(a[5][1]),a[4][0],1); EVAL Glue(Spin(a[6][1]),a[5][0],1); EVAL Glue(Spin(a[7][1]),a[6][0],1); EVAL Glue(Spin(a[4][1]),a[7][0],1); b[0] = a[0][2]; b[1] = a[1][2]; b[2] = a[2][2]; b[3] = a[3][2]; b[4] = a[4][3]; b[5] = a[5][3]; b[6] = a[6][3]; b[7] = a[7][3]; return b; } /* END MakeOctahedronTriang */ PROCEDURE MakeIcosahedronTriang() : ARRAY[0..19] OF Pair == /* Builds a triangualted icosahedron, with twenty tetrahedra, and return one facet-edge pair by original face (i.e. 20). */ VAR a : ARRAY[0..19] OF ARRAY[0..7] OF Pair; b : ARRAY[0..19] OF Pair; /* this variable will be retuned */ { for (i = 0; i <= (19); i++){ a[i] = MakeTetraTopo(1,1); ; } /* inside the first level */ EVAL Glue(Spin(a[1][1]),a[0][0],1); EVAL Glue(Spin(a[2][1]),a[1][0],1); EVAL Glue(Spin(a[3][1]),a[2][0],1); EVAL Glue(Spin(a[4][1]),a[3][0],1); EVAL Glue(Spin(a[0][1]),a[4][0],1); /* between the first and second level */ EVAL Glue(Spin(a[0][3]),a[5][2],1); EVAL Glue(Spin(a[1][3]),a[6][2],1); EVAL Glue(Spin(a[2][3]),a[7][2],1); EVAL Glue(Spin(a[3][3]),a[8][2],1); EVAL Glue(Spin(a[4][3]),a[9][2],1); /* inside the second level */ EVAL Glue(Spin(Enext(a[10][1])), a[5][0], 1); EVAL Glue(Spin(a[6][1]), Enext(a[10][0]),1); EVAL Glue(Spin(Enext(a[11][1])), a[6][0], 1); EVAL Glue(Spin(a[7][1]), Enext(a[11][0]),1); EVAL Glue(Spin(Enext(a[12][1])), a[7][0], 1); EVAL Glue(Spin(a[8][1]), Enext(a[12][0]),1); EVAL Glue(Spin(Enext(a[13][1])), a[8][0], 1); EVAL Glue(Spin(a[9][1]), Enext(a[13][0]),1); EVAL Glue(Spin(Enext(a[14][1])), a[9][0], 1); EVAL Glue(Spin(a[5][1]), Enext(a[14][0]),1); /* between the second and third level */ EVAL Glue (Spin(a[10][3]), a[15][2],1); EVAL Glue (Spin(a[11][3]), a[16][2],1); EVAL Glue (Spin(a[12][3]), a[17][2],1); EVAL Glue (Spin(a[13][3]), a[18][2],1); EVAL Glue (Spin(a[14][3]), a[19][2],1); /* inside the third level */ EVAL Glue(Spin(a[16][1]),a[15][0],1); EVAL Glue(Spin(a[17][1]),a[16][0],1); EVAL Glue(Spin(a[18][1]),a[17][0],1); EVAL Glue(Spin(a[19][1]),a[18][0],1); EVAL Glue(Spin(a[15][1]),a[19][0],1); /* Rescue the pairs to be returned */ for (i = 0; i <= (4); i++){ b[i] = a[i,2]; } for (i = 5; i <= (9); i++){ b[i] = a[i,3]; } for (i = 10; i <= (14); i++){ b[i] = a[i,2]; } for (i = 15; i <= (19); i++){ b[i] = a[i,3]; } /* some assertions */ for (i = 0; i <= (4); i++){ assert(Fnext(b[i+5]) == Spin(b[i])); ; } for (i = 10; i <= (14); i++){ assert(Fnext(b[i+5]) == Spin(b[i])); ; } return b; } /* END MakeIcosahedronTriang */ PROCEDURE MakeDodecahedronTriang() : ARRAY[0..11] OF ARRAY [0..4] OF Pair == /* Builds a triangulated Dodecahedron, with 60 tetrahedra, trough the automatic gluing of tetrahedra. */ TYPE double Row4I = ARRAY[0..3] OF uint ; VAR cell4 : REF ARRAY OF Row4I; tetra : REF ARRAY OF ARRAY [0..3] OF Pair; uint cellnum; dode : ARRAY[0..11] OF ARRAY [0..4] OF Pair; PROCEDURE Gluing(Ti,Tj,Ci,Cj: uint) : Pair == /* Gluing the tetrahedra Ti with Tj through the free faces Ci and Cj respectively. */ { IF /* 1 */ double Ci = 0) AND AND (Cj == 0)){ ??? Oi = cell4[Ti,0]; with ( ??? Di = cell4[Ti,1]; double Oj = cell4[Tj,0]; ??? Dj = cell4[Tj,1] DO if (((Oi==Dj)) AND AND ((Di==Oj))){ EVAL Glue(Clock(tetra[Tj,Cj]), tetra[Ti,Ci], 1); return Clock(tetra[Tj,Cj]); ;} ;}; ELSIF /* 2 */ double Ci = 0) AND AND (Cj == 1)){ ??? Oi = cell4[Ti,0]; with ( ??? Di = cell4[Ti,1]; double Oj = cell4[Tj,0]; ??? Dj = cell4[Tj,1] DO if (((Oi == Oj)) AND AND ((Di == Dj))){ EVAL Glue(Spin(tetra[Tj,Cj]), tetra[Ti,Ci], 1); return Spin(tetra[Tj,Cj]); ;} ;}; ELSIF /* 3 */ double Ci = 0) AND AND (Cj == 2)){ ??? Oi = cell4[Ti,0]; with ( ??? Di = cell4[Ti,1]; ??? Oj = cell4[Tj,2]; ??? Dj = cell4[Tj,3]; { if (((Oi == Oj)) AND AND ((Di == Dj))){ EVAL Glue(Spin(tetra[Tj,Cj]), tetra[Ti,Ci], 1); return Spin(tetra[Tj,Cj]); ;} ;} ELSIF /* 4 */ double Ci = 0) AND AND (Cj == 3)){ ??? Oi = cell4[Ti,0]; with ( ??? Di = cell4[Ti,1]; ??? Oj = cell4[Tj,2]; ??? Dj = cell4[Tj,3]; { if (((Oi == Dj)) AND AND ((Di == Oj))){ EVAL Glue(Clock(tetra[Tj,Cj]), tetra[Ti,Ci], 1); return Clock(tetra[Tj,Cj]); ;} ; } ELSIF /* 5 */ double Ci = 1) AND AND (Cj == 0)){ ??? Oi = cell4[Ti,0]; with ( ??? Di = cell4[Ti,1]; ??? Oj = cell4[Tj,0]; ??? Dj = cell4[Tj,1]; { if (((Oi == Oj)) AND AND ((Di == Dj))){ EVAL Glue(Spin(tetra[Tj,Cj]), tetra[Ti,Ci], 1); return Spin(tetra[Tj,Cj]); ;} ; } ELSIF /* 6 */ double Ci = 1) AND AND (Cj == 1)){ ??? Oi = cell4[Ti,0]; with ( ??? Di = cell4[Ti,1]; ??? Oj = cell4[Tj,0]; ??? Dj = cell4[Tj,1]; { if (((Oi==Dj)) AND AND ((Di==Oj))){ EVAL Glue(Clock(tetra[Tj,Cj]), tetra[Ti,Ci], 1); return Clock(tetra[Tj,Cj]); ;} ; } ELSIF /* 7 */ double Ci = 1) AND AND (Cj == 2)){ ??? Oi = cell4[Ti,0]; with ( ??? Di = cell4[Ti,1]; ??? Oj = cell4[Tj,2]; ??? Dj = cell4[Tj,3]; { if (((Oi==Dj)) AND AND ((Di==Oj))){ EVAL Glue(Clock(tetra[Tj,Cj]), tetra[Ti,Ci], 1); return Clock(tetra[Tj,Cj]); ;} ; } ELSIF /* 8 */ double Ci = 1) AND AND (Cj == 3)){ ??? Oi = cell4[Ti,0]; with ( ??? Di = cell4[Ti,1]; ??? Oj = cell4[Tj,2]; ??? Dj = cell4[Tj,3]; { if (((Oi == Oj)) AND AND ((Di == Dj))){ EVAL Glue(Spin(tetra[Tj,Cj]), tetra[Ti,Ci], 1); return Spin(tetra[Tj,Cj]); ;} ; } ELSIF /* 9 */ double Ci = 2) AND AND (Cj == 0)){ ??? Oi = cell4[Ti,2]; with ( ??? Di = cell4[Ti,3]; ??? Oj = cell4[Tj,0]; ??? Dj = cell4[Tj,1]; { if (((Oi==Oj)) AND AND ((Di==Dj))){ EVAL Glue(Spin(tetra[Tj,Cj]), tetra[Ti,Ci], 1); return Spin(tetra[Tj,Cj]); ;} ; } ELSIF /* 10 */ double Ci = 2) AND AND (Cj == 1)){ ??? Oi = cell4[Ti,2]; with ( ??? Di = cell4[Ti,3]; ??? Oj = cell4[Tj,0]; ??? Dj = cell4[Tj,1]; { if (((Oi==Dj)) AND AND ((Di==Oj))){ EVAL Glue(Clock(tetra[Tj,Cj]), tetra[Ti,Ci], 1); return Clock(tetra[Tj,Cj]); ;} ; } ELSIF /* 11 */ double Ci = 2) AND AND (Cj == 2)){ ??? Oi = cell4[Ti,2]; ??? Di = cell4[Ti,3]; ??? Oj = cell4[Tj,2]; ??? Dj = cell4[Tj,3]; { if (((Oi==Dj)) AND AND ((Di==Oj))){ EVAL Glue(Clock(tetra[Tj,Cj]), tetra[Ti,Ci], 1); return Clock(tetra[Tj,Cj]); ;} ;}; ELSIF /* 12 */ double Ci = 2) AND AND (Cj == 3)){ ??? Oi = cell4[Ti,2]; ??? Di == cell4[Ti,3]; ??? Oj = cell4[Tj,2]; ??? Dj = cell4[Tj,3]; { if (((Oi==Oj)) AND AND ((Di==Dj))){ EVAL Glue(Spin(tetra[Tj,Cj]), tetra[Ti,Ci], 1); return Spin(tetra[Tj,Cj]); ;} ;}; ELSIF /* 13 */ double Ci = 3) AND AND (Cj == 0)){ ??? Oi = cell4[Ti,2]; with ( ??? Di = cell4[Ti,3]; ??? Oj = cell4[Tj,0]; ??? Dj = cell4[Tj,1]; { if (((Oi==Dj)) AND AND ((Di==Oj))){ EVAL Glue(Clock(tetra[Tj,Cj]), tetra[Ti,Ci], 1); return Clock(tetra[Tj,Cj]); ;} ; } ELSIF /* 14 */ double Ci = 3) AND AND (Cj == 1)){ ??? Oi = cell4[Ti,2]; with ( ??? Di = cell4[Ti,3]; ??? Oj = cell4[Tj,0]; ??? Dj = cell4[Tj,1]; { if (((Oi==Oj)) AND AND ((Di==Dj))){ EVAL Glue(Spin(tetra[Tj,Cj]), tetra[Ti,Ci], 1); return Spin(tetra[Tj,Cj]); ;} ; } ELSIF /* 15 */ double Ci = 3) AND AND (Cj == 2)){ ??? Oi = cell4[Ti,2]; ??? Di == cell4[Ti,3]; ??? Oj = cell4[Tj,2]; ??? Dj = cell4[Tj,3]; { if (((Oi==Oj)) AND AND ((Di==Dj))){ EVAL Glue(Spin(tetra[Tj,Cj]), tetra[Ti,Ci], 1); return Spin(tetra[Tj,Cj]); ;} ;} ELSIF /* 16 */ double Ci = 3) AND AND (Cj == 3)){ ??? Oi = cell4[Ti,2]; ??? Di == cell4[Ti,3]; ??? Oj = cell4[Tj,2]; ??? Dj = cell4[Tj,3]; { if (((Oi==Dj)) AND AND ((Di==Oj))){ EVAL Glue(Clock(tetra[Tj,Cj]), tetra[Ti,Ci], 1); return Clock(tetra[Tj,Cj]); ;} ;} ; } return tetra[Ti,0]; } /* END Gluing */ PROCEDURE SetCornersTetra(Ti: uint ; row: Row4I) == /* Set the labels "row" in the tetrahedron Ti. */ { ??? a = OrgV(tetra[Ti,0]); with ( double b = OrgV(Clock(tetra[Ti,0])); double c = OrgV(Enext_1(tetra[Ti,1])); double d = OrgV(Enext_1(tetra[Ti,0])) ){ a.num = row[0]; b.num = row[1]; c.num = row[2]; d.num = row[3]; ; } ;} SetCornersTetra; PROCEDURE SaveFreeCorners() == /* Save the free corners on the triangulated dodecahedron, such allow us to use this polyhedron as the unglued scheme for the 3D maps obtained by the gluing of opposite faces on the dode- cahedron. */ PROCEDURE AssignmentOnFace(b: Pair; i: uint) == { if ((Octf.SpinBit(b) == 0)){ dode[i,0] = b; for (j = 1; j <= (4); j++){ dode[i,j] = Clock(Enext(Fnext(Enext_1(dode[i,j-1])))); ; } }else if ((Octf.SpinBit(b) == 1)){ dode[i,0] = b; for (j = 1; j <= (4); j++){ dode[i,j] = Clock(Enext_1(Fnext(Enext(dode[i,j-1])))); ;} ;} ;} AssignmentOnFace; VAR m,n: uint = 0; { for (i = 0; i <= (59); i++){ for (j = 0; j <= (3); j++){ ??? a = tetra[i,j]; with ( double d = Octf.DegreeFaceRing(a) ){ if ((d!=1) AND AND (Ppos(a) == NULL)){ dode[m,n] = tetra[i,j]; INC(n); if ((n == 5)){ INC(m); n = 0;; } ;} ;} ;} ; } ??? a = dode[0,0]; { AssignmentOnFace(a,0);}; dode[1,0] = Spin(Fnext(dode[0,0])); AssignmentOnFace(dode[1,0],1); dode[2,0] = Spin(Fnext(dode[0,4])); AssignmentOnFace(dode[2,0],2); dode[3,0] = Spin(Fnext(dode[0,3])); AssignmentOnFace(dode[3,0],3); dode[4,0] = Spin(Fnext(dode[0,2])); AssignmentOnFace(dode[4,0],4); dode[5,0] = Spin(Fnext(dode[0,1])); AssignmentOnFace(dode[5,0],5); dode[8, 0] = Spin(Fnext(dode[4,2])); AssignmentOnFace(dode[8,0],8); dode[6, 0] = Spin(Fnext(dode[8,3])); AssignmentOnFace(dode[6,0],6); dode[7, 0] = Spin(Fnext(dode[6,1])); AssignmentOnFace(dode[7,0], 7); dode[9, 0] = Spin(Fnext(dode[6,4])); AssignmentOnFace(dode[9,0], 9); dode[10,0] = Spin(Fnext(dode[6,3])); AssignmentOnFace(dode[10,0],10); dode[11,0] = Spin(Fnext(dode[6,2])); AssignmentOnFace(dode[11,0],11); /* the free corners are computed follows the scheme bellow: /|\ /|\ / | \ / | \ / | \ / | \ / 3 | 2 \ / 3 | 2 \ / | \ / | \ /-----|-----\ /-----|-----\ \ 4 / \ 1 / \ 4 / \ 1 / \ / 0 \ / \ / 0 \ / \/_____\/ \/_____\/ <--- ---> double SpinBit = 0 SpinBit == 1 */ ;} SaveFreeCorners; PROCEDURE MustBeGlue(Ti,Tj: Pair) : bool_t == /* Return TRUE if the faces "Ti.facetedge.face" and "Tj.facetedge.face" have coherent orientations and must be glued. */ { ??? a = OrgV(Ti).num; ??? ae = OrgV(Enext(Ti)).num; ??? ae_1 = OrgV(Enext_1(Ti)).num; ??? b = OrgV(Tj).num; ??? be = OrgV(Enext(Tj)).num; ??? be_1 = OrgV(Enext_1(Tj)).num; { if (((a == b) AND AND (ae == be) AND AND (ae_1 == be_1) OR (a == b) AND AND (ae == be_1) AND AND (ae_1 == be))){ return TRUE; } return FALSE ; } ;} MustBeGlue; PROCEDURE EnextK(Ti: Pair; k : uint) : Pair == /* Given a pair "Ti", this procedure return Enext^{k}(Ti). */ { if ((k == 0)){ return Ti }else if ((k == 1)){ return Enext(Ti) }else if ((k == 2)){ return Enext(Enext(Ti)) ; } return Ti; ;} EnextK; VAR poly : REF ARRAY OF ARRAY [0..7] OF Pair; count : uint = 1; faces : REF Pair_vec_t; glues : REF ARRAY OF Row4I; { cellnum = 60; cell4 = NEW(REF ARRAY OF Row4I, cellnum); poly = NEW(REF ARRAY OF ARRAY [0..7] OF Pair,cellnum); tetra = NEW(REF ARRAY OF ARRAY [0..3] OF Pair,cellnum); faces = Pair_vec_new( 4*cellnum); glues = NEW(REF ARRAY OF Row4I, 2*cellnum); /* creating topological tetrahedra */ for (i = 0; i < cellnum; i++){ poly[i] = MakeTetraTopo(1,1); ; } /* creating the tetrahedra */ for (i = 0; i < cellnum; i++){ for (j = 0; j <= (3); j++){ tetra[i,j] = poly[i,j]; assert(Ppos(tetra[i,j]) == NULL); ;} ; } /* cells with corners perfectly assigments */ cell4[ 0]=Row4I{100,500, 0, 1}; cell4[ 1]=Row4I{100,500, 1, 4}; cell4[ 2]=Row4I{100,500, 4, 7}; cell4[ 3]=Row4I{100,500, 7, 2}; cell4[ 4]=Row4I{100,500, 2, 0}; cell4[ 5]=Row4I{500,101, 0, 1}; cell4[ 6]=Row4I{500,101, 1, 5}; cell4[ 7]=Row4I{500,101, 5, 8}; cell4[ 8]=Row4I{500,101, 8, 3}; cell4[ 9]=Row4I{500,101, 3, 0}; cell4[ 10]=Row4I{500,102, 2, 0}; cell4[ 11]=Row4I{500,102, 0, 3}; cell4[ 12]=Row4I{500,102, 3, 9}; cell4[ 13]=Row4I{500,102, 9, 6}; cell4[ 14]=Row4I{500,102, 6, 2}; cell4[ 15]=Row4I{103,500, 6, 2}; cell4[ 16]=Row4I{500,103, 7, 2}; cell4[ 17]=Row4I{103,500, 7, 13}; cell4[ 18]=Row4I{103,500, 13, 12}; cell4[ 19]=Row4I{500,103, 6, 12}; cell4[ 20]=Row4I{500,104, 1, 4}; cell4[ 21]=Row4I{500,104, 4, 10}; cell4[ 22]=Row4I{500,104, 10, 11}; cell4[ 23]=Row4I{500,104, 11, 5}; cell4[ 24]=Row4I{500,104, 5, 1}; cell4[ 25]=Row4I{500,105, 9, 3}; cell4[ 26]=Row4I{500,105, 3, 8}; cell4[ 27]=Row4I{500,105, 8, 14}; cell4[ 28]=Row4I{500,105, 14, 15}; cell4[ 29]=Row4I{500,105, 15, 9}; cell4[ 30]=Row4I{500,106, 12, 6}; cell4[ 31]=Row4I{500,106, 6, 9}; cell4[ 32]=Row4I{500,106, 9, 15}; cell4[ 33]=Row4I{500,106, 15, 18}; cell4[ 34]=Row4I{500,106, 18, 12}; cell4[ 35]=Row4I{500,107, 4, 7}; cell4[ 36]=Row4I{500,107, 7, 13}; cell4[ 37]=Row4I{500,107, 13, 16}; cell4[ 38]=Row4I{500,107, 16, 10}; cell4[ 39]=Row4I{500,107, 10, 4}; cell4[ 40]=Row4I{108,500, 5, 8}; cell4[ 41]=Row4I{108,500, 8, 14}; cell4[ 42]=Row4I{108,500, 14, 17}; cell4[ 43]=Row4I{108,500, 17, 11}; cell4[ 44]=Row4I{108,500, 11, 5}; cell4[ 45]=Row4I{109,500, 18, 12}; cell4[ 46]=Row4I{109,500, 12, 13}; cell4[ 47]=Row4I{109,500, 13, 16}; cell4[ 48]=Row4I{109,500, 16, 19}; cell4[ 49]=Row4I{109,500, 19, 18}; cell4[ 50]=Row4I{500,110, 18, 15}; cell4[ 51]=Row4I{500,110, 15, 14}; cell4[ 52]=Row4I{500,110, 14, 17}; cell4[ 53]=Row4I{500,110, 17, 19}; cell4[ 54]=Row4I{500,110, 19, 18}; cell4[ 55]=Row4I{111,500, 10, 11}; cell4[ 56]=Row4I{111,500, 11, 17}; cell4[ 57]=Row4I{111,500, 17, 19}; cell4[ 58]=Row4I{111,500, 19, 16}; cell4[ 59]=Row4I{111,500, 16, 10}; /* set the labels for each tetrahedra */ for (i = 0; i < cellnum; i++){ SetCornersTetra(i,cell4[i]); ; } /* builds the table of faces for choose which tetrahedra must be gluing. */ for (i = 0; i < cellnum; i++){ for (k = 0; k <= (3); k++){ faces[(4*i)+k] = tetra[i,k]; ;} ; } /* computing which cells must be gluing. */ for (k = 0 TO LAST(faces^)){ for (l = k+1 TO LAST(faces^)){ for (m = 0; m <= (2); m++){ ??? e = EnextK(faces[l],m); with ( ){ if ((MustBeGlue(faces[k],e))){ ??? kc = k MOD 4; ??? kt = k DIV 4; ??? lc = l MOD 4; ??? lt = l DIV 4; { glues[count-1] = Row4I{kt,lt,kc,lc}; INC(count); ;} ;} ;} ;} ;} ; } /* Do the automatic gluing of tetrahedra */ for (i = 0 TO LAST(glues^)){ ??? c = glues[i]; { if ((c!=Row4I{0,0,0,0})){ EVAL Gluing(c[0],c[1],c[2],c[3]); ;} ;} ; } /* setting the origins. */ for (i = 0; i < cellnum; i++){ for (j = 0; j <= (3); j++){ ??? a = tetra[i,j]; with ( double b = Enext(a); double c = Enext_1(a) ){ Triangulation.SetAllOrgs(a,OrgV(a)); Triangulation.SetAllOrgs(b,OrgV(b)); Triangulation.SetAllOrgs(c,OrgV(c)); ;} ;} ; } SaveFreeCorners(); return dode; ;} MakeDodecahedronTriang; { ;} VStar. /**************************************************************************/ /* */ /* Copyright © 2000 Universidade Estadual de Campinas (UNICAMP) */ /* */ /* Authors: */ /* L. Lozada & J. Stolfi - UNICAMP */ /* */ /* This file can be freely used, distributed, and modified, provided that */ /* this copyright and authorship notice is included in every copy or */ /* derived version. */ /* */ /* DISCLAIMER: This software is offered ``as is'', without any guarantee */ /* as to fitness for any particular purpose. Neither the copyright */ /* holder nor the authors or their employers can be held responsible */ /* for any damages that may result from its use. */ /* */ /**************************************************************************/ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // #FILE cur/liboct/JUNK/VStar.h #ifndef VStar_H #define VStar_H [!!! JUNK ??? May have been superseded by Squared.h ??? !!!] /* Procedures to build several faces such as: n-gons, triangles and squares, and complexes : cube, ball, bigcube (3D array of cube) with procedures for the glueing of two such complexes. Created by L. Lozada (see the copyright and authorship futher down). Revisions: 30-08-2000 : Nice version of the "MakeOctahedron" procedure. 19-09-2000 : Added the procedure "MakeDodecahedronTriang". 27-10-2000 : Modified "SetCubePropiertes" procedure. 04-11-2000 : Added the "CubeNegVertices" and "CubeBarycenter" procedures. */ // #INCLUDE FROM Triangulation IMPORT Pair; PROCEDURE MakeTetrahedronTriang(): ARRAY [0..3] OF Pair; /* Builds a Tetrahedral vertex star. */ PROCEDURE MakeOctahedronTriang(): ARRAY [0..7] OF Pair; /* Builds a triangulated Octahedron. If Original==TRUE the procedure empha- size the original elements. */ PROCEDURE MakeIcosahedronTriang(): ARRAY [0..19] OF Pair; /* Builds a triangulated Icosahedron. If Original==TRUE the procedure empha- size the original elements. */ PROCEDURE MakeDodecahedronTriang() : ARRAY[0..11] OF ARRAY [0..4] OF Pair; /* Builds a triangulated Dodecahedron, trough the automatic gluing of tetrahedra. If Original==TRUE the procedure emphasize the original elements. */ ;} VStar. /**************************************************************************/ /* */ /* Copyright © 2000 Universidade Estadual de Campinas (UNICAMP) */ /* */ /* Authors: */ /* L. Lozada & J. Stolfi - UNICAMP */ /* */ /* This file can be freely used, distributed, and modified, provided that */ /* this copyright and authorship notice is included in every copy or */ /* derived version. */ /* */ /* DISCLAIMER: This software is offered ``as is'', without any guarantee */ /* as to fitness for any particular purpose. Neither the copyright */ /* holder nor the authors or their employers can be held responsible */ /* for any damages that may result from its use. */ /* */ /**************************************************************************/ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // #FILE cur/progs-more/JUNK/DeterLR3.c /* See {DeterLR3.h} */ // #INCLUDE // #INCLUDE #include #include VAR a,b,c: r3_t; double d; { a = (r3_t){ 2.0, 3.0, 0.00}; b = (r3_t){ 1.0, 0.0, 2.00}; c = (r3_t){ 1.0, 0.0, 0.00}; d = LR3Extras.Det(a,b,c); Wr.PutText(Stdio.stderr, "Determinante: " \ Fmt.LongReal(d, Fmt.Style.Fix,prec = 2) & "\n"); } DeterLR3. // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~