Index: /issm/trunk/src/c/Bamgx/Mesh2.h =================================================================== --- /issm/trunk/src/c/Bamgx/Mesh2.h (revision 2991) +++ /issm/trunk/src/c/Bamgx/Mesh2.h (revision 2992) @@ -9,6 +9,4 @@ #include "meshtype.h" #include "Metric.h" - -using namespace std; namespace bamg { @@ -808,29 +806,6 @@ }; /*}}}1*/ - - inline void Vertex::Set(const Vertex & rec,const Triangles & ,Triangles & ){ - *this = rec; - } - inline void GeometricalVertex::Set(const GeometricalVertex & rec,const Geometry & ,const Geometry & ){ - *this = rec; - } - inline void VertexOnVertex::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) { - *this = Th.VertexOnBThVertex[i]; - v = ThNew.vertices + Th.Number(v); - - } - Int4 inline Vertex::Optim(int i,int koption){ - Int4 ret=0; - if ( t && (vint >= 0 ) && (vint <3) ) - { - ret = t->Optim(vint,koption); - if(!i) - { - t =0; // for no future optime - vint= 0; } - } - return ret; - } - + + /*INLINE functions{{{1*/ // to sort in descending order template inline void HeapSort(T *c,long n){ @@ -858,19 +833,168 @@ } } - - //From Metric.cpp inline Real8 det3x3(Real8 A[3] ,Real8 B[3],Real8 C[3]){ return A[0]*( B[1]*C[2]-B[2]*C[1]) - - A[1]*( B[0]*C[2]-B[2]*C[0]) - + A[2]*( B[0]*C[1]-B[1]*C[0]); + - A[1]*( B[0]*C[2]-B[2]*C[0]) + + A[2]*( B[0]*C[1]-B[1]*C[0]); } - + inline TriangleAdjacent Adj(const TriangleAdjacent & a) + { return a.Adj();} + + inline TriangleAdjacent Next(const TriangleAdjacent & ta) + { return TriangleAdjacent(ta.t,NextEdge[ta.a]);} + + inline TriangleAdjacent Previous(const TriangleAdjacent & ta) + { return TriangleAdjacent(ta.t,PreviousEdge[ta.a]);} + inline void Adj(GeometricalEdge * & on,int &i) + {int j=i;i=on->DirAdj[i];on=on->Adj[j];} + inline Real4 qualite(const Vertex &va,const Vertex &vb,const Vertex &vc) + { + Real4 ret; + I2 ia=va,ib=vb,ic=vc; + I2 ab=ib-ia,bc=ic-ib,ac=ic-ia; + Icoor2 deta=Det(ab,ac); + if (deta <=0) ret = -1; + else { + Real8 a = sqrt((Real8) (ac,ac)), + b = sqrt((Real8) (bc,bc)), + c = sqrt((Real8) (ab,ab)), + p = a+b+c; + Real8 h= Max(Max(a,b),c),ro=deta/p; + ret = ro/h;} + return ret; + } + Icoor2 inline det(const Vertex & a,const Vertex & b,const Vertex & c){ + register Icoor2 bax = b.i.x - a.i.x ,bay = b.i.y - a.i.y; + register Icoor2 cax = c.i.x - a.i.x ,cay = c.i.y - a.i.y; + return bax*cay - bay*cax; + } + inline TriangleAdjacent FindTriangleAdjacent(Edge &E){ + Vertex * a = E.v[0]; + Vertex * b = E.v[1]; + + Triangle * t = a->t; + int i = a->vint; + TriangleAdjacent ta(t,EdgesVertexTriangle[i][0]); // Previous edge + if (!t || i<0 || i>=3){ + throw ErrorException(__FUNCT__,exprintf("!t || i<0 !! i>=3")); + } + if ( a!=(*t)(i)){ + throw ErrorException(__FUNCT__,exprintf("a!=(*t)(i)")); + } + int k=0; + do { // turn around vertex in direct sens (trigo) + k++; + if (k>=20000){ + throw ErrorException(__FUNCT__,exprintf("k>=20000")); + } + // in no crack => ta.EdgeVertex(1) == a otherwise ??? + if (ta.EdgeVertex(1) == a && ta.EdgeVertex(0) == b) return ta; // find + ta = ta.Adj(); + if (ta.EdgeVertex(0) == a && ta.EdgeVertex(1) == b) return ta; // find + --ta; + } while (t != (Triangle *)ta); + throw ErrorException(__FUNCT__,exprintf("FindTriangleAdjacent: triangle not found")); + return TriangleAdjacent(0,0);//for compiler + } + inline Vertex* TheVertex(Vertex * a){ + // give a unique vertex with smallest number + // in case on crack in mesh + Vertex * r(a), *rr; + Triangle * t = a->t; + int i = a->vint; + TriangleAdjacent ta(t,EdgesVertexTriangle[i][0]); // Previous edge + if (!t || i<0 || i>=3){ + throw ErrorException(__FUNCT__,exprintf("!t || i<0 !! i>=3")); + } + if ( a!=(*t)(i)){ + throw ErrorException(__FUNCT__,exprintf("a!=(*t)(i)")); + } + int k=0; + do { // turn around vertex in direct sens (trigo) + k++; + if (k>=20000){ + throw ErrorException(__FUNCT__,exprintf("k>=20000")); + } + // in no crack => ta.EdgeVertex(1) == a + if ((rr=ta.EdgeVertex(0)) < r) r = rr; + ta = ta.Adj(); + if ((rr=ta.EdgeVertex(1)) < r) r =rr; + --ta; + } while (t != (Triangle*) ta); + return r; + } + + /*}}}1*/ + + /*INLINE functions of CLASS Vertex{{{1*/ + inline void Vertex::Set(const Vertex & rec,const Triangles & ,Triangles & ){ + *this = rec; + } + Int4 inline Vertex::Optim(int i,int koption){ + Int4 ret=0; + if ( t && (vint >= 0 ) && (vint <3) ) + { + ret = t->Optim(vint,koption); + if(!i) + { + t =0; // for no future optime + vint= 0; } + } + return ret; + } + /*}}}1*/ + /*INLINE functions of CLASS GeometricalVertex{{{1*/ + inline void GeometricalVertex::Set(const GeometricalVertex & rec,const Geometry & ,const Geometry & ){ + *this = rec; + } + /*}}}1*/ + /*INLINE functions of CLASS VertexOnVertex{{{1*/ + inline void VertexOnVertex::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) { + *this = Th.VertexOnBThVertex[i]; + v = ThNew.vertices + Th.Number(v); + + } +/*}}}1*/ + /*INLINE functions of CLASS Triangles{{{1*/ inline Triangles::Triangles(Int4 i) :Gh(*new Geometry()),BTh(*this){PreInit(i);} - - TriangleAdjacent CloseBoundaryEdge(I2 ,Triangle *, double &,double &) ; - TriangleAdjacent CloseBoundaryEdgeV2(I2 A,Triangle *t, double &a,double &b); - - Int4 FindTriangle(Triangles &Th, Real8 x, Real8 y, double* a,int & inside); - + inline void Triangles::ReMakeTriangleContainingTheVertex() + { + register Int4 i; + for ( i=0;i=Th.nbt){ - throw ErrorException(__FUNCT__,exprintf("head-Th.triangles<0 || head-Th.triangles>=Th.nbt")); - } - head = ThNew.triangles + Th.Number(head) ; - if (edge-Th.edges<0 || edge-Th.edges>=Th.nbe);{ - throw ErrorException(__FUNCT__,exprintf("edge-Th.edges<0 || edge-Th.edges>=Th.nbe")); - } - edge = ThNew.edges+ Th.Number(edge); - } - inline void GeometricalSubDomain::Set(const GeometricalSubDomain & rec,const Geometry & Gh ,const Geometry & GhNew) - { - *this = rec; - edge = Gh.Number(edge) + GhNew.edges; - } - inline void VertexOnEdge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) - { - *this = Th.VertexOnBThEdge[i]; - v = ThNew.vertices + Th.Number(v); - } - - inline void VertexOnGeom::Set(const VertexOnGeom & rec,const Triangles & Th ,Triangles & ThNew) - { - *this = rec; - mv = ThNew.vertices + Th.Number(mv); - if (gv) - if (abscisse < 0 ) - gv = ThNew.Gh.vertices + Th.Gh.Number(gv); - else - ge = ThNew.Gh.edges + Th.Gh.Number(ge); - - } - inline Real8 Edge::MetricLength() const - { - return LengthInterpole(v[0]->m,v[1]->m,v[1]->r - v[0]->r) ; - } - - inline void Triangles::ReMakeTriangleContainingTheVertex() - { - register Int4 i; - for ( i=0;iTriaAdjTriangles[a]=ta.t; - t->TriaAdjSharedEdge[a]=ta.a|l; - } - if(ta.t) { - ta.t->TriaAdjTriangles[ta.a] = t ; - ta.t->TriaAdjSharedEdge[ta.a] = a| l ; - } - } - - - inline int TriangleAdjacent::Locked() const - { return t->TriaAdjSharedEdge[a] &4;} - inline int TriangleAdjacent::Cracked() const - { return t->TriaAdjSharedEdge[a] &32;} - inline int TriangleAdjacent::GetAllFlag_UnSwap() const - { return t->TriaAdjSharedEdge[a] & 1012;} // take all flag except MarkUnSwap - - inline int TriangleAdjacent::MarkUnSwap() const - { return t->TriaAdjSharedEdge[a] &8;} - - inline void TriangleAdjacent::SetLock(){ t->SetLocked(a);} - - inline void TriangleAdjacent::SetCracked() { t->SetCracked(a);} - - inline TriangleAdjacent TriangleAdjacent::Adj() const - { return t->Adj(a);} - - inline Vertex * TriangleAdjacent::EdgeVertex(const int & i) const - {return t->TriaVertices[VerticesOfTriangularEdge[a][i]]; } - inline Vertex * TriangleAdjacent::OppositeVertex() const - {return t->TriaVertices[bamg::OppositeVertex[a]]; } - inline Icoor2 & TriangleAdjacent::det() const - { return t->det;} - inline TriangleAdjacent Adj(const TriangleAdjacent & a) - { return a.Adj();} - - inline TriangleAdjacent Next(const TriangleAdjacent & ta) - { return TriangleAdjacent(ta.t,NextEdge[ta.a]);} - - inline TriangleAdjacent Previous(const TriangleAdjacent & ta) - { return TriangleAdjacent(ta.t,PreviousEdge[ta.a]);} - - inline void Adj(GeometricalEdge * & on,int &i) - {int j=i;i=on->DirAdj[i];on=on->Adj[j];} - - inline Real4 qualite(const Vertex &va,const Vertex &vb,const Vertex &vc) - { - Real4 ret; - I2 ia=va,ib=vb,ic=vc; - I2 ab=ib-ia,bc=ic-ib,ac=ic-ia; - Icoor2 deta=Det(ab,ac); - if (deta <=0) ret = -1; - else { - Real8 a = sqrt((Real8) (ac,ac)), - b = sqrt((Real8) (bc,bc)), - c = sqrt((Real8) (ab,ab)), - p = a+b+c; - Real8 h= Max(Max(a,b),c),ro=deta/p; - ret = ro/h;} - return ret; - } - - inline Triangle::Triangle(Triangles *Th,Int4 i,Int4 j,Int4 k) { Vertex *v=Th->vertices; @@ -1119,5 +1069,4 @@ det=0; } - inline Triangle::Triangle(Vertex *v0,Vertex *v1,Vertex *v2){ TriaVertices[0]=v0; @@ -1131,5 +1080,4 @@ link=NULL;}; } - inline Real4 Triangle::qualite() { @@ -1137,81 +1085,130 @@ } - - Icoor2 inline det(const Vertex & a,const Vertex & b,const Vertex & c){ - register Icoor2 bax = b.i.x - a.i.x ,bay = b.i.y - a.i.y; - register Icoor2 cax = c.i.x - a.i.x ,cay = c.i.y - a.i.y; - return bax*cay - bay*cax; - } - - - void swap(Triangle *t1,Int1 a1, - Triangle *t2,Int1 a2, - Vertex *s1,Vertex *s2,Icoor2 det1,Icoor2 det2); - - int inline TriangleAdjacent::swap() - { return t->swap(a);} - - int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb , - TriangleAdjacent & tt1,Icoor2 & dets1, - Icoor2 & detsa,Icoor2 & detsb, int & nbswap); - - int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) ; - - inline TriangleAdjacent FindTriangleAdjacent(Edge &E){ - Vertex * a = E.v[0]; - Vertex * b = E.v[1]; - - Triangle * t = a->t; - int i = a->vint; - TriangleAdjacent ta(t,EdgesVertexTriangle[i][0]); // Previous edge - if (!t || i<0 || i>=3){ - throw ErrorException(__FUNCT__,exprintf("!t || i<0 !! i>=3")); - } - if ( a!=(*t)(i)){ - throw ErrorException(__FUNCT__,exprintf("a!=(*t)(i)")); - } - int k=0; - do { // turn around vertex in direct sens (trigo) - k++; - if (k>=20000){ - throw ErrorException(__FUNCT__,exprintf("k>=20000")); - } - // in no crack => ta.EdgeVertex(1) == a otherwise ??? - if (ta.EdgeVertex(1) == a && ta.EdgeVertex(0) == b) return ta; // find - ta = ta.Adj(); - if (ta.EdgeVertex(0) == a && ta.EdgeVertex(1) == b) return ta; // find - --ta; - } while (t != (Triangle *)ta); - throw ErrorException(__FUNCT__,exprintf("FindTriangleAdjacent: triangle not found")); - return TriangleAdjacent(0,0);//for compiler - } - - inline Vertex* TheVertex(Vertex * a){ - // give a unique vertex with smallest number - // in case on crack in mesh - Vertex * r(a), *rr; - Triangle * t = a->t; - int i = a->vint; - TriangleAdjacent ta(t,EdgesVertexTriangle[i][0]); // Previous edge - if (!t || i<0 || i>=3){ - throw ErrorException(__FUNCT__,exprintf("!t || i<0 !! i>=3")); - } - if ( a!=(*t)(i)){ - throw ErrorException(__FUNCT__,exprintf("a!=(*t)(i)")); - } - int k=0; - do { // turn around vertex in direct sens (trigo) - k++; - if (k>=20000){ - throw ErrorException(__FUNCT__,exprintf("k>=20000")); - } - // in no crack => ta.EdgeVertex(1) == a - if ((rr=ta.EdgeVertex(0)) < r) r = rr; - ta = ta.Adj(); - if ((rr=ta.EdgeVertex(1)) < r) r =rr; - --ta; - } while (t != (Triangle*) ta); - return r; - } + /*}}}1*/ + /*INLINE functions of CLASS Edge{{{1*/ + inline void Edge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) + { + *this = Th.edges[i]; + v[0] = ThNew.vertices + Th.Number(v[0]); + v[1] = ThNew.vertices + Th.Number(v[1]); + if (onGeometry) + onGeometry = ThNew.Gh.edges+Th.Gh.Number(onGeometry); + if (adj[0]) adj[0] = ThNew.edges + Th.Number(adj[0]); + if (adj[1]) adj[1] = ThNew.edges + Th.Number(adj[1]); + + } + inline Real8 Edge::MetricLength() const + { + return LengthInterpole(v[0]->m,v[1]->m,v[1]->r - v[0]->r) ; + } + + /*}}}1*/ + /*INLINE functions of CLASS GeometricalEdge{{{1*/ + inline void GeometricalEdge::Set(const GeometricalEdge & rec,const Geometry & Gh ,Geometry & GhNew) + { + *this = rec; + v[0] = GhNew.vertices + Gh.Number(v[0]); + v[1] = GhNew.vertices + Gh.Number(v[1]); + if (Adj[0]) Adj[0] = GhNew.edges + Gh.Number(Adj[0]); + if (Adj[1]) Adj[1] = GhNew.edges + Gh.Number(Adj[1]); + } + /*}}}1*/ + /*INLINE functions of CLASS Curve{{{1*/ + inline void Curve::Set(const Curve & rec,const Geometry & Gh ,Geometry & GhNew) + { + *this = rec; + be = GhNew.edges + Gh.Number(be); + ee = GhNew.edges + Gh.Number(ee); + if(next) next= GhNew.curves + Gh.Number(next); + } + /*}}}1*/ + /*INLINE functions of CLASS SubDomain{{{1*/ + inline void SubDomain::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) + { + *this = Th.subdomains[i]; + if ( head-Th.triangles<0 || head-Th.triangles>=Th.nbt){ + throw ErrorException(__FUNCT__,exprintf("head-Th.triangles<0 || head-Th.triangles>=Th.nbt")); + } + head = ThNew.triangles + Th.Number(head) ; + if (edge-Th.edges<0 || edge-Th.edges>=Th.nbe);{ + throw ErrorException(__FUNCT__,exprintf("edge-Th.edges<0 || edge-Th.edges>=Th.nbe")); + } + edge = ThNew.edges+ Th.Number(edge); + } + /*}}}1*/ + /*INLINE functions of CLASS GeometricalSubDomain{{{1*/ + inline void GeometricalSubDomain::Set(const GeometricalSubDomain & rec,const Geometry & Gh ,const Geometry & GhNew) + { + *this = rec; + edge = Gh.Number(edge) + GhNew.edges; + } + /*}}}1*/ + /*INLINE functions of CLASS VertexOnEdge{{{1*/ + inline void VertexOnEdge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) + { + *this = Th.VertexOnBThEdge[i]; + v = ThNew.vertices + Th.Number(v); + } + /*}}}1*/ + /*INLINE functions of CLASS VertexOnGeom{{{1*/ + inline void VertexOnGeom::Set(const VertexOnGeom & rec,const Triangles & Th ,Triangles & ThNew) + { + *this = rec; + mv = ThNew.vertices + Th.Number(mv); + if (gv) + if (abscisse < 0 ) + gv = ThNew.Gh.vertices + Th.Gh.Number(gv); + else + ge = ThNew.Gh.edges + Th.Gh.Number(ge); + + } + /*}}}1*/ + /*INLINE functions of CLASS TriangleAdjacent{{{1*/ + inline void TriangleAdjacent::SetAdj2(const TriangleAdjacent & ta, int l ){ + //set Adjacent Triangle of a triangle + if(t) { + t->TriaAdjTriangles[a]=ta.t; + t->TriaAdjSharedEdge[a]=ta.a|l; + } + if(ta.t) { + ta.t->TriaAdjTriangles[ta.a] = t ; + ta.t->TriaAdjSharedEdge[ta.a] = a| l ; + } + } + inline int TriangleAdjacent::Locked() const + { return t->TriaAdjSharedEdge[a] &4;} + inline int TriangleAdjacent::Cracked() const + { return t->TriaAdjSharedEdge[a] &32;} + inline int TriangleAdjacent::GetAllFlag_UnSwap() const + { return t->TriaAdjSharedEdge[a] & 1012;} // take all flag except MarkUnSwap + inline int TriangleAdjacent::MarkUnSwap() const + { return t->TriaAdjSharedEdge[a] &8;} + inline void TriangleAdjacent::SetLock(){ t->SetLocked(a);} + inline void TriangleAdjacent::SetCracked() { t->SetCracked(a);} + inline TriangleAdjacent TriangleAdjacent::Adj() const + { return t->Adj(a);} + inline Vertex * TriangleAdjacent::EdgeVertex(const int & i) const + {return t->TriaVertices[VerticesOfTriangularEdge[a][i]]; } + inline Vertex * TriangleAdjacent::OppositeVertex() const + {return t->TriaVertices[bamg::OppositeVertex[a]]; } + inline Icoor2 & TriangleAdjacent::det() const + { return t->det;} + int inline TriangleAdjacent::swap() + { return t->swap(a);} + /*}}}1*/ + + /*Other prototypes {{{1*/ + TriangleAdjacent CloseBoundaryEdge(I2 ,Triangle *, double &,double &) ; + TriangleAdjacent CloseBoundaryEdgeV2(I2 A,Triangle *t, double &a,double &b); + Int4 FindTriangle(Triangles &Th, Real8 x, Real8 y, double* a,int & inside); + void swap(Triangle *t1,Int1 a1, + Triangle *t2,Int1 a2, + Vertex *s1,Vertex *s2,Icoor2 det1,Icoor2 det2); + int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb , + TriangleAdjacent & tt1,Icoor2 & dets1, + Icoor2 & detsa,Icoor2 & detsb, int & nbswap); + int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) ; + /*}}}1*/ + } #endif