diff --git a/contrib/bamg/bamglib/Mesh2.h b/contrib/bamg/bamglib/Mesh2.h
index ce475c00063aa043b4bbc474c11fe77fe3309b34..3c4afa346d185cb181da6631af1d2ce128310c3d 100644
--- a/contrib/bamg/bamglib/Mesh2.h
+++ b/contrib/bamg/bamglib/Mesh2.h
@@ -1,26 +1,26 @@
// -*- Mode : c++ -*-
//
-// SUMMARY :
-// USAGE :
-// ORG :
+// SUMMARY :
+// USAGE :
+// ORG :
// AUTHOR : Frederic Hecht
// E-MAIL : hecht@ann.jussieu.fr
//
/*
-
+
This file is part of Freefem++
-
+
Freefem++ is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.
-
+
Freefem++ is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
-
+
You should have received a copy of the GNU Lesser General Public License
along with Freefem++; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
@@ -74,7 +74,7 @@ inline int BinaryRand(){
return rand() & 16384; // 2^14 (for sun because RAND_MAX is not def in stdlib.h)
#endif
-}
+}
typedef P2<Real8,Real8> R2;
typedef P2xP2<Int2,Int4> I2xI2;
typedef P2<Real4,Real8> R2xR2;
@@ -88,7 +88,7 @@ inline float OppositeAngle(float a)
{return a<0 ? fPi + a :a - fPi ;}
inline double OppositeAngle(double a)
{return a<0 ? Pi + a :a - Pi ;}
-
+
#ifdef DRAWING
extern Real4 xGrafCoef,yGrafCoef,xGrafOffSet,yGrafOffSet;
extern R2 GrafPMin,GrafPMax;
@@ -97,13 +97,13 @@ extern Real8 Grafh;
Icoor2 inline det(const I2 &a,const I2 & b,const I2 &c)
{
- register Icoor2 bax = b.x - a.x ,bay = b.y - a.y;
- register Icoor2 cax = c.x - a.x ,cay = c.y - a.y;
+ register Icoor2 bax = b.x - a.x ,bay = b.y - a.y;
+ register Icoor2 cax = c.x - a.x ,cay = c.y - a.y;
return bax*cay - bay*cax;}
-// def de numerotation dans un triangles
+// def de numerotation dans un triangles
static const Int2 VerticesOfTriangularEdge[3][2] = {{1,2},{2,0},{0,1}};
static const Int2 EdgesVertexTriangle[3][2] = {{1,2},{2,0},{0,1}};
static const Int2 OppositeVertex[3] = {0,1,2};
@@ -133,7 +133,7 @@ const int IsVertexOnEdge = 32;
/////////////////////////////////////////////////////////////////////////////////////
#ifndef NOTFREEFEM
class ErrorMesh : public Error
-{
+{
public:
Triangles *Th;
ErrorMesh(const char * Text,int l,Triangles * TTh=0, const char *t2="") :
@@ -141,17 +141,17 @@ public:
};
#endif
-class Direction { //
+class Direction { //
private:
Icoor1 dir;
public:
Direction(): dir(MaxICoor){}; // no direction set
- Direction(Icoor1 i,Icoor1 j) { Icoor2 n2 = 2*(Abs(i)+Abs(j));
+ Direction(Icoor1 i,Icoor1 j) { Icoor2 n2 = 2*(Abs(i)+Abs(j));
Icoor2 r = MaxICoor* (Icoor2) i;
- Icoor1 r1 = (Icoor1) (2*(r/ n2)); // odd number
+ Icoor1 r1 = (Icoor1) (2*(r/ n2)); // odd number
dir = (j>0) ? r1 : r1+1; // odd -> j>0 even -> j<0
}
- int sens( Icoor1 i,Icoor1 j) { int r =1;
+ int sens( Icoor1 i,Icoor1 j) { int r =1;
if (dir!= MaxICoor) {
Icoor2 x(dir/2),y1(MaxICoor/2-Abs(x)),y(dir%2?-y1:y1);
r = (x*i + y*j) >=0;}
@@ -166,11 +166,11 @@ class Direction { //
rmoveto(D.x,D.y);
}
}
-#endif
-
-
-
-
+#endif
+
+
+
+
};
/////////////////////////////////////////////////////////////////////////////////////
class Vertex {public:
@@ -182,18 +182,18 @@ class Vertex {public:
union {
Triangle * t; // one triangle which contained the vertex
Int4 color;
- Vertex * to;// use in geometry Vertex to now the Mesh Vertex associed
+ Vertex * to;// use in geometry Vertex to now the Mesh Vertex associed
VertexOnGeom * on; // if vint 8; // set with Triangles::SetVertexFieldOn()
Vertex * onbv; // if vint == 16 on Background vertex Triangles::SetVertexFieldOnBTh()
VertexOnEdge * onbe; // if vint == 32 on Background edge
};
Int1 vint; // the vertex number in triangle; varies between 0 and 2 in t
- operator I2 () const {return i;} // operator de cast
- operator const R2 & () const {return r;}// operator de cast
-// operator R2 & () {return r;}// operator de cast
+ operator I2 () const {return i;} // operator de cast
+ operator const R2 & () const {return r;}// operator de cast
+// operator R2 & () {return r;}// operator de cast
Real8 operator()(R2 x) const { return m(x);}
- operator Metric () const {return m;}// operator de cast
- Int4 Optim(int = 1,int =0);
+ operator Metric () const {return m;}// operator de cast
+ Int4 Optim(int = 1,int =0);
// Vertex(){}
// ~Vertex(){}
Real8 Smoothing(Triangles & ,const Triangles & ,Triangle * & ,Real8 =1);
@@ -202,12 +202,12 @@ class Vertex {public:
friend ostream& operator <<(ostream& f, const Vertex & v)
{f << "(" << v.i << "," << v.r << MatVVP2x2(v.m) << ")" ; return f;}
inline void Set(const Vertex & rec,const Triangles &,Triangles &);
-
+
#ifdef DRAWING
void Draw(Int4 =-1) const ;
void MoveTo() const { rmoveto(r.x,r.y); }
void LineTo() const { rlineto(r.x,r.y); }
-#endif
+#endif
};
double QuadQuality(const Vertex &,const Vertex &,const Vertex &,const Vertex &);
@@ -221,21 +221,21 @@ class TriangleAdjacent {
return f;}
public:
- Triangle * t; // le triangle
+ Triangle * t; // le triangle
int a; // le numero de l arete
-
+
TriangleAdjacent(Triangle * tt,int aa): t(tt),a(aa &3) {};
TriangleAdjacent() {};
-
+
operator Triangle * () const {return t;}
operator Triangle & () const {return *t;}
operator int() const {return a;}
- TriangleAdjacent & operator++()
+ TriangleAdjacent & operator++()
{
a= NextEdge[a];
return *this;}
TriangleAdjacent operator--()
- {
+ {
a= PreviousEdge[a];
return *this;}
inline TriangleAdjacent Adj() const ;
@@ -251,7 +251,7 @@ public:
inline void SetMarkUnSwap();
inline void SetCracked();
inline int Cracked() const ;
-};// end of Vertex class
+};// end of Vertex class
/////////////////////////////////////////////////////////////////////////////////////
@@ -262,35 +262,35 @@ class Edge { public:
Vertex & operator[](int i){return *v[i];};
Vertex * operator()(int i){return v[i];};
- void ReNumbering(Vertex *vb,Vertex *ve, Int4 *renu)
+ void ReNumbering(Vertex *vb,Vertex *ve, Int4 *renu)
{
if (v[0] >=vb && v[0] <ve) v[0] = vb + renu[v[0]-vb];
if (v[1] >=vb && v[1] <ve) v[1] = vb + renu[v[1]-vb];
}
const Vertex & operator[](int i) const { return *v[i];};
- R2 operator()(double t) const; // return the point
+ R2 operator()(double t) const; // return the point
// on the curve edge a t in [0:1]
- Edge * adj[2]; // the 2 adj edges if on the same curve
- int Intersection(const Edge & e) const {
- if (!(adj[0]==&e || adj[1]==&e))
- cerr << "Bug : Intersection " << (void*) &e << " "
+ Edge * adj[2]; // the 2 adj edges if on the same curve
+ int Intersection(const Edge & e) const {
+ if (!(adj[0]==&e || adj[1]==&e))
+ cerr << "Bug : Intersection " << (void*) &e << " "
<< adj[0] << " " << adj[1] << endl;
assert(adj[0]==&e || adj[1]==&e);
return adj[0]==&e ? 0 : 1;}
- Real8 MetricLength() const ;
+ Real8 MetricLength() const ;
inline void Set(const Triangles &,Int4,Triangles &);
-
+
#ifdef DRAWING
void Draw(Int4 = -1) const ;
#endif
-}; // end of Edge class
+}; // end of Edge class
/////////////////////////////////////////////////////////////////////////////////////
class GeometricalVertex :public Vertex {
int cas;
friend class Geometry;
- GeometricalVertex * link; // link all the same GeometricalVertex circular (Crack)
+ GeometricalVertex * link; // link all the same GeometricalVertex circular (Crack)
public:
int Corner() const {return cas&4;}
int Required()const {return cas&6;}// a corner is required
@@ -299,10 +299,10 @@ public:
void Set(){cas=0;}
GeometricalVertex() :cas(0), link(this) {};
GeometricalVertex * The() { assert(link); return link;}// return a unique vertices
- int IsThe() const { return link == this;}
-
+ int IsThe() const { return link == this;}
+
inline void Set(const GeometricalVertex & rec,const Geometry & Gh ,const Geometry & GhNew);
- inline friend ostream& operator <<(ostream& f, const GeometricalVertex & s)
+ inline friend ostream& operator <<(ostream& f, const GeometricalVertex & s)
{ f << s.r << "," << s.cas << ".";return f; }
};
@@ -312,20 +312,20 @@ class GeometricalEdge {
GeometricalVertex * v[2];
Int4 ref;
Int4 CurveNumber;
- R2 tg[2]; // the 2 tangente
- // if tg[0] =0 => no continuite
- GeometricalEdge * Adj [2];
+ R2 tg[2]; // the 2 tangente
+ // if tg[0] =0 => no continuite
+ GeometricalEdge * Adj [2];
int SensAdj[2];
// private:
int flag ;
- public:
+ public:
GeometricalEdge * link; // if Cracked() or Equi()
-// end of data
-
+// end of data
+
GeometricalVertex & operator[](int i){return *v[i];};
const GeometricalVertex & operator[](int i) const { return *v[i];};
- GeometricalVertex * operator()(int i){return v[i];};
+ GeometricalVertex * operator()(int i){return v[i];};
// inline void Set(const Geometry &,Int4,Geometry &);
R2 F(Real8 theta) const ; // parametrization of the curve edge
@@ -340,7 +340,7 @@ class GeometricalEdge {
int Mark()const {return flag &16;}
int Required() { return flag & 64;}
void SetCracked() { flag |= 1;}
- void SetDup() { flag |= 32;} // not a real edge
+ void SetDup() { flag |= 32;} // not a real edge
void SetEqui() { flag |= 2;}
void SetTgA() { flag|=4;}
void SetTgB() { flag|=8;}
@@ -348,26 +348,26 @@ class GeometricalEdge {
void SetUnMark() { flag &= 1007 /* 1023-16*/;}
void SetRequired() { flag |= 64;}
void SetReverseEqui() {flag |= 128;}
-
+
inline void Set(const GeometricalEdge & rec,const Geometry & Th ,Geometry & ThNew);
-#ifdef DRAWING
+#ifdef DRAWING
void Draw(Int4 =-1);
#endif
-
+
};
-
+
class Curve {public:
GeometricalEdge * be,*ee; // begin et end edge
int kb,ke; // begin vetex and end vertex
Curve *next; // next curve equi to this
- bool master; // true => of equi curve point on this curve
+ bool master; // true => of equi curve point on this curve
inline void Set(const Curve & rec,const Geometry & Th ,Geometry & ThNew);
- Curve() : be(0),ee(0),kb(0),ke(0),next(0),master(true) {}
- void Reverse() { Exchange(be,ee); Exchange(kb,ke);} // revese the sens of the curse
+ Curve() : be(0),ee(0),kb(0),ke(0),next(0),master(true) {}
+ void Reverse() { Exchange(be,ee); Exchange(kb,ke);} // revese the sens of the curse
};
-
-
+
+
/////////////////////////////////////////////////////////////////////////////////////
class Triangle {
@@ -377,11 +377,11 @@ class Triangle {
private: // les arete sont opposes a un sommet
Vertex * ns [3]; // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer
- Triangle * at [3]; // nu triangle adjacent
+ Triangle * at [3]; // nu triangle adjacent
Int1 aa[3]; // les nu des arete dans le triangles (mod 4)
- public:
+ public:
Icoor2 det; // determinant du triangle (2 fois l aire des vertices entieres)
- union {
+ union {
Triangle * link ;
Int4 color;
};
@@ -395,67 +395,67 @@ class Triangle {
inline int In(Vertex *v) const { return ns[0]==v || ns[1]==v || ns[2]==v ;}
TriangleAdjacent FindBoundaryEdge(int ) const;
- void ReNumbering(Triangle *tb,Triangle *te, Int4 *renu)
+ void ReNumbering(Triangle *tb,Triangle *te, Int4 *renu)
{
if (link >=tb && link <te) link = tb + renu[link -tb];
if (at[0] >=tb && at[0] <te) at[0] = tb + renu[at[0]-tb];
if (at[1] >=tb && at[1] <te) at[1] = tb + renu[at[1]-tb];
- if (at[2] >=tb && at[2] <te) at[2] = tb + renu[at[2]-tb];
+ if (at[2] >=tb && at[2] <te) at[2] = tb + renu[at[2]-tb];
}
- void ReNumbering(Vertex *vb,Vertex *ve, Int4 *renu)
+ void ReNumbering(Vertex *vb,Vertex *ve, Int4 *renu)
{
if (ns[0] >=vb && ns[0] <ve) ns[0] = vb + renu[ns[0]-vb];
if (ns[1] >=vb && ns[1] <ve) ns[1] = vb + renu[ns[1]-vb];
- if (ns[2] >=vb && ns[2] <ve) ns[2] = vb + renu[ns[2]-vb];
+ if (ns[2] >=vb && ns[2] <ve) ns[2] = vb + renu[ns[2]-vb];
}
const Vertex & operator[](int i) const {return *ns[i];};
Vertex & operator[](int i) {return *ns[i];};
-
+
const Vertex * operator()(int i) const {return ns[i];};
Vertex * & operator()(int i) {return ns[i];};
-
- TriangleAdjacent Adj(int i) const // triangle adjacent + arete
+
+ TriangleAdjacent Adj(int i) const // triangle adjacent + arete
{ return TriangleAdjacent(at[i],aa[i]&3);};
- Triangle * TriangleAdj(int i) const
- {return at[i&3];} // triangle adjacent + arete
- Int1 NuEdgeTriangleAdj(int i) const
- {return aa[i&3]&3;} // Number of the adjacent edge in adj tria
+ Triangle * TriangleAdj(int i) const
+ {return at[i&3];} // triangle adjacent + arete
+ Int1 NuEdgeTriangleAdj(int i) const
+ {return aa[i&3]&3;} // Number of the adjacent edge in adj tria
inline Real4 qualite() ;
-
- void SetAdjAdj(Int1 a)
+
+ void SetAdjAdj(Int1 a)
{ a &= 3;
register Triangle *tt=at[a];
aa [a] &= 55; // remove MarkUnSwap
register Int1 aatt = aa[a] & 3;
- if(tt){
+ if(tt){
tt->at[aatt]=this;
- tt->aa[aatt]=a + (aa[a] & 60 ) ;}// Copy all the mark
+ tt->aa[aatt]=a + (aa[a] & 60 ) ;}// Copy all the mark
}
-
+
void SetAdj2(Int1 a,Triangle *t,Int1 aat)
{ at[a]=t;aa[a]=aat;
if(t) {t->at[aat]=this;t->aa[aat]=a;}
}
-
+
void SetTriangleContainingTheVertex()
- {
+ {
if (ns[0]) (ns[0]->t=this,ns[0]->vint=0);
if (ns[1]) (ns[1]->t=this,ns[1]->vint=1);
if (ns[2]) (ns[2]->t=this,ns[2]->vint=2);
}
-
+
int swap(Int2 a1,int=0);
Int4 Optim(Int2 a,int =0);
- int Locked(int a)const { return aa[a]&4;}
- int Hidden(int a)const { return aa[a]&16;}
+ int Locked(int a)const { return aa[a]&4;}
+ int Hidden(int a)const { return aa[a]&16;}
int Cracked(int a) const { return aa[a] & 32;}
- // for optimisation
+ // for optimisation
int GetAllflag(int a){return aa[a] & 1020;}
void SetAllFlag(int a,int f){aa[a] = (aa[a] &3) + (1020 & f);}
@@ -467,7 +467,7 @@ class Triangle {
register Triangle * t = at[a];
if(t) t->aa[aa[a] & 3] |=32;
aa[a] |= 32;}
-
+
double QualityQuad(int a,int option=1) const;
Triangle * Quadrangle(Vertex * & v0,Vertex * & v1,Vertex * & v2,Vertex * & v3) const ;
@@ -482,15 +482,15 @@ class Triangle {
aa[a] |=8 ;}
- void SetUnMarkUnSwap(int a){
+ void SetUnMarkUnSwap(int a){
register Triangle * t = at[a];
- t->aa[aa[a] & 3] &=55; // 23 + 32
+ t->aa[aa[a] & 3] &=55; // 23 + 32
aa[a] &=55 ;}
-
-
-#ifdef DEBUG
- void inline checka(Int1 a);
+
+
+#ifdef DEBUG
+ void inline checka(Int1 a);
void inline check();
#endif
@@ -500,7 +500,7 @@ class Triangle {
#endif
-}; // end of Triangle class
+}; // end of Triangle class
@@ -508,7 +508,7 @@ class Triangle {
class ListofIntersectionTriangles {
/////////////////////////////////////////////////////////////////////////////////////
class IntersectionTriangles {
-public:
+public:
Triangle *t;
Real8 bary[3]; // use if t != 0
R2 x;
@@ -524,13 +524,13 @@ class SegInterpolation {
Real8 sBegin,sEnd; // abscisse of the seg on edge parameter
Real8 lBegin,lEnd; // length abscisse set in ListofIntersectionTriangles::Length
int last;// last index in ListofIntersectionTriangles for this Sub seg of edge
- R2 F(Real8 s){
+ R2 F(Real8 s){
Real8 c01=lEnd-lBegin, c0=(lEnd-s)/c01, c1=(s-lBegin)/c01;
assert(lBegin<= s && s <=lEnd);
return e->F(sBegin*c0+sEnd*c1);}
};
-
- int MaxSize; //
+
+ int MaxSize; //
int Size; //
Real8 len; //
int state;
@@ -543,60 +543,60 @@ class SegInterpolation {
operator int&() {return Size;}
ListofIntersectionTriangles(int n=256,int m=16)
: MaxSize(n), Size(0), len(-1),state(-1),lIntTria(new IntersectionTriangles[n]) ,
- NbSeg(0), MaxNbSeg(m), lSegsI(new SegInterpolation[m])
- { if (verbosity>9)
+ NbSeg(0), MaxNbSeg(m), lSegsI(new SegInterpolation[m])
+ { if (verbosity>9)
cout << " construct ListofIntersectionTriangles"
<< MaxSize << " " << MaxNbSeg<< endl;};
~ListofIntersectionTriangles(){
if (lIntTria) delete [] lIntTria,lIntTria=0;
- if (lSegsI) delete [] lSegsI,lSegsI=0;}
+ if (lSegsI) delete [] lSegsI,lSegsI=0;}
void init(){state=0;len=0;Size=0;}
-
+
int NewItem(Triangle * tt,Real8 d0,Real8 d1,Real8 d2);
int NewItem(R2,const Metric & );
- void NewSubSeg(GeometricalEdge *e,Real8 s0,Real8 s1)
- {
+ void NewSubSeg(GeometricalEdge *e,Real8 s0,Real8 s1)
+ {
if (NbSeg>=MaxNbSeg) {
int mneo= MaxNbSeg;
MaxNbSeg *= 2;
- if (verbosity>3)
- cout <<" reshape lSegsI from " << mneo << " to "
+ if (verbosity>3)
+ cout <<" reshape lSegsI from " << mneo << " to "
<< MaxNbSeg <<endl;
SegInterpolation * lEn = new SegInterpolation[MaxNbSeg];
assert(lSegsI && NbSeg < MaxNbSeg);
- for (int i=0;i< NbSeg;i++)
- lEn[i] = lSegsI[MaxNbSeg]; // copy old to new
+ for (int i=0;i< NbSeg;i++)
+ lEn[i] = lSegsI[MaxNbSeg]; // copy old to new
delete [] lSegsI; // remove old
- lSegsI = lEn;
+ lSegsI = lEn;
}
- if (NbSeg)
+ if (NbSeg)
lSegsI[NbSeg-1].last=Size;
lSegsI[NbSeg].e=e;
lSegsI[NbSeg].sBegin=s0;
- lSegsI[NbSeg].sEnd=s1;
- NbSeg++;
+ lSegsI[NbSeg].sEnd=s1;
+ NbSeg++;
}
-
+
// void CopyMetric(int i,int j){ lIntTria[j].m=lIntTria[i].m;}
// void CopyMetric(const Metric & mm,int j){ lIntTria[j].m=mm;}
- void ReShape() {
+ void ReShape() {
register int newsize = MaxSize*2;
IntersectionTriangles * nw = new IntersectionTriangles[newsize];
assert(nw);
for (int i=0;i<MaxSize;i++) // recopy
- nw[i] = lIntTria[i];
+ nw[i] = lIntTria[i];
if(verbosity>3)
- cout << " ListofIntersectionTriangles ReShape MaxSize "
- << MaxSize << " -> "
+ cout << " ListofIntersectionTriangles ReShape MaxSize "
+ << MaxSize << " -> "
<< newsize << endl;
- MaxSize = newsize;
+ MaxSize = newsize;
delete [] lIntTria;// remove old
lIntTria = nw; // copy pointer
}
-
- void SplitEdge(const Triangles & ,const R2 &,const R2 &,int nbegin=0);
- Real8 Length();
+
+ void SplitEdge(const Triangles & ,const R2 &,const R2 &,int nbegin=0);
+ Real8 Length();
Int4 NewPoints(Vertex *,Int4 & nbv,Int4 nbvx);
};
@@ -608,35 +608,35 @@ public:
int sens; // -1 or 1
Int4 ref;
inline void Set(const GeometricalSubDomain &,const Geometry & ,const Geometry &);
-
+
};
/////////////////////////////////////////////////////////////////////////////////////
class SubDomain {
public:
Triangle * head;
- Int4 ref;
+ Int4 ref;
int sens; // -1 or 1
- Edge * edge; // to geometrical
+ Edge * edge; // to geometrical
inline void Set(const Triangles &,Int4,Triangles &);
-
+
};
/////////////////////////////////////////////////////////////////////////////////////
class VertexOnGeom { public:
Vertex * mv;
- Real8 abscisse;
- union{
- GeometricalVertex * gv; // if abscisse <0;
+ Real8 abscisse;
+ union{
+ GeometricalVertex * gv; // if abscisse <0;
GeometricalEdge * ge; // if abscisse in [0..1]
};
- inline void Set(const VertexOnGeom&,const Triangles &,Triangles &);
+ inline void Set(const VertexOnGeom&,const Triangles &,Triangles &);
int OnGeomVertex()const {return this? abscisse <0 :0;}
int OnGeomEdge() const {return this? abscisse >=0 :0;}
- VertexOnGeom(): mv(0),abscisse(0){gv=0;}
+ VertexOnGeom(): mv(0),abscisse(0){gv=0;}
VertexOnGeom(Vertex & m,GeometricalVertex &g) : mv(&m),abscisse(-1){gv=&g;}
- // cout << " mv = " <<mv << " gv = " << gv << endl;}
+ // cout << " mv = " <<mv << " gv = " << gv << endl;}
VertexOnGeom(Vertex & m,GeometricalEdge &g,Real8 s) : mv(&m),abscisse(s){ge=&g;}
- //cout << &g << " " << ge << endl;}
+ //cout << &g << " " << ge << endl;}
operator Vertex * () const {return mv;}
operator GeometricalVertex * () const {return gv;}
operator GeometricalEdge * () const {return ge;}
@@ -650,7 +650,7 @@ class VertexOnGeom { public:
// else f << *vog.ge << " ;; " ;
return f;}
inline void Set(const Triangles &,Int4,Triangles &);
-
+
};
/////////////////////////////////////////////////////////////////////////////////////
class VertexOnVertex {public:
@@ -667,59 +667,59 @@ class VertexOnEdge {public:
Real8 abcisse;
VertexOnEdge( Vertex * w, Edge *bw,Real8 s) :v(w),be(bw),abcisse(s) {}
VertexOnEdge(){}
- inline void Set(const Triangles &,Int4,Triangles &);
- void SetOnBTh(){v->onbe=this;v->vint=IsVertexOnEdge;}
+ inline void Set(const Triangles &,Int4,Triangles &);
+ void SetOnBTh(){v->onbe=this;v->vint=IsVertexOnEdge;}
Vertex & operator[](int i) const { return (*be)[i];}
operator Real8 () const { return abcisse;}
- operator Vertex * () const { return v;}
- operator Edge * () const { return be;}
+ operator Vertex * () const { return v;}
+ operator Edge * () const { return be;}
};
inline TriangleAdjacent FindTriangleAdjacent(Edge &E);
- inline Vertex * TheVertex(Vertex * a); // for remove crak in mesh
+ inline Vertex * TheVertex(Vertex * a); // for remove crak in mesh
/////////////////////////////////////////////////////////////////////////////////////
-
-class CrackedEdge { // a small class to store on crack an uncrack information
+
+class CrackedEdge { // a small class to store on crack an uncrack information
friend class Triangles;
- friend ostream& operator <<(ostream& f, const Triangles & Th) ;
+ friend ostream& operator <<(ostream& f, const Triangles & Th) ;
class CrackedTriangle {
friend class Triangles;
friend class CrackedEdge;
- friend ostream& operator <<(ostream& f, const Triangles & Th) ;
+ friend ostream& operator <<(ostream& f, const Triangles & Th) ;
Triangle * t; // edge of triangle t
int i; // edge number of in triangle
- Edge *edge; // the 2 edge
- Vertex *New[2]; // new vertex number
- CrackedTriangle() : t(0),i(0),edge(0) { New[0]=New[1]=0;}
- CrackedTriangle(Edge * a) : t(0),i(0),edge(a) { New[0]=New[1]=0;}
- void Crack(){
- Triangle & T(*t);
+ Edge *edge; // the 2 edge
+ Vertex *New[2]; // new vertex number
+ CrackedTriangle() : t(0),i(0),edge(0) { New[0]=New[1]=0;}
+ CrackedTriangle(Edge * a) : t(0),i(0),edge(a) { New[0]=New[1]=0;}
+ void Crack(){
+ Triangle & T(*t);
int i0=VerticesOfTriangularEdge[i][0];
int i1=VerticesOfTriangularEdge[i][0];
assert(New[0] && New[1]);
T(i0) = New[0];
- T(i1) = New[1];}
- void UnCrack(){
- Triangle & T(*t);
+ T(i1) = New[1];}
+ void UnCrack(){
+ Triangle & T(*t);
int i0=VerticesOfTriangularEdge[i][0];
int i1=VerticesOfTriangularEdge[i][0];
assert(New[0] && New[1]);
T(i0) = TheVertex(T(i0));
- T(i1) = TheVertex(T(i1));}
+ T(i1) = TheVertex(T(i1));}
void Set() {
TriangleAdjacent ta ( FindTriangleAdjacent(*edge));
t = ta;
i = ta;
-
+
New[0] = ta.EdgeVertex(0);
New[1] = ta.EdgeVertex(1);
- // warning the ref
-
- }
-
+ // warning the ref
+
+ }
+
}; // end of class CrackedTriangle
- public:
- CrackedTriangle a,b;
+ public:
+ CrackedTriangle a,b;
CrackedEdge() :a(),b() {}
CrackedEdge(Edge * start, Int4 i,Int4 j) : a(start+i),b(start+j) {};
CrackedEdge(Edge * e0, Edge * e1 ) : a(e0),b(e1) {};
@@ -731,35 +731,35 @@ class CrackedEdge { // a small class to store on crack an uncrack information
/////////////////////////////////////////////////////////////////////////////////////
-class Triangles {
+class Triangles {
public:
enum TypeFileMesh {
AutoMesh=0,BDMesh=1,NOPOMesh=2,amMesh=3,am_fmtMesh=4,amdbaMesh=5,
ftqMesh=6,mshMesh=7};
- int static counter; // to kown the number of mesh in memory
- int OnDisk; // true if on disk
+ int static counter; // to kown the number of mesh in memory
+ int OnDisk; // true if on disk
Geometry & Gh; // Geometry
- Triangles & BTh; // Background Mesh Bth==*this =>no background
-
+ Triangles & BTh; // Background Mesh Bth==*this =>no background
+
Int4 NbRef; // counter of ref on the this class if 0 we can delete
Int4 nbvx,nbtx; // nombre max de sommets , de triangles
-
- Int4 nt,nbv,nbt,nbiv,nbe; // nb of legal triangles, nb of vertex, of triangles,
+
+ Int4 nt,nbv,nbt,nbiv,nbe; // nb of legal triangles, nb of vertex, of triangles,
// of initial vertices, of edges with reference,
- Int4 NbOfQuad; // nb of quadrangle
+ Int4 NbOfQuad; // nb of quadrangle
- Int4 NbSubDomains; //
+ Int4 NbSubDomains; //
Int4 NbOutT; // Nb of oudeside triangle
Int4 NbOfTriangleSearchFind;
Int4 NbOfSwapTriangle;
char * name, *identity;
Vertex * vertices; // data of vertices des sommets
-
+
Int4 NbVerticesOnGeomVertex;
VertexOnGeom * VerticesOnGeomVertex;
-
+
Int4 NbVerticesOnGeomEdge;
VertexOnGeom * VerticesOnGeomEdge;
@@ -769,36 +769,36 @@ public:
Int4 NbVertexOnBThEdge;
VertexOnEdge *VertexOnBThEdge;
-
+
Int4 NbCrackedVertices;
-
+
Int4 NbCrackedEdges;
CrackedEdge *CrackedEdges;
-
-
+
+
R2 pmin,pmax; // extrema
Real8 coefIcoor; // coef to integer Icoor1;
Triangle * triangles;
- Edge * edges;
+ Edge * edges;
QuadTree *quadtree;
Vertex ** ordre;
SubDomain * subdomains;
ListofIntersectionTriangles lIntTria;
// end of variable
-
+
Triangles(Int4 i);//:BTh(*this),Gh(*new Geometry()){PreInit(i);}
-
- ~Triangles();
+
+ ~Triangles();
Triangles(const char * ,Real8=-1) ;
-
+
Triangles(Int4 nbvx,Triangles & BT,int keepBackVertices=1)
:Gh(BT.Gh),BTh(BT) {
try {GeomToTriangles1(nbvx,keepBackVertices);}
catch(...) { this->~Triangles(); throw; } }
-
+
Triangles(Int4 nbvx,Geometry & G)
:Gh(G),BTh(*this){
try { GeomToTriangles0(nbvx);}
@@ -806,13 +806,13 @@ public:
Triangles(Triangles &,Geometry * pGh=0,Triangles* pBTh=0,Int4 nbvxx=0 ); // COPY OPERATEUR
// Triangles(Triangles &){ cerr << " BUG call copy opretor of Triangles" << endl;MeshError(111);}
Triangles(const Triangles &,const int *flag,const int *bb); // truncature
-
+
void SetIntCoor(const char * from =0);
// void RandomInit();
// void CubeInit(int ,int);
-
+
Real8 MinimalHmin() {return 2.0/coefIcoor;}
Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);}
const Vertex & operator[] (Int4 i) const { return vertices[i];};
@@ -832,7 +832,7 @@ public:
void FindSubDomain(int );
Int4 ConsRefTriangle(Int4 *) const;
void ShowHistogram() const;
- void ShowRegulaty() const; // Add FH avril 2007
+ void ShowRegulaty() const; // Add FH avril 2007
// void ConsLinkTriangle();
void ReMakeTriangleContainingTheVertex();
@@ -842,14 +842,14 @@ public:
void MaxSubDivision(Real8 maxsubdiv);
void WriteMetric(ostream &,int iso) ;
Edge** MakeGeometricalEdgeToEdge();
- void SetVertexFieldOn();
+ void SetVertexFieldOn();
void SetVertexFieldOnBTh();
Int4 SplitInternalEdgeWithBorderVertices();
void MakeQuadrangles(double costheta);
int SplitElement(int choice);
void MakeQuadTree();
void NewPoints( Triangles &,int KeepBackVertex =1 );
- Int4 InsertNewPoints(Int4 nbvold,Int4 & NbTSwap) ;
+ Int4 InsertNewPoints(Int4 nbvold,Int4 & NbTSwap) ;
void NewPointsOld( Triangles & );
void NewPoints(int KeepBackVertex=1){ NewPoints(*this,KeepBackVertex);}
void ReNumberingTheTriangleBySubDomain(bool justcompress=false);
@@ -858,11 +858,11 @@ public:
Metric MetricAt (const R2 &) const;
GeometricalEdge * ProjectOnCurve( Edge & AB, Vertex & A, Vertex & B,Real8 theta,
Vertex & R,VertexOnEdge & BR,VertexOnGeom & GR);
-
-
+
+
void WriteElements(ostream& f,Int4 * reft ,Int4 nbInT) const;
-
+
Int4 Number(const Triangle & t) const { return &t - triangles;}
Int4 Number(const Triangle * t) const { return t - triangles;}
Int4 Number(const Vertex & t) const { return &t - vertices;}
@@ -874,7 +874,7 @@ public:
return t - triangles;
// else return t - OutSidesTriangles;
}
-
+
Vertex * NearestVertex(Icoor1 i,Icoor1 j) ;
Triangle * FindTriangleContening(const I2 & ,Icoor2 [3],Triangle *tstart=0) const;
void Write(const char * filename,const TypeFileMesh type = AutoMesh);
@@ -901,23 +901,23 @@ public:
const double power=1.0,
const int choise=0);
void IntersectGeomMetric(const Real8 err,const int iso);
-
-
+
+
int isCracked() const {return NbCrackedVertices != 0;}
int Crack();
int UnCrack();
-
+
#ifdef DEBUG
- void inline Check();
+ void inline Check();
#endif
#ifdef DRAWING
void Draw() const ;
void InitDraw() const ;
void inquire() ;
#endif
- friend ostream& operator <<(ostream& f, const Triangles & Th);
+ friend ostream& operator <<(ostream& f, const Triangles & Th);
void Write(const char * filename);
- void ConsGeometry(Real8 =-1.0,int *equiedges=0); // construct a geometry if no geo
+ void ConsGeometry(Real8 =-1.0,int *equiedges=0); // construct a geometry if no geo
void FillHoleInMesh() ;
int CrackMesh();
private:
@@ -928,44 +928,44 @@ public:
void Write_nop5(OFortranUnFormattedFile * f,
Int4 &lnop5,Int4 &nef,Int4 &lgpdn,Int4 ndsr) const ;
-
+
}; // End Class Triangles
/////////////////////////////////////////////////////////////////////////////////////
-class Geometry {
+class Geometry {
public:
- int OnDisk;
+ int OnDisk;
Int4 NbRef; // counter of ref on the this class if 0 we can delete
char *name;
Int4 nbvx,nbtx; // nombre max de sommets , de Geometry
Int4 nbv,nbt,nbiv,nbe; // nombre de sommets, de Geometry, de sommets initiaux,
- Int4 NbSubDomains; //
+ Int4 NbSubDomains; //
Int4 NbEquiEdges;
- Int4 NbCrackedEdges;
+ Int4 NbCrackedEdges;
// Int4 nbtf;// de triangle frontiere
Int4 NbOfCurves;
int empty(){return (nbv ==0) && (nbt==0) && (nbe==0) && (NbSubDomains==0); }
- GeometricalVertex * vertices; // data of vertices des sommets
- Triangle * triangles;
+ GeometricalVertex * vertices; // data of vertices des sommets
+ Triangle * triangles;
GeometricalEdge * edges;
QuadTree *quadtree;
GeometricalSubDomain *subdomains;
Curve *curves;
- ~Geometry();
- Geometry(const Geometry & Gh); //Copy Operator
- Geometry(int nbg,const Geometry ** ag); // intersection operator
-
+ ~Geometry();
+ Geometry(const Geometry & Gh); //Copy Operator
+ Geometry(int nbg,const Geometry ** ag); // intersection operator
+
R2 pmin,pmax; // extrema
Real8 coefIcoor; // coef to integer Icoor1;
Real8 MaximalAngleOfCorner;
-
-// end of data
-
+// end of data
+
+
I2 toI2(const R2 & P) const {
return I2( (Icoor1) (coefIcoor*(P.x-pmin.x))
,(Icoor1) (coefIcoor*(P.y-pmin.y)) );}
-
+
Real8 MinimalHmin() {return 2.0/coefIcoor;}
Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);}
void ReadGeometry(const char * ) ;
@@ -980,29 +980,29 @@ public:
const GeometricalVertex & operator[] (Int4 i) const { return vertices[i];};
GeometricalVertex & operator[](Int4 i) { return vertices[i];};
const GeometricalEdge & operator() (Int4 i) const { return edges[i];};
- GeometricalEdge & operator()(Int4 i) { return edges[i];};
+ GeometricalEdge & operator()(Int4 i) { return edges[i];};
Int4 Number(const GeometricalVertex & t) const { return &t - vertices;}
Int4 Number(const GeometricalVertex * t) const { return t - vertices;}
Int4 Number(const GeometricalEdge & t) const { return &t - edges;}
Int4 Number(const GeometricalEdge * t) const { return t - edges;}
Int4 Number(const Curve * c) const { return c - curves;}
-
+
void UnMarkEdges() {
for (Int4 i=0;i<nbe;i++) edges[i].SetUnMark();}
GeometricalEdge * ProjectOnCurve(const Edge & ,Real8,Vertex &,VertexOnGeom &) const ;
GeometricalEdge * Contening(const R2 P, GeometricalEdge * start) const;
- friend ostream& operator <<(ostream& f, const Geometry & Gh);
+ friend ostream& operator <<(ostream& f, const Geometry & Gh);
void Write(const char * filename);
-
+
#ifdef DEBUG
- void inline Check();
+ void inline Check();
#endif
#ifdef DRAWING
void Draw() const ;
void InitDraw() const ;
#endif
-
+
}; // End Class Geometry
/////////////////////////////////////////////////////////////////////////////////////
@@ -1044,18 +1044,18 @@ inline Triangle * Triangle::Quadrangle(Vertex * & v0,Vertex * & v1,Vertex * &
}
inline double Triangle::QualityQuad(int a,int option) const
-{ // first do the logique part
+{ // first do the logique part
double q;
if (!link || aa[a] &4)
q= -1;
else {
Triangle * t = at[a];
- if (t-this<0) q= -1;// because we do 2 times
+ if (t-this<0) q= -1;// because we do 2 times
else if (!t->link ) q= -1;
else if (aa[0] & 16 || aa[1] & 16 || aa[2] & 16 || t->aa[0] & 16 || t->aa[1] & 16 || t->aa[2] & 16 )
q= -1;
- else if(option)
- {
+ else if(option)
+ {
const Vertex & v2 = *ns[VerticesOfTriangularEdge[a][0]];
const Vertex & v0 = *ns[VerticesOfTriangularEdge[a][1]];
const Vertex & v1 = *ns[OppositeEdge[a]];
@@ -1069,43 +1069,43 @@ inline double Triangle::QualityQuad(int a,int option) const
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 Edge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
- {
+ {
*this = Th.edges[i];
- v[0] = ThNew.vertices + Th.Number(v[0]);
+ v[0] = ThNew.vertices + Th.Number(v[0]);
v[1] = ThNew.vertices + Th.Number(v[1]);
- if (on)
+ if (on)
on = ThNew.Gh.edges+Th.Gh.Number(on);
if (adj[0]) adj[0] = ThNew.edges + Th.Number(adj[0]);
if (adj[1]) adj[1] = ThNew.edges + Th.Number(adj[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]);
+ 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]);
}
-
+
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);
+ be = GhNew.edges + Gh.Number(be);
+ ee = GhNew.edges + Gh.Number(ee);
+ if(next) next= GhNew.curves + Gh.Number(next);
}
inline void Triangle::Set(const Triangle & rec,const Triangles & Th ,Triangles & ThNew)
- {
+ {
*this = rec;
if ( ns[0] ) ns[0] = ThNew.vertices + Th.Number(ns[0]);
if ( ns[1] ) ns[1] = ThNew.vertices + Th.Number(ns[1]);
@@ -1117,8 +1117,8 @@ inline void Triangle::Set(const Triangle & rec,const Triangles & Th ,Triangles &
link = ThNew.triangles + Th.Number(link);
}
inline void VertexOnVertex::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
-{
- *this = Th.VertexOnBThVertex[i];
+{
+ *this = Th.VertexOnBThVertex[i];
v = ThNew.vertices + Th.Number(v);
}
@@ -1126,8 +1126,8 @@ inline void SubDomain::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
{
*this = Th.subdomains[i];
assert( head - Th.triangles >=0 && head - Th.triangles < Th.nbt);
- head = ThNew.triangles + Th.Number(head) ;
- assert(edge - Th.edges >=0 && edge - Th.edges < Th.nbe);
+ head = ThNew.triangles + Th.Number(head) ;
+ assert(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)
@@ -1137,68 +1137,69 @@ inline void SubDomain::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
}
inline void VertexOnEdge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
{
- *this = Th.VertexOnBThEdge[i];
+ *this = Th.VertexOnBThEdge[i];
v = ThNew.vertices + Th.Number(v);
}
inline void VertexOnGeom::Set(const VertexOnGeom & rec,const Triangles & Th ,Triangles & ThNew)
{
- *this = rec;
+ *this = rec;
mv = ThNew.vertices + Th.Number(mv);
- if (gv)
+ 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;i<nbv;i++)
+ for ( i=0;i<nbv;i++)
{
vertices[i].vint = 0;
vertices[i].t=0;
}
- for ( i=0;i<nbt;i++)
+ for ( i=0;i<nbt;i++)
triangles[i].SetTriangleContainingTheVertex();
}
inline void Triangles::UnMarkUnSwapTriangle()
{
register Int4 i;
- for ( i=0;i<nbt;i++)
+ for ( i=0;i<nbt;i++)
for(int j=0;j<3;j++)
triangles[i].SetUnMarkUnSwap(j);
}
inline void Triangles::SetVertexFieldOn()
{
- for (Int4 i=0;i<nbv;i++)
+ for (Int4 i=0;i<nbv;i++)
vertices[i].on=0;
- for (Int4 j=0;j<NbVerticesOnGeomVertex;j++ )
+ for (Int4 j=0;j<NbVerticesOnGeomVertex;j++ )
VerticesOnGeomVertex[j].SetOn();
- for (Int4 k=0;k<NbVerticesOnGeomEdge;k++ )
+ for (Int4 k=0;k<NbVerticesOnGeomEdge;k++ )
VerticesOnGeomEdge[k].SetOn();
- }
+ }
inline void Triangles::SetVertexFieldOnBTh()
{
- for (Int4 i=0;i<nbv;i++)
+ for (Int4 i=0;i<nbv;i++)
vertices[i].on=0;
- for (Int4 j=0;j<NbVertexOnBThVertex;j++ )
+ for (Int4 j=0;j<NbVertexOnBThVertex;j++ )
VertexOnBThVertex[j].SetOnBTh();
- for (Int4 k=0;k<NbVertexOnBThEdge;k++ )
+ for (Int4 k=0;k<NbVertexOnBThEdge;k++ )
VertexOnBThEdge[k].SetOnBTh();
-
- }
+
+ }
inline void TriangleAdjacent::SetAdj2(const TriangleAdjacent & ta, int l )
-{ // set du triangle adjacent
+{ // set du triangle adjacent
if(t) {
t->at[a]=ta.t;
t->aa[a]=ta.a|l;}
@@ -1235,18 +1236,18 @@ inline Icoor2 & TriangleAdjacent::det() const
inline TriangleAdjacent Adj(const TriangleAdjacent & a)
{ return a.Adj();}
-inline TriangleAdjacent Next(const TriangleAdjacent & ta)
+inline TriangleAdjacent Next(const TriangleAdjacent & ta)
{ return TriangleAdjacent(ta.t,NextEdge[ta.a]);}
-inline TriangleAdjacent Previous(const TriangleAdjacent & ta)
+inline TriangleAdjacent Previous(const TriangleAdjacent & ta)
{ return TriangleAdjacent(ta.t,PreviousEdge[ta.a]);}
-
-inline void Adj(GeometricalEdge * & on,int &i)
+
+inline void Adj(GeometricalEdge * & on,int &i)
{int j=i;i=on->SensAdj[i];on=on->Adj[j];}
-
+
inline Real4 qualite(const Vertex &va,const Vertex &vb,const Vertex &vc)
{
- Real4 ret;
+ Real4 ret;
I2 ia=va,ib=vb,ic=vc;
I2 ab=ib-ia,bc=ic-ib,ac=ic-ia;
Icoor2 deta=Det(ab,ac);
@@ -1284,7 +1285,7 @@ inline Triangle::Triangle(Vertex *v0,Vertex *v1,Vertex *v2){
if (v0) det=0;
else {
det=-1;
- link=NULL;};
+ link=NULL;};
}
inline Real4 Triangle::qualite()
@@ -1293,14 +1294,14 @@ inline Real4 Triangle::qualite()
}
Int4 inline Vertex::Optim(int i,int koption)
-{
+{
Int4 ret=0;
if ( t && (vint >= 0 ) && (vint <3) )
{
ret = t->Optim(vint,koption);
- if(!i)
+ if(!i)
{
- t =0; // for no future optime
+ t =0; // for no future optime
vint= 0; }
}
return ret;
@@ -1308,8 +1309,8 @@ Int4 inline Vertex::Optim(int i,int koption)
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;
+ 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;}
@@ -1330,12 +1331,12 @@ int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb ,
int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) ;
-// inline bofbof FH
+// inline bofbof FH
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
@@ -1344,18 +1345,18 @@ inline TriangleAdjacent FindTriangleAdjacent(Edge &E)
int k=0;
do { // turn around vertex in direct sens (trigo)
k++;assert(k< 20000);
- // in no crack => ta.EdgeVertex(1) == a otherwise ???
- if (ta.EdgeVertex(1) == a && ta.EdgeVertex(0) == b) return ta; // find
+ // 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
+ if (ta.EdgeVertex(0) == a && ta.EdgeVertex(1) == b) return ta; // find
--ta;
} while (t != (Triangle *)ta);
assert(0);
- return TriangleAdjacent(0,0);// error
+ return TriangleAdjacent(0,0);// error
}
-
+
inline Vertex * TheVertex(Vertex * a) // give a unique vertex with smallest number
-{ // in case on crack in mesh
+{ // in case on crack in mesh
Vertex * r(a), *rr;
Triangle * t = a->t;
int i = a->vint;
@@ -1370,7 +1371,7 @@ inline Vertex * TheVertex(Vertex * a) // give a unique vertex with smallest numb
ta = ta.Adj();
if ((rr=ta.EdgeVertex(1)) < r) r =rr;
--ta;
- } while (t != (Triangle*) ta);
+ } while (t != (Triangle*) ta);
return r;
}
@@ -1391,28 +1392,28 @@ void inline Triangle::checka(Int1 a) {
assert(a < 3 && a >= 0 );
Triangle *t1=this,*t2=at[a];
Int2 a1=a,a2=aa[a]%4;
-
+
assert(a2 < 3 && a2 >= 0 );
if (t2 && ( ((*t1).ns[VerticesOfTriangularEdge[a1][0]] != (*t2).ns[VerticesOfTriangularEdge[a2][1]])
|| ((*t1).ns[VerticesOfTriangularEdge[a1][1]] != (*t2).ns[VerticesOfTriangularEdge[a2][0]])))
{
- if (CurrentTh) cerr << " In Triangles beetween Triangle " << CurrentTh->Number(t1) << " and "
+ if (CurrentTh) cerr << " In Triangles beetween Triangle " << CurrentTh->Number(t1) << " and "
<< CurrentTh->Number(t2) << endl;
cerr << "---- t1="<< t1 << " " << a1 <<", t2="<< t2 << " " << a2 << endl;
- cerr <<"t1="<< t1 << " " << a1 << " " << t1->ns[VerticesOfTriangularEdge[a1][0]]
+ cerr <<"t1="<< t1 << " " << a1 << " " << t1->ns[VerticesOfTriangularEdge[a1][0]]
<< " " << t1->ns[VerticesOfTriangularEdge[a1][1]] <<endl;
if (CurrentTh)
cerr <<"t1="<< t1 << " " << a1 << " " << CurrentTh->Number(t1->ns[VerticesOfTriangularEdge[a1][0]])
<< " " << CurrentTh->Number(t1->ns[VerticesOfTriangularEdge[a1][1]]) <<endl;
- if (t2) cerr <<"t2="<< t2 << " " << a2 << " "
- << t2->ns[VerticesOfTriangularEdge[a2][0]]
+ if (t2) cerr <<"t2="<< t2 << " " << a2 << " "
+ << t2->ns[VerticesOfTriangularEdge[a2][0]]
<< " " << t2->ns[VerticesOfTriangularEdge[a2][1]] <<endl;
if (t2 &&CurrentTh)
- cerr <<"t2="<< t2 << " " << a2 << " "
+ cerr <<"t2="<< t2 << " " << a2 << " "
<< CurrentTh->Number(t2->ns[VerticesOfTriangularEdge[a2][0]])
<< " " << CurrentTh->Number(t2->ns[VerticesOfTriangularEdge[a2][1]]) <<endl;
- assert(0);
- }
+ assert(0);
+ }
if (t2) assert(t1->aa[a1]/4 == t2->aa[a2]/4); // lock compatibite
}
@@ -1431,15 +1432,15 @@ void inline Triangle::check() {
{ if (CurrentTh) cerr << " In Triangles " << CurrentTh->Number(this) << endl;
cerr << " det = " << det << " and " << infv << endl;
MeshError(5);
- }
-
- if (det >=0)
+ }
+
+ if (det >=0)
if( det != (det2=bamg::det(*ns[0],*ns[1],*ns[2])))
{ // penthickness(4);Draw();
- if (CurrentTh) cerr << " In Triangles" << CurrentTh->Number(this)
+ if (CurrentTh) cerr << " In Triangles" << CurrentTh->Number(this)
<< endl;
cerr << *ns[0] << *ns[1] << " " << *ns[2] << " " << endl;
- cerr << " Bug in triangle " << this
+ cerr << " Bug in triangle " << this
<< ":" << det << " != " << det2 << endl;
MeshError(5);
}
@@ -1463,7 +1464,7 @@ void inline Triangle::check() {
-#ifdef DRAWING
+#ifdef DRAWING
extern Real4 xGrafCoef,yGrafCoef,xGrafOffSet,yGrafOffSet; // R2 -> I2 transform
extern R2 Gpmin,Gpmax;
//extern Real8 Gh;
diff --git a/contrib/rtree/rtree.h b/contrib/rtree/rtree.h
index ffd63ed751fc39ae7ee545f23a0cae32dae7b4b1..138de134039a6fea75da30c609323cfd70949d31 100644
--- a/contrib/rtree/rtree.h
+++ b/contrib/rtree/rtree.h
@@ -12,7 +12,7 @@
#define ASSERT assert // RTree uses ASSERT( condition )
#ifndef Min
- #define Min std::min
+ #define Min std::min
#endif //Min
#ifndef Max
#define Max std::max
@@ -48,11 +48,11 @@ class RTFileStream; // File I/O helper class, look below for implementation and
/// Instead of using a callback function for returned results, I recommend and efficient pre-sized, grow-only memory
/// array similar to MFC CArray or STL Vector for returning search query result.
///
-template<class DATATYPE, class ELEMTYPE, int NUMDIMS,
+template<class DATATYPE, class ELEMTYPE, int NUMDIMS,
class ELEMTYPEREAL = ELEMTYPE, int TMAXNODES = 8, int TMINNODES = TMAXNODES / 2>
class RTree
{
-protected:
+protected:
struct Node; // Fwd decl. Used by other internal structs and iterator
@@ -71,19 +71,19 @@ public:
RTree();
virtual ~RTree();
-
+
/// Insert entry
/// \param a_min Min of bounding rect
/// \param a_max Max of bounding rect
/// \param a_dataId Positive Id of data. Maybe zero, but negative numbers not allowed.
void Insert(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMDIMS], const DATATYPE& a_dataId);
-
+
/// Remove entry
/// \param a_min Min of bounding rect
/// \param a_max Max of bounding rect
/// \param a_dataId Positive Id of data. Maybe zero, but negative numbers not allowed.
void Remove(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMDIMS], const DATATYPE& a_dataId);
-
+
/// Find all within search rectangle
/// \param a_min Min of search bounding rect
/// \param a_max Max of search bounding rect
@@ -91,7 +91,7 @@ public:
/// \param a_context User context to pass as parameter to a_resultCallback
/// \return Returns the number of entries found
int Search(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMDIMS], bool a_resultCallback(DATATYPE a_data, void* a_context), void* a_context);
-
+
/// Remove all entries from tree
void RemoveAll();
@@ -103,7 +103,7 @@ public:
/// Load tree contents from stream
bool Load(RTFileStream& a_stream);
-
+
/// Save tree contents to file
bool Save(const char* a_fileName);
/// Save tree contents to stream
@@ -113,21 +113,21 @@ public:
class Iterator
{
private:
-
+
enum { MAX_STACK = 32 }; // Max stack size. Allows almost n^32 where n is number of branches in node
-
+
struct StackElement
{
Node* m_node;
int m_branchIndex;
};
-
+
public:
-
+
Iterator() { Init(); }
~Iterator() { }
-
+
/// Is iterator invalid
bool IsNull() { return (m_tos <= 0); }
@@ -140,7 +140,7 @@ public:
ASSERT(IsNotNull());
StackElement& curTos = m_stack[m_tos - 1];
return curTos.m_node->m_branch[curTos.m_branchIndex].m_data;
- }
+ }
/// Access the current data element. Caller must be sure iterator is not NULL first.
const DATATYPE& operator*() const
@@ -148,7 +148,7 @@ public:
ASSERT(IsNotNull());
StackElement& curTos = m_stack[m_tos - 1];
return curTos.m_node->m_branch[curTos.m_branchIndex].m_data;
- }
+ }
/// Find the next data element
bool operator++() { return FindNextData(); }
@@ -159,7 +159,7 @@ public:
ASSERT(IsNotNull());
StackElement& curTos = m_stack[m_tos - 1];
Branch& curBranch = curTos.m_node->m_branch[curTos.m_branchIndex];
-
+
for(int index = 0; index < NUMDIMS; ++index)
{
a_min[index] = curBranch.m_rect.m_min[index];
@@ -168,7 +168,7 @@ public:
}
private:
-
+
/// Reset iterator
void Init() { m_tos = 0; }
@@ -205,7 +205,7 @@ public:
// Since cur node is not a leaf, push first of next level to get deeper into the tree
Node* nextLevelnode = curTos.m_node->m_branch[curTos.m_branchIndex].m_child;
Push(nextLevelnode, 0);
-
+
// If we pushed on a new leaf, exit as the data is ready at TOS
if(nextLevelnode->IsLeaf())
{
@@ -223,7 +223,7 @@ public:
++m_tos;
ASSERT(m_tos <= MAX_STACK);
}
-
+
/// Pop element off iteration stack (For internal use only)
StackElement& Pop()
{
@@ -234,7 +234,7 @@ public:
StackElement m_stack[MAX_STACK]; ///< Stack as we are doing iteration instead of recursion
int m_tos; ///< Top Of Stack index
-
+
friend class RTree; // Allow hiding of non-public functions while allowing manipulation by logical owner
};
@@ -259,7 +259,7 @@ public:
}
first = first->m_branch[0].m_child;
}
- }
+ }
/// Get Next for iteration
void GetNext(Iterator& a_it) { ++a_it; }
@@ -275,8 +275,8 @@ protected:
/// Minimal bounding rectangle (n-dimensional)
struct Rect
{
- ELEMTYPE m_min[NUMDIMS]; ///< Min dimensions of bounding box
- ELEMTYPE m_max[NUMDIMS]; ///< Max dimensions of bounding box
+ ELEMTYPE m_min[NUMDIMS]; ///< Min dimensions of bounding box
+ ELEMTYPE m_max[NUMDIMS]; ///< Max dimensions of bounding box
};
/// May be data or may be another subtree
@@ -297,12 +297,12 @@ protected:
{
bool IsInternalNode() { return (m_level > 0); } // Not a leaf, but a internal node
bool IsLeaf() { return (m_level == 0); } // A leaf, contains data
-
+
int m_count; ///< Count
int m_level; ///< Leaf is zero, others positive
Branch m_branch[MAXNODES]; ///< Branch
};
-
+
/// A link list of nodes for reinsertion after a delete operation
struct ListNode
{
@@ -325,8 +325,8 @@ protected:
int m_branchCount;
Rect m_coverSplit;
ELEMTYPEREAL m_coverSplitArea;
- };
-
+ };
+
Node* AllocNode();
void FreeNode(Node* a_node);
void InitNode(Node* a_node);
@@ -361,7 +361,7 @@ protected:
bool SaveRec(Node* a_node, RTFileStream& a_stream);
bool LoadRec(Node* a_node, RTFileStream& a_stream);
-
+
Node* m_root; ///< Root of tree
ELEMTYPEREAL m_unitSphereVolume; ///< Unit sphere constant for required number of dimensions
};
@@ -375,7 +375,7 @@ class RTFileStream
public:
-
+
RTFileStream()
{
m_file = NULL;
@@ -464,7 +464,7 @@ RTREE_QUAL::RTree()
3.298509f, 2.550164f, 1.884104f, // Dimension 9,10,11
1.335263f, 0.910629f, 0.599265f, // Dimension 12,13,14
0.381443f, 0.235331f, 0.140981f, // Dimension 15,16,17
- 0.082146f, 0.046622f, 0.025807f, // Dimension 18,19,20
+ 0.082146f, 0.046622f, 0.025807f, // Dimension 18,19,20
};
m_root = AllocNode();
@@ -491,13 +491,13 @@ void RTREE_QUAL::Insert(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMD
#endif //_DEBUG
Rect rect;
-
+
for(int axis=0; axis<NUMDIMS; ++axis)
{
rect.m_min[axis] = a_min[axis];
rect.m_max[axis] = a_max[axis];
}
-
+
InsertRect(&rect, a_dataId, &m_root, 0);
}
@@ -513,7 +513,7 @@ void RTREE_QUAL::Remove(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMD
#endif //_DEBUG
Rect rect;
-
+
for(int axis=0; axis<NUMDIMS; ++axis)
{
rect.m_min[axis] = a_min[axis];
@@ -535,7 +535,7 @@ int RTREE_QUAL::Search(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMDI
#endif //_DEBUG
Rect rect;
-
+
for(int axis=0; axis<NUMDIMS; ++axis)
{
rect.m_min[axis] = a_min[axis];
@@ -556,7 +556,7 @@ int RTREE_QUAL::Count()
{
int count = 0;
CountRec(m_root, count);
-
+
return count;
}
@@ -591,7 +591,7 @@ bool RTREE_QUAL::Load(const char* a_fileName)
}
bool result = Load(stream);
-
+
stream.Close();
return result;
@@ -630,13 +630,13 @@ bool RTREE_QUAL::Load(RTFileStream& a_stream)
bool result = false;
// Test if header was valid and compatible
- if( (dataFileId == _dataFileId)
- && (dataSize == _dataSize)
- && (dataNumDims == _dataNumDims)
- && (dataElemSize == _dataElemSize)
- && (dataElemRealSize == _dataElemRealSize)
- && (dataMaxNodes == _dataMaxNodes)
- && (dataMinNodes == _dataMinNodes)
+ if( (dataFileId == _dataFileId)
+ && (dataSize == _dataSize)
+ && (dataNumDims == _dataNumDims)
+ && (dataElemSize == _dataElemSize)
+ && (dataElemRealSize == _dataElemRealSize)
+ && (dataMaxNodes == _dataMaxNodes)
+ && (dataMinNodes == _dataMinNodes)
)
{
// Recursively load tree
@@ -722,7 +722,7 @@ bool RTREE_QUAL::Save(RTFileStream& a_stream)
// Recursively save tree
bool result = SaveRec(m_root, a_stream);
-
+
return result;
}
@@ -799,7 +799,7 @@ void RTREE_QUAL::RemoveAllRec(Node* a_node)
RemoveAllRec(a_node->m_branch[index].m_child);
}
}
- FreeNode(a_node);
+ FreeNode(a_node);
}
@@ -941,7 +941,7 @@ bool RTREE_QUAL::InsertRect(Rect* a_rect, const DATATYPE& a_id, Node** a_root, i
{
ASSERT(a_rect->m_min[index] <= a_rect->m_max[index]);
}
-#endif //_DEBUG
+#endif //_DEBUG
Node* newRoot;
Node* newNode;
@@ -970,11 +970,11 @@ RTREE_TEMPLATE
typename RTREE_QUAL::Rect RTREE_QUAL::NodeCover(Node* a_node)
{
ASSERT(a_node);
-
+
int firstTime = true;
Rect rect;
InitRect(&rect);
-
+
for(int index = 0; index < a_node->m_count; ++index)
{
if(firstTime)
@@ -987,7 +987,7 @@ typename RTREE_QUAL::Rect RTREE_QUAL::NodeCover(Node* a_node)
rect = CombineRect(&rect, &(a_node->m_branch[index].m_rect));
}
}
-
+
return rect;
}
@@ -1012,7 +1012,7 @@ bool RTREE_QUAL::AddBranch(Branch* a_branch, Node* a_node, Node** a_newNode)
else
{
ASSERT(a_newNode);
-
+
SplitNode(a_node, a_branch, a_newNode);
return true;
}
@@ -1029,7 +1029,7 @@ void RTREE_QUAL::DisconnectBranch(Node* a_node, int a_index)
// Remove element by swapping with the last element to prevent gaps in array
a_node->m_branch[a_index] = a_node->m_branch[a_node->m_count - 1];
-
+
--a_node->m_count;
}
@@ -1043,7 +1043,7 @@ RTREE_TEMPLATE
int RTREE_QUAL::PickBranch(Rect* a_rect, Node* a_node)
{
ASSERT(a_rect && a_node);
-
+
bool firstTime = true;
ELEMTYPEREAL increase;
ELEMTYPEREAL bestIncr = (ELEMTYPEREAL)-1;
@@ -1121,7 +1121,7 @@ void RTREE_QUAL::SplitNode(Node* a_node, Branch* a_branch, Node** a_newNode)
*a_newNode = AllocNode();
(*a_newNode)->m_level = a_node->m_level = level;
LoadNodes(a_node, *a_newNode, parVars);
-
+
ASSERT((a_node->m_count + (*a_newNode)->m_count) == parVars->m_total);
}
@@ -1131,16 +1131,16 @@ RTREE_TEMPLATE
ELEMTYPEREAL RTREE_QUAL::RectVolume(Rect* a_rect)
{
ASSERT(a_rect);
-
+
ELEMTYPEREAL volume = (ELEMTYPEREAL)1;
for(int index=0; index<NUMDIMS; ++index)
{
volume *= a_rect->m_max[index] - a_rect->m_min[index];
}
-
+
ASSERT(volume >= (ELEMTYPEREAL)0);
-
+
return volume;
}
@@ -1150,18 +1150,18 @@ RTREE_TEMPLATE
ELEMTYPEREAL RTREE_QUAL::RectSphericalVolume(Rect* a_rect)
{
ASSERT(a_rect);
-
+
ELEMTYPEREAL sumOfSquares = (ELEMTYPEREAL)0;
ELEMTYPEREAL radius;
- for(int index=0; index < NUMDIMS; ++index)
+ for(int index=0; index < NUMDIMS; ++index)
{
ELEMTYPEREAL halfExtent = ((ELEMTYPEREAL)a_rect->m_max[index] - (ELEMTYPEREAL)a_rect->m_min[index]) * 0.5f;
sumOfSquares += halfExtent * halfExtent;
}
radius = (ELEMTYPEREAL)sqrt(sumOfSquares);
-
+
// Pow maybe slow, so test for common dims like 2,3 and just use x*x, x*x*x.
if(NUMDIMS == 3)
{
@@ -1186,7 +1186,7 @@ ELEMTYPEREAL RTREE_QUAL::CalcRectVolume(Rect* a_rect)
return RectSphericalVolume(a_rect); // Slower but helps certain merge cases
#else // RTREE_USE_SPHERICAL_VOLUME
return RectVolume(a_rect); // Faster but can cause poor merges
-#endif // RTREE_USE_SPHERICAL_VOLUME
+#endif // RTREE_USE_SPHERICAL_VOLUME
}
@@ -1198,7 +1198,7 @@ void RTREE_QUAL::GetBranches(Node* a_node, Branch* a_branch, PartitionVars* a_pa
ASSERT(a_branch);
ASSERT(a_node->m_count == MAXNODES);
-
+
// Load the branch buffer
for(int index=0; index < MAXNODES; ++index)
{
@@ -1234,10 +1234,10 @@ RTREE_TEMPLATE
void RTREE_QUAL::ChoosePartition(PartitionVars* a_parVars, int a_minFill)
{
ASSERT(a_parVars);
-
+
ELEMTYPEREAL biggestDiff;
int group, chosen, betterGroup;
-
+
InitParVars(a_parVars, a_parVars->m_branchCount, a_minFill);
PickSeeds(a_parVars);
@@ -1303,7 +1303,7 @@ void RTREE_QUAL::ChoosePartition(PartitionVars* a_parVars, int a_minFill)
}
ASSERT((a_parVars->m_count[0] + a_parVars->m_count[1]) == a_parVars->m_total);
- ASSERT((a_parVars->m_count[0] >= a_parVars->m_minFill) &&
+ ASSERT((a_parVars->m_count[0] >= a_parVars->m_minFill) &&
(a_parVars->m_count[1] >= a_parVars->m_minFill));
}
@@ -1319,7 +1319,7 @@ void RTREE_QUAL::LoadNodes(Node* a_nodeA, Node* a_nodeB, PartitionVars* a_parVar
for(int index=0; index < a_parVars->m_total; ++index)
{
ASSERT(a_parVars->m_partition[index] == 0 || a_parVars->m_partition[index] == 1);
-
+
if(a_parVars->m_partition[index] == 0)
{
AddBranch(&a_parVars->m_branchBuf[index], a_nodeA, NULL);
@@ -1353,7 +1353,7 @@ void RTREE_QUAL::InitParVars(PartitionVars* a_parVars, int a_maxRects, int a_min
RTREE_TEMPLATE
void RTREE_QUAL::PickSeeds(PartitionVars* a_parVars)
{
- int seed0, seed1;
+ int seed0 = 0, seed1 = 0;
ELEMTYPEREAL worst, waste;
ELEMTYPEREAL area[MAXNODES+1];
@@ -1433,19 +1433,19 @@ bool RTREE_QUAL::RemoveRect(Rect* a_rect, const DATATYPE& a_id, Node** a_root)
a_root,
tempNode->m_level);
}
-
+
ListNode* remLNode = reInsertList;
reInsertList = reInsertList->m_next;
-
+
FreeNode(remLNode->m_node);
FreeListNode(remLNode);
}
-
+
// Check for redundant root (not leaf, 1 child) and eliminate
if((*a_root)->m_count == 1 && (*a_root)->IsInternalNode())
{
tempNode = (*a_root)->m_branch[0].m_child;
-
+
ASSERT(tempNode);
FreeNode(*a_root);
*a_root = tempNode;
@@ -1569,9 +1569,9 @@ bool RTREE_QUAL::Search(Node* a_node, Rect* a_rect, int& a_foundCount, bool a_re
if(Overlap(a_rect, &a_node->m_branch[index].m_rect))
{
DATATYPE& id = a_node->m_branch[index].m_data;
-
+
// NOTE: There are different ways to return results. Here's where to modify
- if(&a_resultCallback)
+ //if(&a_resultCallback)
{
++a_foundCount;
if(!a_resultCallback(id, a_context))