![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 3230
- Timestamp:
- 03/09/10 08:46:44 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 4 edited
-
Mesh2.h (modified) (26 diffs)
-
meshtype.h (modified) (1 diff)
-
objects/Triangles.cpp (modified) (4 diffs)
-
objects/Vertex.cpp (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/Mesh2.h
r3148 r3230 23 23 /*CLASS: DoubleAndInt4 {{{1*/ 24 24 class DoubleAndInt4 { 25 public: double q; Int4 i3j; 26 int operator<(DoubleAndInt4 a){return q > a.q;} 25 //class used by Triangles::MakeQuadrangles 26 27 public: 28 double q; 29 Int4 i3j; 30 31 //Operators 32 int operator<(DoubleAndInt4 a){return q > a.q;} 27 33 }; 28 34 /*}}}1*/ … … 31 37 private: 32 38 Icoor1 dir; 33 public: 39 40 public: 41 //Methods 34 42 Direction(): dir(MaxICoor){}; // no direction set 35 Direction(Icoor1 i,Icoor1 j) { Icoor2 n2 = 2*(Abs(i)+Abs(j)); 36 Icoor2 r = MaxICoor* (Icoor2) i; 43 Direction(Icoor1 i,Icoor1 j) { 44 Icoor2 n2 = 2*(Abs(i)+Abs(j)); 45 Icoor2 r = MaxICoor* (Icoor2) i; 37 46 Icoor1 r1 = (Icoor1) (2*(r/ n2)); // odd number 38 dir = (j>0) ? r1 : r1+1; // odd -> j>0 even -> j<0 39 } 40 int sens( Icoor1 i,Icoor1 j) { int r =1; 47 dir = (j>0) ? r1 : r1+1; // odd-> j>0 even-> j<0 48 } 49 int sens(Icoor1 i,Icoor1 j) { 50 int r =1; 41 51 if (dir!= MaxICoor) { 42 52 Icoor2 x(dir/2),y1(MaxICoor/2-Abs(x)),y(dir%2?-y1:y1); 43 r = (x*i + y*j) >=0;} 44 return r;} 53 r = (x*i + y*j) >=0; 54 } 55 return r; 56 } 45 57 }; 46 58 /*}}}1*/ 47 59 /*CLASS: Vertex{{{1*/ 48 60 class Vertex { 49 public: 50 I2 i; // allow to use i.x, and i.y in long int (beware of scale and centering) 51 R2 r; // allow to use r.x, and r.y in double 61 62 public: 63 I2 i; // integer coordinates 64 R2 r; // real coordinates 52 65 Metric m; 53 66 Int4 ReferenceNumber; 54 67 Direction DirOfSearch; 68 Int1 vint; // the vertex number in triangle; varies between 0 and 2 in t 55 69 union { 56 Triangle* t; // one triangle which containedthe vertex70 Triangle* t; // one triangle which is containing the vertex 57 71 Int4 color; 58 Vertex * to;// use in geometry Vertex to now the Mesh Vertex associed59 VertexOnGeom * onGeometry; // if vint8; // set with Triangles::SetVertexFieldOn()60 Vertex * onBackgroundVertex;// if vint == 16 on Background vertex Triangles::SetVertexFieldOnBTh()61 VertexOnEdge * onBackgroundEdge;// if vint == 32 on Background edge72 Vertex* to; // use in geometry Vertex to now the Mesh Vertex associed 73 VertexOnGeom* onGeometry; // if vint == 8; // set with Triangles::SetVertexFieldOn() 74 Vertex* onBackgroundVertex; // if vint == 16 on Background vertex Triangles::SetVertexFieldOnBTh() 75 VertexOnEdge* onBackgroundEdge; // if vint == 32 on Background edge 62 76 }; 63 Int1 vint; // the vertex number in triangle; varies between 0 and 2 in t 64 operator I2 () const {return i;} // operator de cast 65 operator const R2 & () const {return r;}// operator de cast 66 // operator R2 & () {return r;}// operator de cast 67 Real8 operator()(R2 x) const { return m(x);} 68 operator Metric () const {return m;}// operator de cast 69 Int4 Optim(int = 1,int =0); 70 // Vertex(){} 71 // ~Vertex(){} 72 Real8 Smoothing(Triangles & ,const Triangles & ,Triangle * & ,Real8 =1); 73 void MetricFromHessian(const double Hxx,const double Hyx, const double Hyy, const double smin,const double smax,const double s,const double err,BamgOpts* bamgopts); 74 void Echo(); 75 int ref() const { return ReferenceNumber;} 76 77 inline void Set(const Vertex & rec,const Triangles &,Triangles &); 77 78 //Operators 79 operator I2() const {return i;} // Cast operator 80 operator const R2 & () const {return r;} // Cast operator 81 operator Metric () const {return m;} // Cast operator 82 Real8 operator()(R2 x) const { return m(x);} // Get x in the metric m 83 84 //methods (No constructor and no destructors...) 85 Int4 Optim(int =1,int =0); 86 Real8 Smoothing(Triangles & ,const Triangles & ,Triangle * & ,Real8 =1); 87 void MetricFromHessian(const double Hxx,const double Hyx, const double Hyy, const double smin,const double smax,const double s,const double err,BamgOpts* bamgopts); 88 void Echo(); 89 int ref() const { return ReferenceNumber;} 90 91 //inline functions 92 inline void Set(const Vertex & rec,const Triangles &,Triangles &); 78 93 }; 79 94 /*}}}1*/ … … 87 102 int a; //Edge number 88 103 104 //Constructors 105 TriangleAdjacent() {}; 89 106 TriangleAdjacent(Triangle * tt,int aa): t(tt),a(aa &3) {}; 90 TriangleAdjacent() {}; 91 107 108 //Operators 92 109 operator Triangle * () const {return t;} 93 110 operator Triangle & () const {return *t;} 94 111 operator int() const {return a;} 95 TriangleAdjacent & operator++(){ 96 a= NextEdge[a]; 97 return *this; 98 } 99 TriangleAdjacent operator--(){ 100 a= PreviousEdge[a]; 101 return *this; 102 } 112 TriangleAdjacent & operator++(){ a= NextEdge[a]; return *this; } 113 TriangleAdjacent operator--(){ a= PreviousEdge[a]; return *this; } 114 115 //Methods 116 int swap(); 117 118 //Inline methods 103 119 inline TriangleAdjacent Adj() const ; 104 int swap();105 120 inline void SetAdj2(const TriangleAdjacent& , int =0); 106 121 inline Vertex * EdgeVertex(const int &) const ; … … 118 133 /*CLASS: Edge{{{1*/ 119 134 class Edge { 120 public: 121 Vertex * v[2]; 122 Int4 ref; 123 GeometricalEdge* onGeometry; 124 Vertex & operator[](int i){return *v[i];}; 125 Vertex * operator()(int i){return v[i];}; 126 127 void ReNumbering(Vertex *vb,Vertex *ve, Int4 *renu) 128 { 129 if (v[0] >=vb && v[0] <ve) v[0] = vb + renu[v[0]-vb]; 130 if (v[1] >=vb && v[1] <ve) v[1] = vb + renu[v[1]-vb]; 131 } 132 133 const Vertex & operator[](int i) const { return *v[i];}; 134 R2 operator()(double t) const; // return the point 135 // on the curve edge a t in [0:1] 136 Edge * adj[2]; // the 2 adj edges if on the same curve 137 138 int Intersection(const Edge & e) const { 139 if (!(adj[0]==&e || adj[1]==&e)){ 140 throw ErrorException(__FUNCT__,exprintf("Intersection bug")); 141 } 142 if (adj[0]!=&e && adj[1]!=&e){ 143 throw ErrorException(__FUNCT__,exprintf("adj[0]!=&e && adj[1]!=&e")); 144 } 145 return adj[0]==&e ? 0 : 1; 146 } 147 Real8 MetricLength() const ; 148 inline void Set(const Triangles &,Int4,Triangles &); 135 136 public: 137 Vertex* v[2]; 138 Int4 ref; 139 GeometricalEdge* onGeometry; 140 Edge* adj[2]; // the 2 adj edges if on the same curve 141 142 //Operators 143 Vertex & operator[](int i){return *v[i];}; 144 Vertex * operator()(int i){return v[i];}; 145 R2 operator()(double t) const; // return the point 146 const Vertex & operator[](int i) const { return *v[i];}; 147 148 //Methods 149 void ReNumbering(Vertex *vb,Vertex *ve, Int4 *renu){ 150 if (v[0] >=vb && v[0] <ve) v[0] = vb + renu[v[0]-vb]; 151 if (v[1] >=vb && v[1] <ve) v[1] = vb + renu[v[1]-vb]; 152 } 153 int Intersection(const Edge & e) const { 154 if (!(adj[0]==&e || adj[1]==&e)){ 155 throw ErrorException(__FUNCT__,exprintf("Intersection bug")); 156 } 157 if (adj[0]!=&e && adj[1]!=&e){ 158 throw ErrorException(__FUNCT__,exprintf("adj[0]!=&e && adj[1]!=&e")); 159 } 160 return adj[0]==&e ? 0 : 1; 161 } 162 163 //Inline methods 164 inline void Set(const Triangles &,Int4,Triangles &); 149 165 150 166 }; … … 152 168 inline TriangleAdjacent FindTriangleAdjacent(Edge &E); 153 169 /*CLASS: GeometricalVertex{{{1*/ 154 class GeometricalVertex :public Vertex { 155 int cas; 156 friend class Geometry; 157 GeometricalVertex* link; // link all the same GeometricalVertex circular (Crack) 158 public: 159 int Corner() const {return cas&4;} 160 int Required()const {return cas&6;}// a corner is required 161 void SetCorner(){ cas |= 4;} 162 void SetRequired(){ cas |= 2;} 163 void Set(){cas=0;} 164 GeometricalVertex() :cas(0), link(this) {}; 165 166 GeometricalVertex* The() { 167 if (!link){ 168 throw ErrorException(__FUNCT__,exprintf("!link")); 169 } 170 return link; 171 }// return a unique vertex 172 173 int IsThe() const { return link == this;} 174 175 inline void Set(const GeometricalVertex & rec,const Geometry & Gh ,const Geometry & GhNew); 170 class GeometricalVertex : public Vertex { 171 172 public: 173 friend class Geometry; 174 175 int cas; 176 GeometricalVertex* link; // link all the same GeometricalVertex circular (Crack) 177 178 //Methods 179 int Corner() const {return cas&4;} 180 int Required()const {return cas&6;}// a corner is required 181 int IsThe() const { return link == this;} 182 void SetCorner(){ cas |= 4;} 183 void SetRequired(){ cas |= 2;} 184 void Set(){cas=0;} 185 GeometricalVertex() :cas(0), link(this) {}; 186 GeometricalVertex* The() { 187 if (!link){ throw ErrorException(__FUNCT__,exprintf("!link"));} 188 return link; 189 }// return a unique vertex 190 191 //Inline methods 192 inline void Set(const GeometricalVertex & rec,const Geometry & Gh ,const Geometry & GhNew); 176 193 }; 177 194 /*}}}1*/ 178 195 /*CLASS: GeometricalEdge{{{1*/ 179 196 class GeometricalEdge { 197 180 198 public: 181 199 GeometricalVertex* v[2]; 182 200 Int4 ref; 183 201 Int4 CurveNumber; 184 R2 tg[2]; // the 2 tangentes 185 // if tg[0] =0 => no continuity 202 R2 tg[2]; // the 2 tangentes (tg[0] =0 => no continuity) 186 203 GeometricalEdge* Adj[2]; 187 204 int DirAdj[2]; 188 205 int flag ; 189 206 GeometricalEdge* link; // if Cracked() or Equi() 190 // end of data 191 192 GeometricalVertex & operator[](int i){return *v[i];};207 208 //Operators 209 GeometricalVertex & operator[](int i){return *v[i];}; 193 210 const GeometricalVertex & operator[](int i) const { return *v[i];}; 194 GeometricalVertex * operator()(int i){return v[i];};195 // inline void Set(const Geometry &,Int4,Geometry &); 196 211 GeometricalVertex * operator()(int i){return v[i];}; 212 213 //Methods 197 214 R2 F(Real8 theta) const ; // parametrization of the curve edge 198 215 Real8 R1tg(Real8 theta,R2 &t) const ; // 1/radius of curvature + tangente 199 int Cracked() const {return flag & 1;}200 int Dup() const { return flag & 32;}201 int Equi()const {return flag & 2;}202 int ReverseEqui()const {return flag & 128;}203 int TgA()const {return flag &4;}204 int TgB()const {return flag &8;}205 int Tg(int i) const{return i==0 ? TgA() : TgB();}206 int Mark()const {return flag &16;}207 int Required() { return flag & 64;}216 int Cracked() const {return flag & 1;} 217 int Dup() const { return flag & 32;} 218 int Equi()const {return flag & 2;} 219 int ReverseEqui()const {return flag & 128;} 220 int TgA()const {return flag &4;} 221 int TgB()const {return flag &8;} 222 int Tg(int i) const{return i==0 ? TgA() : TgB();} 223 int Mark()const {return flag &16;} 224 int Required() { return flag & 64;} 208 225 void SetCracked() { flag |= 1;} 209 226 void SetDup() { flag |= 32;} // not a real edge … … 216 233 void SetReverseEqui() {flag |= 128;} 217 234 235 //Inline methods 218 236 inline void Set(const GeometricalEdge & rec,const Geometry & Th ,Geometry & ThNew); 219 220 237 }; 221 238 /*}}}1*/ 222 239 /*CLASS: Curve{{{1*/ 223 class Curve {public: 224 GeometricalEdge * be,*ee; // begin et end edge 225 int kb,ke; // begin vetex and end vertex 226 Curve *next; // next curve equi to this 227 bool master; // true => of equi curve point on this curve 228 inline void Set(const Curve & rec,const Geometry & Th ,Geometry & ThNew); 229 Curve() : be(0),ee(0),kb(0),ke(0),next(0),master(true) {} 230 void Reverse() { Exchange(be,ee); Exchange(kb,ke);} // revese the sens of the curse 240 class Curve { 241 public: 242 GeometricalEdge * be,*ee; // begin et end edge 243 int kb,ke; // begin vetex and end vertex 244 Curve *next; // next curve equi to this 245 bool master; // true => of equi curve point on this curve 246 247 //Methods 248 Curve() : be(0),ee(0),kb(0),ke(0),next(0),master(true) {} 249 void Reverse() { Exchange(be,ee); Exchange(kb,ke);} // revese the sens of the curse 250 251 //Inline methods 252 inline void Set(const Curve & rec,const Geometry & Th ,Geometry & ThNew); 231 253 }; 232 254 /*}}}1*/ 233 255 /*CLASS: Triangle{{{1*/ 234 256 class Triangle { 257 235 258 friend class TriangleAdjacent; 236 259 237 private: // les arete sont opposes a un sommet 238 Vertex* TriaVertices[3]; // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer 239 Triangle* TriaAdjTriangles[3]; // 3 pointers toward the adjacent triangles 240 Int1 TriaAdjSharedEdge[3]; // number of the edges in the adjacent triangles the edge number 1 is the edge number TriaAdjSharedEdge[1] in the Adjacent triangle 1 260 private: 261 Vertex* TriaVertices[3]; // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer 262 Triangle* TriaAdjTriangles[3]; // 3 pointers toward the adjacent triangles 263 Int1 TriaAdjSharedEdge[3]; // number of the edges in the adjacent triangles the edge number 1 is the edge number TriaAdjSharedEdge[1] in the Adjacent triangle 1 264 241 265 public: 242 Icoor2 det; // determinant du triangle (2 fois l aire des vertices entieres) 243 union { 244 Triangle * link ; 245 Int4 color; 246 }; 247 void Echo(); 248 void SetDet() { 249 if(TriaVertices[0] && TriaVertices[1] && TriaVertices[2]) det = bamg::det(*TriaVertices[0],*TriaVertices[1],*TriaVertices[2]); 250 else det = -1; } 251 Triangle() {} 252 Triangle(Triangles *Th,Int4 i,Int4 j,Int4 k); 253 Triangle(Vertex *v0,Vertex *v1,Vertex *v2); 254 inline void Set(const Triangle &,const Triangles &,Triangles &); 255 inline int In(Vertex *v) const { return TriaVertices[0]==v || TriaVertices[1]==v || TriaVertices[2]==v ;} 256 TriangleAdjacent FindBoundaryEdge(int i) const; 257 258 void ReNumbering(Triangle *tb,Triangle *te, Int4 *renu){ 259 if (link >=tb && link <te) link = tb + renu[link -tb]; 260 if (TriaAdjTriangles[0] >=tb && TriaAdjTriangles[0] <te) TriaAdjTriangles[0] = tb + renu[TriaAdjTriangles[0]-tb]; 261 if (TriaAdjTriangles[1] >=tb && TriaAdjTriangles[1] <te) TriaAdjTriangles[1] = tb + renu[TriaAdjTriangles[1]-tb]; 262 if (TriaAdjTriangles[2] >=tb && TriaAdjTriangles[2] <te) TriaAdjTriangles[2] = tb + renu[TriaAdjTriangles[2]-tb]; 263 } 264 void ReNumbering(Vertex *vb,Vertex *ve, Int4 *renu){ 265 if (TriaVertices[0] >=vb && TriaVertices[0] <ve) TriaVertices[0] = vb + renu[TriaVertices[0]-vb]; 266 if (TriaVertices[1] >=vb && TriaVertices[1] <ve) TriaVertices[1] = vb + renu[TriaVertices[1]-vb]; 267 if (TriaVertices[2] >=vb && TriaVertices[2] <ve) TriaVertices[2] = vb + renu[TriaVertices[2]-vb]; 268 } 269 270 const Vertex & operator[](int i) const {return *TriaVertices[i];}; 271 Vertex & operator[](int i) {return *TriaVertices[i];}; 272 273 const Vertex * operator()(int i) const {return TriaVertices[i];}; 274 Vertex * & operator()(int i) {return TriaVertices[i];}; 275 276 TriangleAdjacent Adj(int i) const // triangle adjacent + arete 277 { return TriangleAdjacent(TriaAdjTriangles[i],TriaAdjSharedEdge[i]&3);}; 278 279 Triangle * TriangleAdj(int i) const 280 {return TriaAdjTriangles[i&3];} // triangle adjacent + arete 281 Int1 NuEdgeTriangleAdj(int i) const 282 {return TriaAdjSharedEdge[i&3]&3;} // Number of the adjacent edge in adj tria 283 284 inline Real4 qualite() ; 285 286 void SetAdjAdj(Int1 a) 287 { a &= 3; 288 register Triangle *tt=TriaAdjTriangles[a]; 289 TriaAdjSharedEdge [a] &= 55; // remove MarkUnSwap 290 register Int1 aatt = TriaAdjSharedEdge[a] & 3; 291 if(tt){ 292 tt->TriaAdjTriangles[aatt]=this; 293 tt->TriaAdjSharedEdge[aatt]=a + (TriaAdjSharedEdge[a] & 60 ) ;}// Copy all the mark 294 } 295 296 void SetAdj2(Int1 a,Triangle *t,Int1 aat){ 297 //TriaAdjTriangles = pointer toward the adjacent triangles of this 298 //TriaAdjSharedEdge = position of the edge in the triangle (mod 3) 299 TriaAdjTriangles[a]=t; //the adjacent triangle to the edge a is t 300 TriaAdjSharedEdge[a]=aat; //position of the edge in the adjacent triangle 301 //if t!=NULL add adjacent triangle to t (this) 302 if(t) { 303 t->TriaAdjTriangles[aat]=this; 304 t->TriaAdjSharedEdge[aat]=a; 305 } 306 } 307 308 void SetTriangleContainingTheVertex() { 309 if (TriaVertices[0]) (TriaVertices[0]->t=this,TriaVertices[0]->vint=0); 310 if (TriaVertices[1]) (TriaVertices[1]->t=this,TriaVertices[1]->vint=1); 311 if (TriaVertices[2]) (TriaVertices[2]->t=this,TriaVertices[2]->vint=2); 312 } 313 314 int swap(Int2 a1,int=0); 315 Int4 Optim(Int2 a,int =0); 316 317 int Locked(int a)const { return TriaAdjSharedEdge[a]&4;} 318 int Hidden(int a)const { return TriaAdjSharedEdge[a]&16;} 319 int Cracked(int a) const { return TriaAdjSharedEdge[a] & 32;} 320 // for optimisation 321 int GetAllflag(int a){return TriaAdjSharedEdge[a] & 1020;} 322 void SetAllFlag(int a,int f){TriaAdjSharedEdge[a] = (TriaAdjSharedEdge[a] &3) + (1020 & f);} 323 324 void SetHidden(int a){ 325 326 //Get Adjacent Triangle number a 327 register Triangle* t = TriaAdjTriangles[a]; 328 329 //if it exist 330 //C|=D -> C=(C|D) bitwise inclusive OR 331 if(t) t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=16; 332 TriaAdjSharedEdge[a] |= 16; 333 } 334 void SetCracked(int a){ 335 register Triangle * t = TriaAdjTriangles[a]; 336 if(t) t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=32; 337 TriaAdjSharedEdge[a] |= 32; 338 } 339 340 double QualityQuad(int a,int option=1) const; 341 Triangle * Quadrangle(Vertex * & v0,Vertex * & v1,Vertex * & v2,Vertex * & v3) const ; 342 343 void SetLocked(int a){ 344 //mark the edge as on Boundary 345 register Triangle * t = TriaAdjTriangles[a]; 346 t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=4; 347 TriaAdjSharedEdge[a] |= 4; 348 } 349 350 void SetMarkUnSwap(int a){ 351 register Triangle * t = TriaAdjTriangles[a]; 352 t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=8; 353 TriaAdjSharedEdge[a] |=8 ; 354 } 355 356 void SetUnMarkUnSwap(int a){ 357 register Triangle * t = TriaAdjTriangles[a]; 358 t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] &=55; // 23 + 32 359 TriaAdjSharedEdge[a] &=55 ; 360 } 266 Icoor2 det; // determinant du triangle (2 fois l aire des vertices entieres) 267 union { 268 Triangle * link ; 269 Int4 color; 270 }; 271 272 //Constructors/Destructors 273 Triangle() {} 274 Triangle(Triangles *Th,Int4 i,Int4 j,Int4 k); 275 Triangle(Vertex *v0,Vertex *v1,Vertex *v2); 276 277 //Operators 278 const Vertex & operator[](int i) const {return *TriaVertices[i];}; 279 Vertex & operator[](int i) {return *TriaVertices[i];}; 280 const Vertex * operator()(int i) const {return TriaVertices[i];}; 281 Vertex * & operator()(int i) {return TriaVertices[i];}; 282 283 //Methods 284 void Echo(); 285 int swap(Int2 a1,int=0); 286 Int4 Optim(Int2 a,int =0); 287 int Locked(int a)const { return TriaAdjSharedEdge[a]&4;} 288 int Hidden(int a)const { return TriaAdjSharedEdge[a]&16;} 289 int Cracked(int a) const { return TriaAdjSharedEdge[a] & 32;} 290 int GetAllflag(int a){return TriaAdjSharedEdge[a] & 1020;} 291 void SetAllFlag(int a,int f){TriaAdjSharedEdge[a] = (TriaAdjSharedEdge[a] &3) + (1020 & f);} 292 double QualityQuad(int a,int option=1) const; 293 Int1 NuEdgeTriangleAdj(int i) const {return TriaAdjSharedEdge[i&3]&3;} // Number of the adjacent edge in adj tria 294 TriangleAdjacent FindBoundaryEdge(int i) const; 295 TriangleAdjacent Adj(int i) const {return TriangleAdjacent(TriaAdjTriangles[i],TriaAdjSharedEdge[i]&3);}; 296 Triangle* TriangleAdj(int i) const {return TriaAdjTriangles[i&3];} 297 Triangle* Quadrangle(Vertex * & v0,Vertex * & v1,Vertex * & v2,Vertex * & v3) const ; 298 void ReNumbering(Triangle *tb,Triangle *te, Int4 *renu){ 299 if (link >=tb && link <te) link = tb + renu[link -tb]; 300 if (TriaAdjTriangles[0] >=tb && TriaAdjTriangles[0] <te) TriaAdjTriangles[0] = tb + renu[TriaAdjTriangles[0]-tb]; 301 if (TriaAdjTriangles[1] >=tb && TriaAdjTriangles[1] <te) TriaAdjTriangles[1] = tb + renu[TriaAdjTriangles[1]-tb]; 302 if (TriaAdjTriangles[2] >=tb && TriaAdjTriangles[2] <te) TriaAdjTriangles[2] = tb + renu[TriaAdjTriangles[2]-tb]; 303 } 304 void ReNumbering(Vertex *vb,Vertex *ve, Int4 *renu){ 305 if (TriaVertices[0] >=vb && TriaVertices[0] <ve) TriaVertices[0] = vb + renu[TriaVertices[0]-vb]; 306 if (TriaVertices[1] >=vb && TriaVertices[1] <ve) TriaVertices[1] = vb + renu[TriaVertices[1]-vb]; 307 if (TriaVertices[2] >=vb && TriaVertices[2] <ve) TriaVertices[2] = vb + renu[TriaVertices[2]-vb]; 308 } 309 void SetAdjAdj(Int1 a){ 310 a &= 3; 311 register Triangle *tt=TriaAdjTriangles[a]; 312 TriaAdjSharedEdge [a] &= 55; // remove MarkUnSwap 313 register Int1 aatt = TriaAdjSharedEdge[a] & 3; 314 if(tt){ 315 tt->TriaAdjTriangles[aatt]=this; 316 tt->TriaAdjSharedEdge[aatt]=a + (TriaAdjSharedEdge[a] & 60 ) ;}// Copy all the mark 317 } 318 void SetAdj2(Int1 a,Triangle *t,Int1 aat){ 319 TriaAdjTriangles[a]=t; //the adjacent triangle to the edge a is t 320 TriaAdjSharedEdge[a]=aat; //position of the edge in the adjacent triangle 321 if(t) { //if t!=NULL add adjacent triangle to t (this) 322 t->TriaAdjTriangles[aat]=this; 323 t->TriaAdjSharedEdge[aat]=a; 324 } 325 } 326 void SetTriangleContainingTheVertex() { 327 if (TriaVertices[0]) (TriaVertices[0]->t=this,TriaVertices[0]->vint=0); 328 if (TriaVertices[1]) (TriaVertices[1]->t=this,TriaVertices[1]->vint=1); 329 if (TriaVertices[2]) (TriaVertices[2]->t=this,TriaVertices[2]->vint=2); 330 } 331 void SetHidden(int a){ 332 //Get Adjacent Triangle number a 333 register Triangle* t = TriaAdjTriangles[a]; 334 //if it exist 335 //C|=D -> C=(C|D) bitwise inclusive OR 336 if(t) t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=16; 337 TriaAdjSharedEdge[a] |= 16; 338 } 339 void SetCracked(int a){ 340 register Triangle * t = TriaAdjTriangles[a]; 341 if(t) t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=32; 342 TriaAdjSharedEdge[a] |= 32; 343 } 344 345 void SetLocked(int a){ 346 //mark the edge as on Boundary 347 register Triangle * t = TriaAdjTriangles[a]; 348 t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=4; 349 TriaAdjSharedEdge[a] |= 4; 350 } 351 void SetMarkUnSwap(int a){ 352 register Triangle * t = TriaAdjTriangles[a]; 353 t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=8; 354 TriaAdjSharedEdge[a] |=8 ; 355 } 356 void SetUnMarkUnSwap(int a){ 357 register Triangle * t = TriaAdjTriangles[a]; 358 t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] &=55; // 23 + 32 359 TriaAdjSharedEdge[a] &=55 ; 360 } 361 void SetDet() { 362 if(TriaVertices[0] && TriaVertices[1] && TriaVertices[2]) det = bamg::det(*TriaVertices[0],*TriaVertices[1],*TriaVertices[2]); 363 else det = -1; } 364 365 //Inline methods 366 inline Real4 qualite() ; 367 inline void Set(const Triangle &,const Triangles &,Triangles &); 368 inline int In(Vertex *v) const { return TriaVertices[0]==v || TriaVertices[1]==v || TriaVertices[2]==v ;} 369 361 370 362 371 }; // end of Triangle class … … 382 391 Real8 lBegin,lEnd; // length abscisse set in ListofIntersectionTriangles::Length 383 392 int last;// last index in ListofIntersectionTriangles for this Sub seg of edge 393 394 //Methods 384 395 R2 F(Real8 s){ 385 396 Real8 c01=lEnd-lBegin, c0=(lEnd-s)/c01, c1=(s-lBegin)/c01; … … 387 398 throw ErrorException(__FUNCT__,exprintf("lBegin>s || s>lEnd")); 388 399 } 389 return e->F(sBegin*c0+sEnd*c1);} 400 return e->F(sBegin*c0+sEnd*c1); 401 } 390 402 }; 391 403 392 int MaxSize; 393 int Size; 394 Real8 len; 395 int state; 396 IntersectionTriangles * lIntTria; 397 int NbSeg; 398 int MaxNbSeg; 399 SegInterpolation * lSegsI; 400 public: 401 IntersectionTriangles & operator[](int i) {return lIntTria[i];} 402 operator int&() {return Size;} 403 404 ListofIntersectionTriangles(int n=256,int m=16) 405 : MaxSize(n), Size(0), len(-1),state(-1),lIntTria(new IntersectionTriangles[n]) , 406 NbSeg(0), MaxNbSeg(m), lSegsI(new SegInterpolation[m]){ 407 long int verbosity=0; 408 if (verbosity>9) printf(" construct ListofIntersectionTriangles %i %i\n",MaxSize,MaxNbSeg); 409 }; 410 411 ~ListofIntersectionTriangles(){ 412 if (lIntTria) delete [] lIntTria,lIntTria=0; 413 if (lSegsI) delete [] lSegsI,lSegsI=0;} 414 void init(){state=0;len=0;Size=0;} 415 416 int NewItem(Triangle * tt,Real8 d0,Real8 d1,Real8 d2); 417 int NewItem(R2,const Metric & ); 418 void NewSubSeg(GeometricalEdge *e,Real8 s0,Real8 s1){ 404 public: 405 406 int MaxSize; 407 int Size; 408 Real8 len; 409 int state; 410 IntersectionTriangles * lIntTria; 411 int NbSeg; 412 int MaxNbSeg; 413 SegInterpolation * lSegsI; 414 415 //Constructors/Destructors 416 ListofIntersectionTriangles(int n=256,int m=16) 417 : MaxSize(n), Size(0), len(-1),state(-1),lIntTria(new IntersectionTriangles[n]) , 418 NbSeg(0), MaxNbSeg(m), lSegsI(new SegInterpolation[m]){ 419 long int verbosity=0; 420 if (verbosity>9) printf(" construct ListofIntersectionTriangles %i %i\n",MaxSize,MaxNbSeg); 421 }; 422 ~ListofIntersectionTriangles(){ 423 if (lIntTria) delete [] lIntTria,lIntTria=0; 424 if (lSegsI) delete [] lSegsI,lSegsI=0; 425 } 426 427 //Operators 428 IntersectionTriangles & operator[](int i) {return lIntTria[i];} 429 operator int&() {return Size;} 430 431 //Methods 432 void init(){state=0;len=0;Size=0;} 433 int NewItem(Triangle * tt,Real8 d0,Real8 d1,Real8 d2); 434 int NewItem(R2,const Metric & ); 435 void SplitEdge(const Triangles & ,const R2 &,const R2 &,int nbegin=0); 436 Real8 Length(); 437 Int4 NewPoints(Vertex *,Int4 & nbv,Int4 nbvx); 438 void NewSubSeg(GeometricalEdge *e,Real8 s0,Real8 s1){ 419 439 long int verbosity=0; 420 421 440 if (NbSeg>=MaxNbSeg) { 422 441 int mneo= MaxNbSeg; … … 439 458 lSegsI[NbSeg].sEnd=s1; 440 459 NbSeg++; 441 } 442 443 void ReShape() { 460 } 461 void ReShape(){ 444 462 register int newsize = MaxSize*2; 445 463 IntersectionTriangles* nw = new IntersectionTriangles[newsize]; 446 if (!nw){ 447 throw ErrorException(__FUNCT__,exprintf("!nw")); 448 } 464 if (!nw){ throw ErrorException(__FUNCT__,exprintf("!nw"));} 449 465 // recopy 450 466 for (int i=0;i<MaxSize;i++) nw[i] = lIntTria[i]; … … 456 472 } 457 473 458 void SplitEdge(const Triangles & ,const R2 &,const R2 &,int nbegin=0);459 Real8 Length();460 Int4 NewPoints(Vertex *,Int4 & nbv,Int4 nbvx);461 474 }; 462 475 /*}}}1*/ … … 467 480 int sens; // -1 or 1 468 481 Int4 ref; 482 483 //Inline methods 469 484 inline void Set(const GeometricalSubDomain &,const Geometry & ,const Geometry &); 470 485 … … 478 493 int sens; // -1 or 1 479 494 Edge * edge; // to geometrical 495 496 //Inline methods 480 497 inline void Set(const Triangles &,Int4,Triangles &); 481 482 498 }; 483 499 /*}}}1*/ 484 500 /*CLASS: VertexOnGeom{{{1*/ 485 501 class VertexOnGeom{ 486 public: 487 488 Vertex * mv; 489 Real8 abscisse; 490 union{ 491 GeometricalVertex * gv; // if abscisse <0; 492 GeometricalEdge * ge; // if abscisse in [0..1] 493 }; 494 inline void Set(const VertexOnGeom&,const Triangles &,Triangles &); 495 int OnGeomVertex()const {return this? abscisse <0 :0;} 496 int OnGeomEdge() const {return this? abscisse >=0 :0;} 497 VertexOnGeom(): mv(0),abscisse(0){gv=0;} 498 VertexOnGeom(Vertex & m,GeometricalVertex &g) : mv(&m),abscisse(-1){gv=&g;} 499 VertexOnGeom(Vertex & m,GeometricalEdge &g,Real8 s) : mv(&m),abscisse(s){ge=&g;} 500 operator Vertex * () const {return mv;} 501 operator GeometricalVertex * () const {return gv;} 502 operator GeometricalEdge * () const {return ge;} 503 operator const Real8 & () const {return abscisse;} 504 int IsRequiredVertex(){ return this? (( abscisse<0 ? (gv?gv->Required():0):0 )) : 0;} 505 void SetOn(){mv->onGeometry=this;mv->vint=IsVertexOnGeom;} 506 inline void Set(const Triangles &,Int4,Triangles &); 502 503 public: 504 505 Vertex* mv; 506 Real8 abscisse; 507 union{ 508 GeometricalVertex * gv; // if abscisse <0; 509 GeometricalEdge * ge; // if abscisse in [0..1] 510 }; 511 512 //Constructors/Destructors 513 VertexOnGeom(): mv(0),abscisse(0){gv=0;} 514 VertexOnGeom(Vertex & m,GeometricalVertex &g) : mv(&m),abscisse(-1){gv=&g;} 515 VertexOnGeom(Vertex & m,GeometricalEdge &g,Real8 s) : mv(&m),abscisse(s){ge=&g;} 516 517 //Operators 518 operator Vertex * () const {return mv;} 519 operator GeometricalVertex * () const {return gv;} 520 operator GeometricalEdge * () const {return ge;} 521 operator const Real8 & () const {return abscisse;} 522 523 //Methods 524 int OnGeomVertex()const {return this? abscisse <0 :0;} 525 int OnGeomEdge() const {return this? abscisse >=0 :0;} 526 int IsRequiredVertex(){ return this? (( abscisse<0 ? (gv?gv->Required():0):0 )) : 0;} 527 void SetOn(){mv->onGeometry=this;mv->vint=IsVertexOnGeom;} 528 529 //Inline methods 530 inline void Set(const Triangles &,Int4,Triangles &); 531 inline void Set(const VertexOnGeom&,const Triangles &,Triangles &); 507 532 508 533 }; … … 510 535 /*CLASS: VertexOnVertex{{{1*/ 511 536 class VertexOnVertex { 512 public: 513 Vertex * v, *bv; 537 538 public: 539 Vertex* v; 540 Vertex* bv; 541 542 //Constructors 514 543 VertexOnVertex(Vertex * w,Vertex *bw) :v(w),bv(bw){} 515 544 VertexOnVertex() {}; 545 546 //Methods 547 void SetOnBTh(){v->onBackgroundVertex=bv;v->vint=IsVertexOnVertex;} 548 549 //Inline methods 516 550 inline void Set(const Triangles &,Int4,Triangles &); 517 void SetOnBTh(){v->onBackgroundVertex=bv;v->vint=IsVertexOnVertex;}518 551 }; 519 552 /*}}}1*/ 520 553 /*CLASS: VertexOnEdge{{{1*/ 521 554 class VertexOnEdge { 522 public: 523 Vertex * v; 524 Edge * be; 525 Real8 abcisse; 526 VertexOnEdge( Vertex * w, Edge *bw,Real8 s) :v(w),be(bw),abcisse(s) {} 527 VertexOnEdge(){} 528 inline void Set(const Triangles &,Int4,Triangles &); 529 void SetOnBTh(){v->onBackgroundEdge=this;v->vint=IsVertexOnEdge;} 530 Vertex & operator[](int i) const { return (*be)[i];} 531 operator Real8 () const { return abcisse;} 532 operator Vertex * () const { return v;} 533 operator Edge * () const { return be;} 555 556 public: 557 Vertex* v; 558 Edge* be; 559 Real8 abcisse; 560 561 //Constructors 562 VertexOnEdge( Vertex * w, Edge *bw,Real8 s) :v(w),be(bw),abcisse(s) {} 563 VertexOnEdge(){} 564 565 //Operators 566 operator Real8 () const { return abcisse;} 567 operator Vertex* () const { return v;} 568 operator Edge* () const { return be;} 569 Vertex & operator[](int i) const { return (*be)[i];} 570 571 //Methods 572 void SetOnBTh(){v->onBackgroundEdge=this;v->vint=IsVertexOnEdge;} 573 574 //Inline methods 575 inline void Set(const Triangles &,Int4,Triangles &); 534 576 }; 535 577 /*}}}1*/ 536 578 /*CLASS: CrackedEdge{{{1*/ 537 579 class CrackedEdge { 580 538 581 friend class Triangles; 539 582 … … 541 584 friend class Triangles; 542 585 friend class CrackedEdge; 543 Triangle 544 int i; // edge number of in triangle545 Edge *edge; // the 2 edge586 Triangle* t; // edge of triangle t 587 int i; // edge number of in triangle 588 Edge *edge; // the 2 edge 546 589 Vertex *New[2]; // new vertex number 590 591 //Constructors 547 592 CrackedTriangle() : t(0),i(0),edge(0) { New[0]=New[1]=0;} 548 593 CrackedTriangle(Edge * a) : t(0),i(0),edge(a) { New[0]=New[1]=0;} 594 595 //Methods 549 596 void Crack(){ 550 597 Triangle & T(*t); … … 555 602 } 556 603 T(i0) = New[0]; 557 T(i1) = New[1];} 558 void UnCrack(){ 559 Triangle & T(*t); 560 int i0=VerticesOfTriangularEdge[i][0]; 561 int i1=VerticesOfTriangularEdge[i][0]; 562 if (!New[0] && !New[1]){ 563 throw ErrorException(__FUNCT__,exprintf("!New[0] && !New[1]")); 564 } 565 T(i0) = TheVertex(T(i0)); 566 T(i1) = TheVertex(T(i1));} 567 void Set() { 568 TriangleAdjacent ta ( FindTriangleAdjacent(*edge)); 569 t = ta; 570 i = ta; 571 572 New[0] = ta.EdgeVertex(0); 573 New[1] = ta.EdgeVertex(1); 574 // warning the ref 575 } 604 T(i1) = New[1]; 605 } 606 void UnCrack(){ 607 Triangle & T(*t); 608 int i0=VerticesOfTriangularEdge[i][0]; 609 int i1=VerticesOfTriangularEdge[i][0]; 610 if (!New[0] && !New[1]){ 611 throw ErrorException(__FUNCT__,exprintf("!New[0] && !New[1]")); 612 } 613 T(i0) = TheVertex(T(i0)); 614 T(i1) = TheVertex(T(i1)); 615 } 616 void Set() { 617 TriangleAdjacent ta ( FindTriangleAdjacent(*edge)); 618 t = ta; 619 i = ta; 620 621 New[0] = ta.EdgeVertex(0); 622 New[1] = ta.EdgeVertex(1); 623 // warning the ref 624 } 576 625 }; 577 626 578 627 public: 579 CrackedTriangle a,b; 580 CrackedEdge() :a(),b() {} 581 CrackedEdge(Edge * start, Int4 i,Int4 j) : a(start+i),b(start+j) {}; 582 CrackedEdge(Edge * e0, Edge * e1 ) : a(e0),b(e1) {}; 583 584 void Crack() { a.Crack(); b.Crack();} 585 void UnCrack() { a.UnCrack(); b.UnCrack();} 586 void Set() { a.Set(), b.Set();} 628 CrackedTriangle a,b; 629 630 //Constructors 631 CrackedEdge() :a(),b() {} 632 CrackedEdge(Edge * start, Int4 i,Int4 j) : a(start+i),b(start+j) {}; 633 CrackedEdge(Edge * e0, Edge * e1 ) : a(e0),b(e1) {}; 634 635 //Methods 636 void Crack() { a.Crack(); b.Crack();} 637 void UnCrack() { a.UnCrack(); b.UnCrack();} 638 void Set() { a.Set(), b.Set();} 587 639 }; 588 640 /*}}}1*/ … … 593 645 Geometry & Gh; // Geometry 594 646 Triangles & BTh; // Background Mesh Bth==*this =>no background 595 596 Int4 NbRef; // counter of ref on the this class if 0 we can delete 647 Int4 NbRef; // counter of ref on the this class if 0 we can delete 597 648 Int4 nbvx,nbtx; // nombre max de sommets , de triangles 598 599 Int4 nt,nbv,nbt,nbiv,nbe; // nb of legal triangles, nb of vertex, of triangles, 600 // of initial vertices, of edges with reference, 601 Int4 NbOfQuad; // nb of quadrangle 602 603 Int4 NbSubDomains; // 604 Int4 NbOutT; // Nb of oudeside triangle 649 Int4 nt,nbv,nbt,nbiv,nbe; // nb of legal triangles, nb of vertex, of triangles, of initial vertices, of edges with reference, 650 Int4 NbOfQuad; // nb of quadrangle 651 Int4 NbSubDomains; 652 Int4 NbOutT; // Nb of oudeside triangle 605 653 Int4 NbOfTriangleSearchFind; 606 654 Int4 NbOfSwapTriangle; 607 Vertex * vertices; // data of vertices des sommets 608 655 Vertex* vertices; 609 656 Int4 NbVerticesOnGeomVertex; 610 657 VertexOnGeom * VerticesOnGeomVertex; 611 612 658 Int4 NbVerticesOnGeomEdge; 613 659 VertexOnGeom * VerticesOnGeomEdge; 614 615 660 Int4 NbVertexOnBThVertex; 616 661 VertexOnVertex *VertexOnBThVertex; 617 618 662 Int4 NbVertexOnBThEdge; 619 663 VertexOnEdge *VertexOnBThEdge; 620 621 664 Int4 NbCrackedVertices; 622 623 665 Int4 NbCrackedEdges; 624 666 CrackedEdge *CrackedEdges; 625 626 R2 pmin,pmax; // extrema 627 Real8 coefIcoor; // coef to integer Icoor1; 628 629 Triangle * triangles; 630 Edge * edges; 631 667 R2 pmin,pmax; // extrema 668 Real8 coefIcoor; // coef to integer Icoor1; 669 Triangle* triangles; 670 Edge* edges; 632 671 QuadTree *quadtree; 633 Vertex 634 SubDomain 672 Vertex** ordre; 673 SubDomain* subdomains; 635 674 ListofIntersectionTriangles lIntTria; 636 // end of variable 637 638 Triangles(Int4 i);//:BTh(*this),Gh(*new Geometry()){PreInit(i);} 639 640 ~Triangles(); 641 Triangles(const char * ,Real8=-1) ; 675 676 //Constructors/Destructors 642 677 Triangles(BamgMesh* bamgmesh,BamgOpts* bamgopts); 643 678 Triangles(double* index,double* x,double* y,int nods,int nels); 644 679 Triangles(Triangles &,Geometry * pGh=0,Triangles* pBTh=0,Int4 nbvxx=0 ); //copy operator 680 Triangles(const Triangles &,const int *flag,const int *bb,BamgOpts* bamgopts); // truncature 645 681 Triangles(Int4 nbvx,Triangles & BT,BamgOpts* bamgopts,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) { 646 682 try {GeomToTriangles1(nbvx,bamgopts,keepBackVertices);} 647 683 catch(...) { this->~Triangles(); throw; } 648 684 } 649 650 //Genetare mesh from geometry651 685 Triangles(Int4 nbvx,Geometry & G,BamgOpts* bamgopts) :Gh(G),BTh(*this){ 652 686 try { GeomToTriangles0(nbvx,bamgopts);} 653 687 catch(...) { this->~Triangles(); throw; } 654 688 } 655 656 Triangles(Triangles &,Geometry * pGh=0,Triangles* pBTh=0,Int4 nbvxx=0 ); // COPY OPERATEUR 657 Triangles(const Triangles &,const int *flag,const int *bb,BamgOpts* bamgopts); // truncature 658 659 void SetIntCoor(const char * from =0); 660 661 Real8 MinimalHmin() {return 2.0/coefIcoor;} 662 Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);} 689 ~Triangles(); 690 691 //Operators 663 692 const Vertex & operator[] (Int4 i) const { return vertices[i];}; 664 693 Vertex & operator[](Int4 i) { return vertices[i];}; 665 694 const Triangle & operator() (Int4 i) const { return triangles[i];}; 666 695 Triangle & operator()(Int4 i) { return triangles[i];}; 696 697 //Methods 698 void SetIntCoor(const char * from =0); 699 Real8 MinimalHmin() {return 2.0/coefIcoor;} 700 Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);} 667 701 I2 toI2(const R2 & P) const { 668 702 return I2( (Icoor1) (coefIcoor*(P.x-pmin.x)) … … 674 708 void Insert(); 675 709 void ForceBoundary(); 676 //void Heap();677 710 void Renumber(BamgOpts* bamgopts); 678 711 void FindSubDomain(int OutSide=0); 679 712 Int4 ConsRefTriangle(Int4 *) const; 680 713 void ShowHistogram() const; 681 void ShowRegulaty() const; // Add FH avril 2007 682 714 void ShowRegulaty() const; 683 715 void ReMakeTriangleContainingTheVertex(); 684 716 void UnMarkUnSwapTriangle(); … … 691 723 Int4 SplitInternalEdgeWithBorderVertices(); 692 724 void MakeQuadrangles(double costheta); 693 int SplitElement(int choice);725 int SplitElement(int choice); 694 726 void MakeQuadTree(); 695 727 void NewPoints(Triangles &,BamgOpts* bamgopts,int KeepVertices=1); … … 699 731 void SmoothingVertex(int =3,Real8=0.3); 700 732 Metric MetricAt (const R2 &) const; 701 GeometricalEdge * ProjectOnCurve( Edge & AB, Vertex & A, Vertex & B,Real8 theta, 702 Vertex & R,VertexOnEdge & BR,VertexOnGeom & GR); 703 733 GeometricalEdge* ProjectOnCurve( Edge & AB, Vertex & A, Vertex & B,Real8 theta, Vertex & R,VertexOnEdge & BR,VertexOnGeom & GR); 704 734 Int4 Number(const Triangle & t) const { return &t - triangles;} 705 735 Int4 Number(const Triangle * t) const { return t - triangles;} … … 709 739 Int4 Number(const Edge * t) const { return t - edges;} 710 740 Int4 Number2(const Triangle * t) const { 711 // if(t>= triangles && t < triangles + nbt )712 741 return t - triangles; 713 // else return t - OutSidesTriangles; 714 } 715 716 Vertex * NearestVertex(Icoor1 i,Icoor1 j) ; 717 Triangle * FindTriangleContaining(const I2 & ,Icoor2 [3],Triangle *tstart=0) const; 718 742 } 743 Vertex* NearestVertex(Icoor1 i,Icoor1 j) ; 744 Triangle* FindTriangleContaining(const I2 & ,Icoor2 [3],Triangle *tstart=0) const; 719 745 void ReadMesh(double* index,double* x,double* y,int nods,int nels); 720 746 void ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts); 721 747 void WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts); 722 723 748 void ReadMetric(BamgOpts* bamgopts,const Real8 hmin,const Real8 hmax,const Real8 coef); 724 749 void WriteMetric(BamgOpts* bamgopts); … … 727 752 void BuildMetric1(BamgOpts* bamgopts); 728 753 void IntersectGeomMetric(BamgOpts* bamgopts); 729 730 int isCracked() const {return NbCrackedVertices != 0;} 731 int Crack(); 732 int UnCrack(); 733 754 int isCracked() const {return NbCrackedVertices != 0;} 755 int Crack(); 756 int UnCrack(); 734 757 void BuildGeometryFromMesh(BamgOpts* bamgopts=NULL); 735 758 void FillHoleInMesh() ; 736 759 int CrackMesh(); 760 737 761 private: 738 762 void GeomToTriangles1(Int4 nbvx,BamgOpts* bamgopts,int KeepVertices=1);// the real constructor mesh adaption … … 744 768 /*CLASS: Geometry{{{1*/ 745 769 class Geometry { 746 public: 747 Int4 NbRef; // counter of ref on the this class if 0 we can delete748 749 Int4 nbvx,nbtx; // nombre max de sommets , de Geometry750 Int4 nbv,nbt,nbiv,nbe; // n ombre de sommets, de Geometry, de sommets initiaux,770 771 public: 772 Int4 NbRef; // counter of ref on the this class if 0 we can delete 773 Int4 nbvx,nbtx; // maximum number of vertices 774 Int4 nbv,nbt,nbiv,nbe; // number of vertices 751 775 Int4 NbSubDomains; // 752 // Int4 nbtf;// de triangle frontiere753 776 Int4 NbOfCurves; 754 int empty(){return (nbv ==0) && (nbt==0) && (nbe==0) && (NbSubDomains==0); }755 777 GeometricalVertex* vertices; 756 Triangle * triangles; 757 GeometricalEdge * edges; 758 QuadTree *quadtree; 759 GeometricalSubDomain *subdomains; 760 Curve *curves; 761 ~Geometry(); 762 Geometry(const Geometry & Gh); //Copy Operator 763 Geometry(int nbg,const Geometry ** ag); // intersection operator 764 778 Triangle* triangles; 779 GeometricalEdge* edges; 780 QuadTree* quadtree; 781 GeometricalSubDomain* subdomains; 782 Curve* curves; 765 783 R2 pmin,pmax; // extrema 766 784 Real8 coefIcoor; // coef to integer Icoor1; 767 785 Real8 MaximalAngleOfCorner; 768 786 787 //Constructor/Destructors 788 ~Geometry(); 789 Geometry(const Geometry & Gh); //Copy Operator 790 Geometry(int nbg,const Geometry** ag); // intersection operator 791 792 //Operators 793 const GeometricalVertex & operator[] (Int4 i) const { return vertices[i];}; 794 GeometricalVertex & operator[](Int4 i) { return vertices[i];}; 795 const GeometricalEdge & operator() (Int4 i) const { return edges[i];}; 796 GeometricalEdge & operator()(Int4 i) { return edges[i];}; 797 798 //Methods 799 int empty(){return (nbv ==0) && (nbt==0) && (nbe==0) && (NbSubDomains==0); } 769 800 void Echo(); 770 771 801 I2 toI2(const R2 & P) const { 772 802 return I2( (Icoor1) (coefIcoor*(P.x-pmin.x)) 773 803 ,(Icoor1) (coefIcoor*(P.y-pmin.y)) );} 774 775 804 Real8 MinimalHmin() {return 2.0/coefIcoor;} 776 805 Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);} … … 780 809 void AfterRead(); 781 810 Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts) {EmptyGeometry();ReadGeometry(bamggeom,bamgopts);AfterRead();} 782 783 const GeometricalVertex & operator[] (Int4 i) const { return vertices[i];};784 GeometricalVertex & operator[](Int4 i) { return vertices[i];};785 const GeometricalEdge & operator() (Int4 i) const { return edges[i];};786 GeometricalEdge & operator()(Int4 i) { return edges[i];};787 811 Int4 Number(const GeometricalVertex & t) const { return &t - vertices;} 788 812 Int4 Number(const GeometricalVertex * t) const { return t - vertices;} … … 790 814 Int4 Number(const GeometricalEdge * t) const { return t - edges;} 791 815 Int4 Number(const Curve * c) const { return c - curves;} 792 793 void UnMarkEdges() { 794 for (Int4 i=0;i<nbe;i++) edges[i].SetUnMark();} 795 816 void UnMarkEdges() {for (Int4 i=0;i<nbe;i++) edges[i].SetUnMark();} 796 817 GeometricalEdge * ProjectOnCurve(const Edge & ,Real8,Vertex &,VertexOnGeom &) const ; 797 818 GeometricalEdge * Contening(const R2 P, GeometricalEdge * start) const; … … 989 1010 /*}}}1*/ 990 1011 /*INLINE functions of CLASS Triangles{{{1*/ 991 inline Triangles::Triangles(Int4 i) :Gh(*new Geometry()),BTh(*this){PreInit(i);}992 1012 inline void Triangles::ReMakeTriangleContainingTheVertex(){ 993 1013 register Int4 i; … … 1127 1147 1128 1148 } 1129 inline Real8 Edge::MetricLength() const1130 {1131 return LengthInterpole(v[0]->m,v[1]->m,v[1]->r - v[0]->r) ;1132 }1133 1149 1134 1150 /*}}}1*/ -
issm/trunk/src/c/Bamgx/meshtype.h
r2985 r3230 38 38 static const Int2 PreviousVertex[3] = {2,0,1}; 39 39 #if LONG_BIT > 63 40 const Icoor1 MaxICoor = 1073741823; // 2^30-1 40 const Icoor1 MaxICoor = 1073741823; // 2^30-1 =111...111 (29 times) 41 41 #else 42 42 const Icoor1 MaxICoor = 8388608; // 2^23 -
issm/trunk/src/c/Bamgx/objects/Triangles.cpp
r3217 r3230 141 141 } 142 142 /*}}}1*/ 143 /*FUNCTION Triangles::~Triangles(){{{1*/144 Triangles::~Triangles() {145 long int verbosity=2;146 //if(vertices) delete [] vertices; //TEST crash if not commented147 if(edges) delete [] edges;148 if(triangles) delete [] triangles;149 if(quadtree) delete quadtree;150 //if(ordre) delete [] ordre; //TEST crash if not commented151 if( subdomains) delete [] subdomains;152 if (VerticesOnGeomEdge) delete [] VerticesOnGeomEdge;153 if (VerticesOnGeomVertex) delete [] VerticesOnGeomVertex;154 if (VertexOnBThVertex) delete [] VertexOnBThVertex;155 if (VertexOnBThEdge) delete [] VertexOnBThEdge;156 157 if (&Gh) {158 if (Gh.NbRef>0) Gh.NbRef--;159 else if (Gh.NbRef==0) delete &Gh;160 }161 if (&BTh && (&BTh != this)) {162 if (BTh.NbRef>0) BTh.NbRef--;163 else if (BTh.NbRef==0) delete &BTh;164 }165 PreInit(0); // set all to zero166 }167 /*}}}1*/168 143 /*FUNCTION Triangles::Triangles(Triangles & Th,Geometry * pGh,Triangles * pBth,Int4 nbvxx) COPY{{{1*/ 169 144 Triangles::Triangles(Triangles & Th,Geometry * pGh,Triangles * pBth,Int4 nbvxx) … … 231 206 232 207 } 208 /*}}}1*/ 209 /*FUNCTION Triangles::~Triangles(){{{1*/ 210 Triangles::~Triangles() { 211 long int verbosity=2; 212 //if(vertices) delete [] vertices; //TEST crash if not commented 213 if(edges) delete [] edges; 214 if(triangles) delete [] triangles; 215 if(quadtree) delete quadtree; 216 //if(ordre) delete [] ordre; //TEST crash if not commented 217 if( subdomains) delete [] subdomains; 218 if (VerticesOnGeomEdge) delete [] VerticesOnGeomEdge; 219 if (VerticesOnGeomVertex) delete [] VerticesOnGeomVertex; 220 if (VertexOnBThVertex) delete [] VertexOnBThVertex; 221 if (VertexOnBThEdge) delete [] VertexOnBThEdge; 222 223 if (&Gh) { 224 if (Gh.NbRef>0) Gh.NbRef--; 225 else if (Gh.NbRef==0) delete &Gh; 226 } 227 if (&BTh && (&BTh != this)) { 228 if (BTh.NbRef>0) BTh.NbRef--; 229 else if (BTh.NbRef==0) delete &BTh; 230 } 231 PreInit(0); // set all to zero 232 } 233 233 /*}}}1*/ 234 234 … … 2190 2190 TriangleAdjacent ta(0,0); 2191 2191 for (Int4 i = 0; i < nbe; i++){ 2192 2193 //Force edge i 2192 2194 nbswp = ForceEdge(edges[i][0],edges[i][1],ta); 2193 2194 if ( nbswp < 0) k++; 2195 if (nbswp<0) k++; 2195 2196 else Nbswap += nbswp; 2197 2196 2198 if (nbswp) nbfe++; 2197 2199 if ( nbswp < 0 && k < 5){ 2198 throw ErrorException(__FUNCT__,exprintf("Missing Edge %i, v0=%i,v1=%i",i ,Number(edges[i][0]),Number(edges[i][1]))); 2200 for (Int4 j = 0; j < nbe; j++){ 2201 printf("Edge %i: %i %i\n",j,Number(edges[j][0]),Number(edges[j][1])); 2202 } 2203 throw ErrorException(__FUNCT__,exprintf("Missing Edge %i, v0=%i,v1=%i",i,Number(edges[i][0]),Number(edges[i][1]))); 2199 2204 } 2200 2205 if ( nbswp >=0 && edges[i].onGeometry->Cracked()) … … 3154 3159 FindSubDomain(); 3155 3160 3156 // NewPointsOld(*this) ;3157 3161 if (verbose>3) printf(" Inserting internal points\n"); 3158 3162 NewPoints(*this,bamgopts,0) ; -
issm/trunk/src/c/Bamgx/objects/Vertex.cpp
r3126 r3230 20 20 /*FUNCTION Vertex::Smoothing{{{1*/ 21 21 Real8 Vertex::Smoothing(Triangles &Th,const Triangles &BTh,Triangle* &tstart ,Real8 omega){ 22 22 23 register Vertex* s=this; 23 24 Vertex &vP = *s,vPsave=vP; … … 43 44 j = NextEdge[jc]; 44 45 if (k>=2000){ 45 throw ErrorException(__FUNCT__,exprintf("k>=2000 "));46 throw ErrorException(__FUNCT__,exprintf("k>=2000 (Maximum number of iterations reached)")); 46 47 } 47 48 } while ( tbegin != tria);
Note:
See TracChangeset
for help on using the changeset viewer.
