source: issm/trunk/src/c/Bamgx/Mesh2.h@ 2985

Last change on this file since 2985 was 2985, checked in by Mathieu Morlighem, 15 years ago

minor
File size: 40.1 KB
Line 
1#ifndef _MESH2_H_
2#define _MESH2_H_
3
4#include "../objects/objects.h"
5#include "../shared/shared.h"
6#include "../include/macros.h"
7#include "../toolkits/toolkits.h"
8
9//From MeshIo
10#include <cstdio>
11#include <cstring>
12#include <cstdlib>
13#include <cctype>
14#include <stdlib.h>
15#include <math.h>
16
17#include "meshtype.h"
18#include "Metric.h"
19
20using namespace std;
21
22namespace bamg {
23
24 //classes
25 class Geometry;
26 class Triangles;
27 class Triangle;
28 class QuadTree;
29 class GeometricalEdge;
30 class VertexOnGeom;
31 class VertexOnEdge;
32
33 /*CLASS: DoubleAndInt4 {{{1*/
34 class DoubleAndInt4 {
35 public: double q; Int4 i3j;
36 int operator<(DoubleAndInt4 a){return q > a.q;}
37 };
38 /*}}}1*/
39 /*CLASS: Direction{{{1*/
40 class Direction { //
41 private:
42 Icoor1 dir;
43 public:
44 Direction(): dir(MaxICoor){}; // no direction set
45 Direction(Icoor1 i,Icoor1 j) { Icoor2 n2 = 2*(Abs(i)+Abs(j));
46 Icoor2 r = MaxICoor* (Icoor2) i;
47 Icoor1 r1 = (Icoor1) (2*(r/ n2)); // odd number
48 dir = (j>0) ? r1 : r1+1; // odd -> j>0 even -> j<0
49 }
50 int sens( Icoor1 i,Icoor1 j) { int r =1;
51 if (dir!= MaxICoor) {
52 Icoor2 x(dir/2),y1(MaxICoor/2-Abs(x)),y(dir%2?-y1:y1);
53 r = (x*i + y*j) >=0;}
54 return r;}
55 };
56 /*}}}1*/
57 /*CLASS: Vertex{{{1*/
58 class Vertex {
59 public:
60 I2 i; // allow to use i.x, and i.y in long int (beware of scale and centering)
61 R2 r; // allow to use r.x, and r.y in double
62 Metric m;
63 Int4 ReferenceNumber;
64 Direction DirOfSearch;
65 union {
66 Triangle* t; // one triangle which contained the vertex
67 Int4 color;
68 Vertex * to;// use in geometry Vertex to now the Mesh Vertex associed
69 VertexOnGeom * onGeometry; // if vint 8; // set with Triangles::SetVertexFieldOn()
70 Vertex * onBackgroundVertex; // if vint == 16 on Background vertex Triangles::SetVertexFieldOnBTh()
71 VertexOnEdge * onBackgroundEdge; // if vint == 32 on Background edge
72 };
73 Int1 vint; // the vertex number in triangle; varies between 0 and 2 in t
74 operator I2 () const {return i;} // operator de cast
75 operator const R2 & () const {return r;}// operator de cast
76 // operator R2 & () {return r;}// operator de cast
77 Real8 operator()(R2 x) const { return m(x);}
78 operator Metric () const {return m;}// operator de cast
79 Int4 Optim(int = 1,int =0);
80 // Vertex(){}
81 // ~Vertex(){}
82 Real8 Smoothing(Triangles & ,const Triangles & ,Triangle * & ,Real8 =1);
83 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);
84 void Echo();
85 int ref() const { return ReferenceNumber;}
86
87 inline void Set(const Vertex & rec,const Triangles &,Triangles &);
88 };
89 /*}}}1*/
90 inline Vertex* TheVertex(Vertex * a); // for remove crak in mesh
91 double QuadQuality(const Vertex &,const Vertex &,const Vertex &,const Vertex &);
92 /*CLASS: TriangleAdjacent{{{1*/
93 class TriangleAdjacent {
94
95 public:
96 Triangle* t; //pointer toward the triangle
97 int a; //Edge number
98
99 TriangleAdjacent(Triangle * tt,int aa): t(tt),a(aa &3) {};
100 TriangleAdjacent() {};
101
102 operator Triangle * () const {return t;}
103 operator Triangle & () const {return *t;}
104 operator int() const {return a;}
105 TriangleAdjacent & operator++(){
106 a= NextEdge[a];
107 return *this;
108 }
109 TriangleAdjacent operator--(){
110 a= PreviousEdge[a];
111 return *this;
112 }
113 inline TriangleAdjacent Adj() const ;
114 int swap();
115 inline void SetAdj2(const TriangleAdjacent& , int =0);
116 inline Vertex * EdgeVertex(const int &) const ;
117 inline Vertex * OppositeVertex() const ;
118 inline Icoor2 & det() const;
119 inline int Locked() const ;
120 inline int GetAllFlag_UnSwap() const ;
121 inline void SetLock();
122 inline int MarkUnSwap() const;
123 inline void SetMarkUnSwap();
124 inline void SetCracked();
125 inline int Cracked() const ;
126 };
127 /*}}}1*/
128 /*CLASS: Edge{{{1*/
129 class Edge {
130 public:
131 Vertex * v[2];
132 Int4 ref;
133 GeometricalEdge* onGeometry;
134 Vertex & operator[](int i){return *v[i];};
135 Vertex * operator()(int i){return v[i];};
136
137 void ReNumbering(Vertex *vb,Vertex *ve, Int4 *renu)
138 {
139 if (v[0] >=vb && v[0] <ve) v[0] = vb + renu[v[0]-vb];
140 if (v[1] >=vb && v[1] <ve) v[1] = vb + renu[v[1]-vb];
141 }
142
143 const Vertex & operator[](int i) const { return *v[i];};
144 R2 operator()(double t) const; // return the point
145 // on the curve edge a t in [0:1]
146 Edge * adj[2]; // the 2 adj edges if on the same curve
147
148 int Intersection(const Edge & e) const {
149 if (!(adj[0]==&e || adj[1]==&e)){
150 throw ErrorException(__FUNCT__,exprintf("Intersection bug"));
151 }
152 if (adj[0]!=&e && adj[1]!=&e){
153 throw ErrorException(__FUNCT__,exprintf("adj[0]!=&e && adj[1]!=&e"));
154 }
155 return adj[0]==&e ? 0 : 1;
156 }
157 Real8 MetricLength() const ;
158 inline void Set(const Triangles &,Int4,Triangles &);
159
160 };
161 /*}}}1*/
162 inline TriangleAdjacent FindTriangleAdjacent(Edge &E);
163 /*CLASS: GeometricalVertex{{{1*/
164 class GeometricalVertex :public Vertex {
165 int cas;
166 friend class Geometry;
167 GeometricalVertex* link; // link all the same GeometricalVertex circular (Crack)
168 public:
169 int Corner() const {return cas&4;}
170 int Required()const {return cas&6;}// a corner is required
171 void SetCorner(){ cas |= 4;}
172 void SetRequired(){ cas |= 2;}
173 void Set(){cas=0;}
174 GeometricalVertex() :cas(0), link(this) {};
175
176 GeometricalVertex* The() {
177 if (!link){
178 throw ErrorException(__FUNCT__,exprintf("!link"));
179 }
180 return link;
181 }// return a unique vertex
182
183 int IsThe() const { return link == this;}
184
185 inline void Set(const GeometricalVertex & rec,const Geometry & Gh ,const Geometry & GhNew);
186 };
187 /*}}}1*/
188 /*CLASS: GeometricalEdge{{{1*/
189 class GeometricalEdge {
190 public:
191 GeometricalVertex* v[2];
192 Int4 ref;
193 Int4 CurveNumber;
194 R2 tg[2]; // the 2 tangentes
195 // if tg[0] =0 => no continuity
196 GeometricalEdge* Adj[2];
197 int DirAdj[2];
198 int flag ;
199 GeometricalEdge* link; // if Cracked() or Equi()
200 // end of data
201
202 GeometricalVertex & operator[](int i){return *v[i];};
203 const GeometricalVertex & operator[](int i) const { return *v[i];};
204 GeometricalVertex * operator()(int i){return v[i];};
205 // inline void Set(const Geometry &,Int4,Geometry &);
206
207 R2 F(Real8 theta) const ; // parametrization of the curve edge
208 Real8 R1tg(Real8 theta,R2 &t) const ; // 1/radius of curvature + tangente
209 int Cracked() const {return flag & 1;}
210 int Dup() const { return flag & 32;}
211 int Equi()const {return flag & 2;}
212 int ReverseEqui()const {return flag & 128;}
213 int TgA()const {return flag &4;}
214 int TgB()const {return flag &8;}
215 int Tg(int i) const{return i==0 ? TgA() : TgB();}
216 int Mark()const {return flag &16;}
217 int Required() { return flag & 64;}
218 void SetCracked() { flag |= 1;}
219 void SetDup() { flag |= 32;} // not a real edge
220 void SetEqui() { flag |= 2;}
221 void SetTgA() { flag|=4;}
222 void SetTgB() { flag|=8;}
223 void SetMark() { flag|=16;}
224 void SetUnMark() { flag &= 1007 /* 1023-16*/;}
225 void SetRequired() { flag |= 64;}
226 void SetReverseEqui() {flag |= 128;}
227
228 inline void Set(const GeometricalEdge & rec,const Geometry & Th ,Geometry & ThNew);
229
230 };
231 /*}}}1*/
232 /*CLASS: Curve{{{1*/
233 class Curve {public:
234 GeometricalEdge * be,*ee; // begin et end edge
235 int kb,ke; // begin vetex and end vertex
236 Curve *next; // next curve equi to this
237 bool master; // true => of equi curve point on this curve
238 inline void Set(const Curve & rec,const Geometry & Th ,Geometry & ThNew);
239 Curve() : be(0),ee(0),kb(0),ke(0),next(0),master(true) {}
240 void Reverse() { Exchange(be,ee); Exchange(kb,ke);} // revese the sens of the curse
241 };
242 /*}}}1*/
243 /*CLASS: Triangle{{{1*/
244 class Triangle {
245 friend class TriangleAdjacent;
246
247 private: // les arete sont opposes a un sommet
248 Vertex* TriaVertices[3]; // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer
249 Triangle* TriaAdjTriangles[3]; // 3 pointers toward the adjacent triangles
250 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
251 public:
252 Icoor2 det; // determinant du triangle (2 fois l aire des vertices entieres)
253 union {
254 Triangle * link ;
255 Int4 color;
256 };
257 void Echo();
258 void SetDet() {
259 if(TriaVertices[0] && TriaVertices[1] && TriaVertices[2]) det = bamg::det(*TriaVertices[0],*TriaVertices[1],*TriaVertices[2]);
260 else det = -1; }
261 Triangle() {}
262 Triangle(Triangles *Th,Int4 i,Int4 j,Int4 k);
263 Triangle(Vertex *v0,Vertex *v1,Vertex *v2);
264 inline void Set(const Triangle &,const Triangles &,Triangles &);
265 inline int In(Vertex *v) const { return TriaVertices[0]==v || TriaVertices[1]==v || TriaVertices[2]==v ;}
266 TriangleAdjacent FindBoundaryEdge(int i) const;
267
268 void ReNumbering(Triangle *tb,Triangle *te, Int4 *renu){
269 if (link >=tb && link <te) link = tb + renu[link -tb];
270 if (TriaAdjTriangles[0] >=tb && TriaAdjTriangles[0] <te) TriaAdjTriangles[0] = tb + renu[TriaAdjTriangles[0]-tb];
271 if (TriaAdjTriangles[1] >=tb && TriaAdjTriangles[1] <te) TriaAdjTriangles[1] = tb + renu[TriaAdjTriangles[1]-tb];
272 if (TriaAdjTriangles[2] >=tb && TriaAdjTriangles[2] <te) TriaAdjTriangles[2] = tb + renu[TriaAdjTriangles[2]-tb];
273 }
274 void ReNumbering(Vertex *vb,Vertex *ve, Int4 *renu){
275 if (TriaVertices[0] >=vb && TriaVertices[0] <ve) TriaVertices[0] = vb + renu[TriaVertices[0]-vb];
276 if (TriaVertices[1] >=vb && TriaVertices[1] <ve) TriaVertices[1] = vb + renu[TriaVertices[1]-vb];
277 if (TriaVertices[2] >=vb && TriaVertices[2] <ve) TriaVertices[2] = vb + renu[TriaVertices[2]-vb];
278 }
279
280 const Vertex & operator[](int i) const {return *TriaVertices[i];};
281 Vertex & operator[](int i) {return *TriaVertices[i];};
282
283 const Vertex * operator()(int i) const {return TriaVertices[i];};
284 Vertex * & operator()(int i) {return TriaVertices[i];};
285
286 TriangleAdjacent Adj(int i) const // triangle adjacent + arete
287 { return TriangleAdjacent(TriaAdjTriangles[i],TriaAdjSharedEdge[i]&3);};
288
289 Triangle * TriangleAdj(int i) const
290 {return TriaAdjTriangles[i&3];} // triangle adjacent + arete
291 Int1 NuEdgeTriangleAdj(int i) const
292 {return TriaAdjSharedEdge[i&3]&3;} // Number of the adjacent edge in adj tria
293
294 inline Real4 qualite() ;
295
296 void SetAdjAdj(Int1 a)
297 { a &= 3;
298 register Triangle *tt=TriaAdjTriangles[a];
299 TriaAdjSharedEdge [a] &= 55; // remove MarkUnSwap
300 register Int1 aatt = TriaAdjSharedEdge[a] & 3;
301 if(tt){
302 tt->TriaAdjTriangles[aatt]=this;
303 tt->TriaAdjSharedEdge[aatt]=a + (TriaAdjSharedEdge[a] & 60 ) ;}// Copy all the mark
304 }
305
306 void SetAdj2(Int1 a,Triangle *t,Int1 aat){
307 //TriaAdjTriangles = pointer toward the adjacent triangles of this
308 //TriaAdjSharedEdge = position of the edge in the triangle (mod 3)
309 TriaAdjTriangles[a]=t; //the adjacent triangle to the edge a is t
310 TriaAdjSharedEdge[a]=aat; //position of the edge in the adjacent triangle
311 //if t!=NULL add adjacent triangle to t (this)
312 if(t) {
313 t->TriaAdjTriangles[aat]=this;
314 t->TriaAdjSharedEdge[aat]=a;
315 }
316 }
317
318 void SetTriangleContainingTheVertex() {
319 if (TriaVertices[0]) (TriaVertices[0]->t=this,TriaVertices[0]->vint=0);
320 if (TriaVertices[1]) (TriaVertices[1]->t=this,TriaVertices[1]->vint=1);
321 if (TriaVertices[2]) (TriaVertices[2]->t=this,TriaVertices[2]->vint=2);
322 }
323
324 int swap(Int2 a1,int=0);
325 Int4 Optim(Int2 a,int =0);
326
327 int Locked(int a)const { return TriaAdjSharedEdge[a]&4;}
328 int Hidden(int a)const { return TriaAdjSharedEdge[a]&16;}
329 int Cracked(int a) const { return TriaAdjSharedEdge[a] & 32;}
330 // for optimisation
331 int GetAllflag(int a){return TriaAdjSharedEdge[a] & 1020;}
332 void SetAllFlag(int a,int f){TriaAdjSharedEdge[a] = (TriaAdjSharedEdge[a] &3) + (1020 & f);}
333
334 void SetHidden(int a){
335
336 //Get Adjacent Triangle number a
337 register Triangle* t = TriaAdjTriangles[a];
338
339 //if it exist
340 //C|=D -> C=(C|D) bitwise inclusive OR
341 if(t) t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=16;
342 TriaAdjSharedEdge[a] |= 16;
343 }
344 void SetCracked(int a){
345 register Triangle * t = TriaAdjTriangles[a];
346 if(t) t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=32;
347 TriaAdjSharedEdge[a] |= 32;
348 }
349
350 double QualityQuad(int a,int option=1) const;
351 Triangle * Quadrangle(Vertex * & v0,Vertex * & v1,Vertex * & v2,Vertex * & v3) const ;
352
353 void SetLocked(int a){
354 //mark the edge as on Boundary
355 register Triangle * t = TriaAdjTriangles[a];
356 t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=4;
357 TriaAdjSharedEdge[a] |= 4;
358 }
359
360 void SetMarkUnSwap(int a){
361 register Triangle * t = TriaAdjTriangles[a];
362 t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=8;
363 TriaAdjSharedEdge[a] |=8 ;
364 }
365
366 void SetUnMarkUnSwap(int a){
367 register Triangle * t = TriaAdjTriangles[a];
368 t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] &=55; // 23 + 32
369 TriaAdjSharedEdge[a] &=55 ;
370 }
371
372 }; // end of Triangle class
373 /*}}}1*/
374 /*CLASS: ListofIntersectionTriangles{{{1*/
375 class ListofIntersectionTriangles {
376
377 class IntersectionTriangles {
378 public:
379 Triangle *t;
380 Real8 bary[3]; // use if t != 0
381 R2 x;
382 Metric m;
383 Real8 s;// abscisse curviline
384 Real8 sp; // len of the previous seg in m
385 Real8 sn;// len of the next seg in m
386 };
387
388 class SegInterpolation {
389 public:
390 GeometricalEdge * e;
391 Real8 sBegin,sEnd; // abscisse of the seg on edge parameter
392 Real8 lBegin,lEnd; // length abscisse set in ListofIntersectionTriangles::Length
393 int last;// last index in ListofIntersectionTriangles for this Sub seg of edge
394 R2 F(Real8 s){
395 Real8 c01=lEnd-lBegin, c0=(lEnd-s)/c01, c1=(s-lBegin)/c01;
396 if (lBegin>s || s>lEnd){
397 throw ErrorException(__FUNCT__,exprintf("lBegin>s || s>lEnd"));
398 }
399 return e->F(sBegin*c0+sEnd*c1);}
400 };
401
402 int MaxSize; //
403 int Size; //
404 Real8 len; //
405 int state;
406 IntersectionTriangles * lIntTria;
407 int NbSeg;
408 int MaxNbSeg;
409 SegInterpolation * lSegsI;
410 public:
411 IntersectionTriangles & operator[](int i) {return lIntTria[i];}
412 operator int&() {return Size;}
413
414 ListofIntersectionTriangles(int n=256,int m=16)
415 : MaxSize(n), Size(0), len(-1),state(-1),lIntTria(new IntersectionTriangles[n]) ,
416 NbSeg(0), MaxNbSeg(m), lSegsI(new SegInterpolation[m])
417 {
418 long int verbosity=0;
419 if (verbosity>9) printf(" construct ListofIntersectionTriangles %i %i\n",MaxSize,MaxNbSeg);
420 };
421
422 ~ListofIntersectionTriangles(){
423 if (lIntTria) delete [] lIntTria,lIntTria=0;
424 if (lSegsI) delete [] lSegsI,lSegsI=0;}
425 void init(){state=0;len=0;Size=0;}
426
427 int NewItem(Triangle * tt,Real8 d0,Real8 d1,Real8 d2);
428 int NewItem(R2,const Metric & );
429 void NewSubSeg(GeometricalEdge *e,Real8 s0,Real8 s1)
430 {
431 long int verbosity=0;
432
433 if (NbSeg>=MaxNbSeg) {
434 int mneo= MaxNbSeg;
435 MaxNbSeg *= 2;
436 if (verbosity>3){
437 printf(" reshape lSegsI from %i to %i\n",mneo,MaxNbSeg);
438 }
439 SegInterpolation * lEn = new SegInterpolation[MaxNbSeg];
440 if (!lSegsI || NbSeg>=MaxNbSeg){
441 throw ErrorException(__FUNCT__,exprintf("!lSegsI || NbSeg>=MaxNbSeg"));
442 }
443 for (int i=0;i< NbSeg;i++)
444 lEn[i] = lSegsI[MaxNbSeg]; // copy old to new
445 delete [] lSegsI; // remove old
446 lSegsI = lEn;
447 }
448 if (NbSeg)
449 lSegsI[NbSeg-1].last=Size;
450 lSegsI[NbSeg].e=e;
451 lSegsI[NbSeg].sBegin=s0;
452 lSegsI[NbSeg].sEnd=s1;
453 NbSeg++;
454 }
455
456 // void CopyMetric(int i,int j){ lIntTria[j].m=lIntTria[i].m;}
457 // void CopyMetric(const Metric & mm,int j){ lIntTria[j].m=mm;}
458
459 void ReShape() {
460 register int newsize = MaxSize*2;
461 IntersectionTriangles* nw = new IntersectionTriangles[newsize];
462 if (!nw){
463 throw ErrorException(__FUNCT__,exprintf("!nw"));
464 }
465 for (int i=0;i<MaxSize;i++) // recopy
466 nw[i] = lIntTria[i];
467 long int verbosity=0;
468 if(verbosity>3) printf(" ListofIntersectionTriangles ReShape Maxsize %i -> %i\n",MaxSize,MaxNbSeg);
469 MaxSize = newsize;
470 delete [] lIntTria;// remove old
471 lIntTria = nw; // copy pointer
472 }
473
474 void SplitEdge(const Triangles & ,const R2 &,const R2 &,int nbegin=0);
475 Real8 Length();
476 Int4 NewPoints(Vertex *,Int4 & nbv,Int4 nbvx);
477 };
478 /*}}}1*/
479 /*CLASS: GeometricalSubDomain{{{1*/
480 class GeometricalSubDomain {
481 public:
482 GeometricalEdge *edge;
483 int sens; // -1 or 1
484 Int4 ref;
485 inline void Set(const GeometricalSubDomain &,const Geometry & ,const Geometry &);
486
487 };
488 /*}}}1*/
489 /*CLASS: SubDomain{{{1*/
490 class SubDomain {
491 public:
492 Triangle * head;
493 Int4 ref;
494 int sens; // -1 or 1
495 Edge * edge; // to geometrical
496 inline void Set(const Triangles &,Int4,Triangles &);
497
498 };
499 /*}}}1*/
500 /*CLASS: VertexOnGeom{{{1*/
501 class VertexOnGeom{
502 public:
503
504 Vertex * mv;
505 Real8 abscisse;
506 union{
507 GeometricalVertex * gv; // if abscisse <0;
508 GeometricalEdge * ge; // if abscisse in [0..1]
509 };
510 inline void Set(const VertexOnGeom&,const Triangles &,Triangles &);
511 int OnGeomVertex()const {return this? abscisse <0 :0;}
512 int OnGeomEdge() const {return this? abscisse >=0 :0;}
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 operator Vertex * () const {return mv;}
517 operator GeometricalVertex * () const {return gv;}
518 operator GeometricalEdge * () const {return ge;}
519 operator const Real8 & () const {return abscisse;}
520 int IsRequiredVertex(){ return this? (( abscisse<0 ? (gv?gv->Required():0):0 )) : 0;}
521 void SetOn(){mv->onGeometry=this;mv->vint=IsVertexOnGeom;}
522 inline void Set(const Triangles &,Int4,Triangles &);
523
524 };
525 /*}}}1*/
526 /*CLASS: VertexOnVertex{{{1*/
527 class VertexOnVertex {
528 public:
529 Vertex * v, *bv;
530 VertexOnVertex(Vertex * w,Vertex *bw) :v(w),bv(bw){}
531 VertexOnVertex() {};
532 inline void Set(const Triangles &,Int4,Triangles &);
533 void SetOnBTh(){v->onBackgroundVertex=bv;v->vint=IsVertexOnVertex;}
534 };
535 /*}}}1*/
536 /*CLASS: VertexOnEdge{{{1*/
537 class VertexOnEdge {
538 public:
539 Vertex * v;
540 Edge * be;
541 Real8 abcisse;
542 VertexOnEdge( Vertex * w, Edge *bw,Real8 s) :v(w),be(bw),abcisse(s) {}
543 VertexOnEdge(){}
544 inline void Set(const Triangles &,Int4,Triangles &);
545 void SetOnBTh(){v->onBackgroundEdge=this;v->vint=IsVertexOnEdge;}
546 Vertex & operator[](int i) const { return (*be)[i];}
547 operator Real8 () const { return abcisse;}
548 operator Vertex * () const { return v;}
549 operator Edge * () const { return be;}
550 };
551 /*}}}1*/
552 /*CLASS: CrackedEdge{{{1*/
553 class CrackedEdge {
554 friend class Triangles;
555
556 class CrackedTriangle {
557 friend class Triangles;
558 friend class CrackedEdge;
559 Triangle * t; // edge of triangle t
560 int i; // edge number of in triangle
561 Edge *edge; // the 2 edge
562 Vertex *New[2]; // new vertex number
563 CrackedTriangle() : t(0),i(0),edge(0) { New[0]=New[1]=0;}
564 CrackedTriangle(Edge * a) : t(0),i(0),edge(a) { New[0]=New[1]=0;}
565 void Crack(){
566 Triangle & T(*t);
567 int i0=VerticesOfTriangularEdge[i][0];
568 int i1=VerticesOfTriangularEdge[i][0];
569 if (!New[0] && !New[1]){
570 throw ErrorException(__FUNCT__,exprintf("!New[0] && !New[1]"));
571 }
572 T(i0) = New[0];
573 T(i1) = New[1];}
574 void UnCrack(){
575 Triangle & T(*t);
576 int i0=VerticesOfTriangularEdge[i][0];
577 int i1=VerticesOfTriangularEdge[i][0];
578 if (!New[0] && !New[1]){
579 throw ErrorException(__FUNCT__,exprintf("!New[0] && !New[1]"));
580 }
581 T(i0) = TheVertex(T(i0));
582 T(i1) = TheVertex(T(i1));}
583 void Set() {
584 TriangleAdjacent ta ( FindTriangleAdjacent(*edge));
585 t = ta;
586 i = ta;
587
588 New[0] = ta.EdgeVertex(0);
589 New[1] = ta.EdgeVertex(1);
590 // warning the ref
591 }
592 };
593
594 public:
595 CrackedTriangle a,b;
596 CrackedEdge() :a(),b() {}
597 CrackedEdge(Edge * start, Int4 i,Int4 j) : a(start+i),b(start+j) {};
598 CrackedEdge(Edge * e0, Edge * e1 ) : a(e0),b(e1) {};
599
600 void Crack() { a.Crack(); b.Crack();}
601 void UnCrack() { a.UnCrack(); b.UnCrack();}
602 void Set() { a.Set(), b.Set();}
603 };
604 /*}}}1*/
605 /*CLASS: Triangles{{{1*/
606 class Triangles {
607 public:
608
609 int static counter; // to kown the number of mesh in memory
610 Geometry & Gh; // Geometry
611 Triangles & BTh; // Background Mesh Bth==*this =>no background
612
613 Int4 NbRef; // counter of ref on the this class if 0 we can delete
614 Int4 nbvx,nbtx; // nombre max de sommets , de triangles
615
616 Int4 nt,nbv,nbt,nbiv,nbe; // nb of legal triangles, nb of vertex, of triangles,
617 // of initial vertices, of edges with reference,
618 Int4 NbOfQuad; // nb of quadrangle
619
620 Int4 NbSubDomains; //
621 Int4 NbOutT; // Nb of oudeside triangle
622 Int4 NbOfTriangleSearchFind;
623 Int4 NbOfSwapTriangle;
624 Vertex * vertices; // data of vertices des sommets
625
626 Int4 NbVerticesOnGeomVertex;
627 VertexOnGeom * VerticesOnGeomVertex;
628
629 Int4 NbVerticesOnGeomEdge;
630 VertexOnGeom * VerticesOnGeomEdge;
631
632 Int4 NbVertexOnBThVertex;
633 VertexOnVertex *VertexOnBThVertex;
634
635 Int4 NbVertexOnBThEdge;
636 VertexOnEdge *VertexOnBThEdge;
637
638 Int4 NbCrackedVertices;
639
640 Int4 NbCrackedEdges;
641 CrackedEdge *CrackedEdges;
642
643 R2 pmin,pmax; // extrema
644 Real8 coefIcoor; // coef to integer Icoor1;
645
646 Triangle * triangles;
647 Edge * edges;
648
649 QuadTree *quadtree;
650 Vertex ** ordre;
651 SubDomain * subdomains;
652 ListofIntersectionTriangles lIntTria;
653 // end of variable
654
655 Triangles(Int4 i);//:BTh(*this),Gh(*new Geometry()){PreInit(i);}
656
657 ~Triangles();
658 Triangles(const char * ,Real8=-1) ;
659 Triangles(BamgMesh* bamgmesh,BamgOpts* bamgopts);
660
661 Triangles(Int4 nbvx,Triangles & BT,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) {
662 try {GeomToTriangles1(nbvx,keepBackVertices);}
663 catch(...) { this->~Triangles(); throw; }
664 }
665
666 //Genetare mesh from geometry
667 Triangles(Int4 nbvx,Geometry & G) :Gh(G),BTh(*this){
668 try { GeomToTriangles0(nbvx);}
669 catch(...) { this->~Triangles(); throw; }
670 }
671
672 Triangles(Triangles &,Geometry * pGh=0,Triangles* pBTh=0,Int4 nbvxx=0 ); // COPY OPERATEUR
673 Triangles(const Triangles &,const int *flag,const int *bb,BamgOpts* bamgopts); // truncature
674
675 void SetIntCoor(const char * from =0);
676
677 Real8 MinimalHmin() {return 2.0/coefIcoor;}
678 Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);}
679 const Vertex & operator[] (Int4 i) const { return vertices[i];};
680 Vertex & operator[](Int4 i) { return vertices[i];};
681 const Triangle & operator() (Int4 i) const { return triangles[i];};
682 Triangle & operator()(Int4 i) { return triangles[i];};
683 I2 toI2(const R2 & P) const {
684 return I2( (Icoor1) (coefIcoor*(P.x-pmin.x))
685 ,(Icoor1) (coefIcoor*(P.y-pmin.y)) );}
686 R2 toR2(const I2 & P) const {
687 return R2( (double) P.x/coefIcoor+pmin.x, (double) P.y/coefIcoor+pmin.y);
688 }
689 void Add( Vertex & s,Triangle * t,Icoor2 * =0) ;
690 void Insert();
691 void ForceBoundary();
692 void Heap();
693 void FindSubDomain(int OutSide=0);
694 Int4 ConsRefTriangle(Int4 *) const;
695 void ShowHistogram() const;
696 void ShowRegulaty() const; // Add FH avril 2007
697
698 void ReMakeTriangleContainingTheVertex();
699 void UnMarkUnSwapTriangle();
700 void SmoothMetric(Real8 raisonmax) ;
701 void BoundAnisotropy(Real8 anisomax,double hminaniso= 1e-100) ;
702 void MaxSubDivision(Real8 maxsubdiv);
703 Edge** MakeGeometricalEdgeToEdge();
704 void SetVertexFieldOn();
705 void SetVertexFieldOnBTh();
706 Int4 SplitInternalEdgeWithBorderVertices();
707 void MakeQuadrangles(double costheta);
708 int SplitElement(int choice);
709 void MakeQuadTree();
710 void NewPoints( Triangles &,int KeepVertices =1 );
711 Int4 InsertNewPoints(Int4 nbvold,Int4 & NbTSwap) ;
712 void NewPoints(int KeepVertices=1){ NewPoints(*this,KeepVertices);}
713 void ReNumberingTheTriangleBySubDomain(bool justcompress=false);
714 void ReNumberingVertex(Int4 * renu);
715 void SmoothingVertex(int =3,Real8=0.3);
716 Metric MetricAt (const R2 &) const;
717 GeometricalEdge * ProjectOnCurve( Edge & AB, Vertex & A, Vertex & B,Real8 theta,
718 Vertex & R,VertexOnEdge & BR,VertexOnGeom & GR);
719
720 Int4 Number(const Triangle & t) const { return &t - triangles;}
721 Int4 Number(const Triangle * t) const { return t - triangles;}
722 Int4 Number(const Vertex & t) const { return &t - vertices;}
723 Int4 Number(const Vertex * t) const { return t - vertices;}
724 Int4 Number(const Edge & t) const { return &t - edges;}
725 Int4 Number(const Edge * t) const { return t - edges;}
726 Int4 Number2(const Triangle * t) const {
727 // if(t>= triangles && t < triangles + nbt )
728 return t - triangles;
729 // else return t - OutSidesTriangles;
730 }
731
732 Vertex * NearestVertex(Icoor1 i,Icoor1 j) ;
733 Triangle * FindTriangleContaining(const I2 & ,Icoor2 [3],Triangle *tstart=0) const;
734
735 void ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts);
736 void WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts);
737
738 void ReadMetric(BamgOpts* bamgopts,const Real8 hmin,const Real8 hmax,const Real8 coef);
739 void WriteMetric(BamgOpts* bamgopts);
740 void IntersectConsMetric(BamgOpts* bamgopts);
741 void BuildMetric0(BamgOpts* bamgopts);
742 void BuildMetric1(BamgOpts* bamgopts);
743 void IntersectGeomMetric(BamgOpts* bamgopts);
744
745 int isCracked() const {return NbCrackedVertices != 0;}
746 int Crack();
747 int UnCrack();
748
749 void BuildGeometryFromMesh(BamgOpts* bamgopts);
750 void FillHoleInMesh() ;
751 int CrackMesh();
752 private:
753 void GeomToTriangles1(Int4 nbvx,int KeepVertices=1);// the real constructor mesh adaption
754 void GeomToTriangles0(Int4 nbvx); // the real constructor mesh generator
755 void PreInit(Int4);
756
757 };
758 /*}}}1*/
759 /*CLASS: Geometry{{{1*/
760 class Geometry {
761 public:
762 Int4 NbRef; // counter of ref on the this class if 0 we can delete
763
764 Int4 nbvx,nbtx; // nombre max de sommets , de Geometry
765 Int4 nbv,nbt,nbiv,nbe; // nombre de sommets, de Geometry, de sommets initiaux,
766 Int4 NbSubDomains; //
767 Int4 NbEquiEdges;
768 Int4 NbCrackedEdges;
769 // Int4 nbtf;// de triangle frontiere
770 Int4 NbOfCurves;
771 int empty(){return (nbv ==0) && (nbt==0) && (nbe==0) && (NbSubDomains==0); }
772 GeometricalVertex* vertices;
773 Triangle * triangles;
774 GeometricalEdge * edges;
775 QuadTree *quadtree;
776 GeometricalSubDomain *subdomains;
777 Curve *curves;
778 ~Geometry();
779 Geometry(const Geometry & Gh); //Copy Operator
780 Geometry(int nbg,const Geometry ** ag); // intersection operator
781
782 R2 pmin,pmax; // extrema
783 Real8 coefIcoor; // coef to integer Icoor1;
784 Real8 MaximalAngleOfCorner;
785
786 void Echo();
787
788 I2 toI2(const R2 & P) const {
789 return I2( (Icoor1) (coefIcoor*(P.x-pmin.x))
790 ,(Icoor1) (coefIcoor*(P.y-pmin.y)) );}
791
792 Real8 MinimalHmin() {return 2.0/coefIcoor;}
793 Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);}
794 void ReadGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts);
795 void EmptyGeometry();
796 Geometry() {EmptyGeometry();}// empty Geometry
797 void AfterRead();
798 Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts) {EmptyGeometry();ReadGeometry(bamggeom,bamgopts);AfterRead();}
799
800 const GeometricalVertex & operator[] (Int4 i) const { return vertices[i];};
801 GeometricalVertex & operator[](Int4 i) { return vertices[i];};
802 const GeometricalEdge & operator() (Int4 i) const { return edges[i];};
803 GeometricalEdge & operator()(Int4 i) { return edges[i];};
804 Int4 Number(const GeometricalVertex & t) const { return &t - vertices;}
805 Int4 Number(const GeometricalVertex * t) const { return t - vertices;}
806 Int4 Number(const GeometricalEdge & t) const { return &t - edges;}
807 Int4 Number(const GeometricalEdge * t) const { return t - edges;}
808 Int4 Number(const Curve * c) const { return c - curves;}
809
810 void UnMarkEdges() {
811 for (Int4 i=0;i<nbe;i++) edges[i].SetUnMark();}
812
813 GeometricalEdge * ProjectOnCurve(const Edge & ,Real8,Vertex &,VertexOnGeom &) const ;
814 GeometricalEdge * Contening(const R2 P, GeometricalEdge * start) const;
815 void WriteGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts);
816 };
817 /*}}}1*/
818
819 // to sort in descending order
820 template<class T> inline void HeapSort(T *c,long n){
821 long l,j,r,i;
822 T crit;
823 c--; // on decale de 1 pour que le tableau commence a 1
824 if( n <= 1) return;
825 l = n/2 + 1;
826 r = n;
827 while (1) { // label 2
828 if(l <= 1 ) { // label 20
829 crit = c[r];
830 c[r--] = c[1];
831 if ( r == 1 ) { c[1]=crit; return;}
832 } else crit = c[--l];
833 j=l;
834 while (1) {// label 4
835 i=j;
836 j=2*j;
837 if (j>r) {c[i]=crit;break;} // L8 -> G2
838 if ((j<r) && (c[j] < c[j+1])) j++; // L5
839 if (crit < c[j]) c[i]=c[j]; // L6+1 G4
840 else {c[i]=crit;break;} //L8 -> G2
841 }
842 }
843 }
844
845 //From Metric.cpp
846 inline Real8 det3x3(Real8 A[3] ,Real8 B[3],Real8 C[3]){
847 return A[0]*( B[1]*C[2]-B[2]*C[1])
848 - A[1]*( B[0]*C[2]-B[2]*C[0])
849 + A[2]*( B[0]*C[1]-B[1]*C[0]);
850 }
851
852 inline Triangles::Triangles(Int4 i) :Gh(*new Geometry()),BTh(*this){PreInit(i);}
853
854 TriangleAdjacent CloseBoundaryEdge(I2 ,Triangle *, double &,double &) ;
855 TriangleAdjacent CloseBoundaryEdgeV2(I2 A,Triangle *t, double &a,double &b);
856
857 Int4 FindTriangle(Triangles &Th, Real8 x, Real8 y, double* a,int & inside);
858
859 inline Triangle* Triangle::Quadrangle(Vertex * & v0,Vertex * & v1,Vertex * & v2,Vertex * & v3) const
860 {
861 // return the other triangle of the quad if a quad or 0 if not a quat
862 Triangle * t =0;
863 if (link) {
864 int a=-1;
865 if (TriaAdjSharedEdge[0] & 16 ) a=0;
866 if (TriaAdjSharedEdge[1] & 16 ) a=1;
867 if (TriaAdjSharedEdge[2] & 16 ) a=2;
868 if (a>=0) {
869 t = TriaAdjTriangles[a];
870 // if (t-this<0) return 0;
871 v2 = TriaVertices[VerticesOfTriangularEdge[a][0]];
872 v0 = TriaVertices[VerticesOfTriangularEdge[a][1]];
873 v1 = TriaVertices[OppositeEdge[a]];
874 v3 = t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]];
875 }
876 }
877 return t;
878 }
879
880 inline double Triangle::QualityQuad(int a,int option) const
881 { // first do the logique part
882 double q;
883 if (!link || TriaAdjSharedEdge[a] &4)
884 q= -1;
885 else {
886 Triangle * t = TriaAdjTriangles[a];
887 if (t-this<0) q= -1;// because we do 2 times
888 else if (!t->link ) q= -1;
889 else if (TriaAdjSharedEdge[0] & 16 || TriaAdjSharedEdge[1] & 16 || TriaAdjSharedEdge[2] & 16 || t->TriaAdjSharedEdge[0] & 16 || t->TriaAdjSharedEdge[1] & 16 || t->TriaAdjSharedEdge[2] & 16 )
890 q= -1;
891 else if(option)
892 {
893 const Vertex & v2 = *TriaVertices[VerticesOfTriangularEdge[a][0]];
894 const Vertex & v0 = *TriaVertices[VerticesOfTriangularEdge[a][1]];
895 const Vertex & v1 = *TriaVertices[OppositeEdge[a]];
896 const Vertex & v3 = * t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]];
897 q = QuadQuality(v0,v1,v2,v3); // do the float part
898 }
899 else q= 1;
900 }
901 return q;
902 }
903
904 inline void Vertex::Set(const Vertex & rec,const Triangles & ,Triangles & )
905 {
906 *this = rec;
907 }
908 inline void GeometricalVertex::Set(const GeometricalVertex & rec,const Geometry & ,const Geometry & )
909 {
910 *this = rec;
911 }
912 inline void Edge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
913 {
914 *this = Th.edges[i];
915 v[0] = ThNew.vertices + Th.Number(v[0]);
916 v[1] = ThNew.vertices + Th.Number(v[1]);
917 if (onGeometry)
918 onGeometry = ThNew.Gh.edges+Th.Gh.Number(onGeometry);
919 if (adj[0]) adj[0] = ThNew.edges + Th.Number(adj[0]);
920 if (adj[1]) adj[1] = ThNew.edges + Th.Number(adj[1]);
921
922 }
923 inline void GeometricalEdge::Set(const GeometricalEdge & rec,const Geometry & Gh ,Geometry & GhNew)
924 {
925 *this = rec;
926 v[0] = GhNew.vertices + Gh.Number(v[0]);
927 v[1] = GhNew.vertices + Gh.Number(v[1]);
928 if (Adj[0]) Adj[0] = GhNew.edges + Gh.Number(Adj[0]);
929 if (Adj[1]) Adj[1] = GhNew.edges + Gh.Number(Adj[1]);
930 }
931
932 inline void Curve::Set(const Curve & rec,const Geometry & Gh ,Geometry & GhNew)
933 {
934 *this = rec;
935 be = GhNew.edges + Gh.Number(be);
936 ee = GhNew.edges + Gh.Number(ee);
937 if(next) next= GhNew.curves + Gh.Number(next);
938 }
939
940 inline void Triangle::Set(const Triangle & rec,const Triangles & Th ,Triangles & ThNew)
941 {
942 *this = rec;
943 if ( TriaVertices[0] ) TriaVertices[0] = ThNew.vertices + Th.Number(TriaVertices[0]);
944 if ( TriaVertices[1] ) TriaVertices[1] = ThNew.vertices + Th.Number(TriaVertices[1]);
945 if ( TriaVertices[2] ) TriaVertices[2] = ThNew.vertices + Th.Number(TriaVertices[2]);
946 if(TriaAdjTriangles[0]) TriaAdjTriangles[0] = ThNew.triangles + Th.Number(TriaAdjTriangles[0]);
947 if(TriaAdjTriangles[1]) TriaAdjTriangles[1] = ThNew.triangles + Th.Number(TriaAdjTriangles[1]);
948 if(TriaAdjTriangles[2]) TriaAdjTriangles[2] = ThNew.triangles + Th.Number(TriaAdjTriangles[2]);
949 if (link >= Th.triangles && link < Th.triangles + Th.nbt)
950 link = ThNew.triangles + Th.Number(link);
951 }
952 inline void VertexOnVertex::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
953 {
954 *this = Th.VertexOnBThVertex[i];
955 v = ThNew.vertices + Th.Number(v);
956
957 }
958 inline void SubDomain::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
959 {
960 *this = Th.subdomains[i];
961 if ( head-Th.triangles<0 || head-Th.triangles>=Th.nbt){
962 throw ErrorException(__FUNCT__,exprintf("head-Th.triangles<0 || head-Th.triangles>=Th.nbt"));
963 }
964 head = ThNew.triangles + Th.Number(head) ;
965 if (edge-Th.edges<0 || edge-Th.edges>=Th.nbe);{
966 throw ErrorException(__FUNCT__,exprintf("edge-Th.edges<0 || edge-Th.edges>=Th.nbe"));
967 }
968 edge = ThNew.edges+ Th.Number(edge);
969 }
970 inline void GeometricalSubDomain::Set(const GeometricalSubDomain & rec,const Geometry & Gh ,const Geometry & GhNew)
971 {
972 *this = rec;
973 edge = Gh.Number(edge) + GhNew.edges;
974 }
975 inline void VertexOnEdge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
976 {
977 *this = Th.VertexOnBThEdge[i];
978 v = ThNew.vertices + Th.Number(v);
979 }
980
981 inline void VertexOnGeom::Set(const VertexOnGeom & rec,const Triangles & Th ,Triangles & ThNew)
982 {
983 *this = rec;
984 mv = ThNew.vertices + Th.Number(mv);
985 if (gv)
986 if (abscisse < 0 )
987 gv = ThNew.Gh.vertices + Th.Gh.Number(gv);
988 else
989 ge = ThNew.Gh.edges + Th.Gh.Number(ge);
990
991 }
992 inline Real8 Edge::MetricLength() const
993 {
994 return LengthInterpole(v[0]->m,v[1]->m,v[1]->r - v[0]->r) ;
995 }
996
997 inline void Triangles::ReMakeTriangleContainingTheVertex()
998 {
999 register Int4 i;
1000 for ( i=0;i<nbv;i++)
1001 {
1002 vertices[i].vint = 0;
1003 vertices[i].t=0;
1004 }
1005 for ( i=0;i<nbt;i++)
1006 triangles[i].SetTriangleContainingTheVertex();
1007 }
1008
1009 inline void Triangles::UnMarkUnSwapTriangle()
1010 {
1011 register Int4 i;
1012 for ( i=0;i<nbt;i++)
1013 for(int j=0;j<3;j++)
1014 triangles[i].SetUnMarkUnSwap(j);
1015 }
1016
1017 inline void Triangles::SetVertexFieldOn(){
1018 for (Int4 i=0;i<nbv;i++)
1019 vertices[i].onGeometry=0;
1020 for (Int4 j=0;j<NbVerticesOnGeomVertex;j++ )
1021 VerticesOnGeomVertex[j].SetOn();
1022 for (Int4 k=0;k<NbVerticesOnGeomEdge;k++ )
1023 VerticesOnGeomEdge[k].SetOn();
1024 }
1025 inline void Triangles::SetVertexFieldOnBTh(){
1026 for (Int4 i=0;i<nbv;i++)
1027 vertices[i].onGeometry=0;
1028 for (Int4 j=0;j<NbVertexOnBThVertex;j++ )
1029 VertexOnBThVertex[j].SetOnBTh();
1030 for (Int4 k=0;k<NbVertexOnBThEdge;k++ )
1031 VertexOnBThEdge[k].SetOnBTh();
1032
1033 }
1034
1035 inline void TriangleAdjacent::SetAdj2(const TriangleAdjacent & ta, int l ){
1036 //set Adjacent Triangle of a triangle
1037 if(t) {
1038 t->TriaAdjTriangles[a]=ta.t;
1039 t->TriaAdjSharedEdge[a]=ta.a|l;
1040 }
1041 if(ta.t) {
1042 ta.t->TriaAdjTriangles[ta.a] = t ;
1043 ta.t->TriaAdjSharedEdge[ta.a] = a| l ;
1044 }
1045 }
1046
1047
1048 inline int TriangleAdjacent::Locked() const
1049 { return t->TriaAdjSharedEdge[a] &4;}
1050 inline int TriangleAdjacent::Cracked() const
1051 { return t->TriaAdjSharedEdge[a] &32;}
1052 inline int TriangleAdjacent::GetAllFlag_UnSwap() const
1053 { return t->TriaAdjSharedEdge[a] & 1012;} // take all flag except MarkUnSwap
1054
1055 inline int TriangleAdjacent::MarkUnSwap() const
1056 { return t->TriaAdjSharedEdge[a] &8;}
1057
1058 inline void TriangleAdjacent::SetLock(){ t->SetLocked(a);}
1059
1060 inline void TriangleAdjacent::SetCracked() { t->SetCracked(a);}
1061
1062 inline TriangleAdjacent TriangleAdjacent::Adj() const
1063 { return t->Adj(a);}
1064
1065 inline Vertex * TriangleAdjacent::EdgeVertex(const int & i) const
1066 {return t->TriaVertices[VerticesOfTriangularEdge[a][i]]; }
1067 inline Vertex * TriangleAdjacent::OppositeVertex() const
1068 {return t->TriaVertices[bamg::OppositeVertex[a]]; }
1069 inline Icoor2 & TriangleAdjacent::det() const
1070 { return t->det;}
1071 inline TriangleAdjacent Adj(const TriangleAdjacent & a)
1072 { return a.Adj();}
1073
1074 inline TriangleAdjacent Next(const TriangleAdjacent & ta)
1075 { return TriangleAdjacent(ta.t,NextEdge[ta.a]);}
1076
1077 inline TriangleAdjacent Previous(const TriangleAdjacent & ta)
1078 { return TriangleAdjacent(ta.t,PreviousEdge[ta.a]);}
1079
1080 inline void Adj(GeometricalEdge * & on,int &i)
1081 {int j=i;i=on->DirAdj[i];on=on->Adj[j];}
1082
1083 inline Real4 qualite(const Vertex &va,const Vertex &vb,const Vertex &vc)
1084 {
1085 Real4 ret;
1086 I2 ia=va,ib=vb,ic=vc;
1087 I2 ab=ib-ia,bc=ic-ib,ac=ic-ia;
1088 Icoor2 deta=Det(ab,ac);
1089 if (deta <=0) ret = -1;
1090 else {
1091 Real8 a = sqrt((Real8) (ac,ac)),
1092 b = sqrt((Real8) (bc,bc)),
1093 c = sqrt((Real8) (ab,ab)),
1094 p = a+b+c;
1095 Real8 h= Max(Max(a,b),c),ro=deta/p;
1096 ret = ro/h;}
1097 return ret;
1098 }
1099
1100
1101 inline Triangle::Triangle(Triangles *Th,Int4 i,Int4 j,Int4 k) {
1102 Vertex *v=Th->vertices;
1103 Int4 nbv = Th->nbv;
1104 if (i<0 || j<0 || k<0){
1105 throw ErrorException(__FUNCT__,exprintf("i<0 || j<0 || k<0"));
1106 }
1107 if (i>=nbv || j>=nbv || k>=nbv){
1108 throw ErrorException(__FUNCT__,exprintf("i>=nbv || j>=nbv || k>=nbv"));
1109 }
1110 TriaVertices[0]=v+i;
1111 TriaVertices[1]=v+j;
1112 TriaVertices[2]=v+k;
1113 TriaAdjTriangles[0]=TriaAdjTriangles[1]=TriaAdjTriangles[2]=0;
1114 TriaAdjSharedEdge[0]=TriaAdjSharedEdge[1]=TriaAdjSharedEdge[2]=0;
1115 det=0;
1116 }
1117
1118 inline Triangle::Triangle(Vertex *v0,Vertex *v1,Vertex *v2){
1119 TriaVertices[0]=v0;
1120 TriaVertices[1]=v1;
1121 TriaVertices[2]=v2;
1122 TriaAdjTriangles[0]=TriaAdjTriangles[1]=TriaAdjTriangles[2]=0;
1123 TriaAdjSharedEdge[0]=TriaAdjSharedEdge[1]=TriaAdjSharedEdge[2]=0;
1124 if (v0) det=0;
1125 else {
1126 det=-1;
1127 link=NULL;};
1128 }
1129
1130 inline Real4 Triangle::qualite()
1131 {
1132 return det < 0 ? -1 : bamg::qualite(*TriaVertices[0],*TriaVertices[1],*TriaVertices[2]);
1133 }
1134
1135 Int4 inline Vertex::Optim(int i,int koption)
1136 {
1137 Int4 ret=0;
1138 if ( t && (vint >= 0 ) && (vint <3) )
1139 {
1140 ret = t->Optim(vint,koption);
1141 if(!i)
1142 {
1143 t =0; // for no future optime
1144 vint= 0; }
1145 }
1146 return ret;
1147 }
1148
1149 Icoor2 inline det(const Vertex & a,const Vertex & b,const Vertex & c)
1150 {
1151 register Icoor2 bax = b.i.x - a.i.x ,bay = b.i.y - a.i.y;
1152 register Icoor2 cax = c.i.x - a.i.x ,cay = c.i.y - a.i.y;
1153 return bax*cay - bay*cax;}
1154
1155
1156 void swap(Triangle *t1,Int1 a1,
1157 Triangle *t2,Int1 a2,
1158 Vertex *s1,Vertex *s2,Icoor2 det1,Icoor2 det2);
1159
1160
1161
1162 int inline TriangleAdjacent::swap()
1163 { return t->swap(a);}
1164
1165
1166
1167 int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb ,
1168 TriangleAdjacent & tt1,Icoor2 & dets1,
1169 Icoor2 & detsa,Icoor2 & detsb, int & nbswap);
1170
1171 int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) ;
1172
1173 // inline bofbof FH
1174 inline TriangleAdjacent FindTriangleAdjacent(Edge &E)
1175 {
1176 Vertex * a = E.v[0];
1177 Vertex * b = E.v[1];
1178
1179 Triangle * t = a->t;
1180 int i = a->vint;
1181 TriangleAdjacent ta(t,EdgesVertexTriangle[i][0]); // Previous edge
1182 if (!t || i<0 || i>=3){
1183 throw ErrorException(__FUNCT__,exprintf("!t || i<0 !! i>=3"));
1184 }
1185 if ( a!=(*t)(i)){
1186 throw ErrorException(__FUNCT__,exprintf("a!=(*t)(i)"));
1187 }
1188 int k=0;
1189 do { // turn around vertex in direct sens (trigo)
1190 k++;
1191 if (k>=20000){
1192 throw ErrorException(__FUNCT__,exprintf("k>=20000"));
1193 }
1194 // in no crack => ta.EdgeVertex(1) == a otherwise ???
1195 if (ta.EdgeVertex(1) == a && ta.EdgeVertex(0) == b) return ta; // find
1196 ta = ta.Adj();
1197 if (ta.EdgeVertex(0) == a && ta.EdgeVertex(1) == b) return ta; // find
1198 --ta;
1199 } while (t != (Triangle *)ta);
1200 throw ErrorException(__FUNCT__,exprintf("FindTriangleAdjacent: triangle not found"));
1201 return TriangleAdjacent(0,0);//for compiler
1202 }
1203
1204 inline Vertex * TheVertex(Vertex * a) // give a unique vertex with smallest number
1205 { // in case on crack in mesh
1206 Vertex * r(a), *rr;
1207 Triangle * t = a->t;
1208 int i = a->vint;
1209 TriangleAdjacent ta(t,EdgesVertexTriangle[i][0]); // Previous edge
1210 if (!t || i<0 || i>=3){
1211 throw ErrorException(__FUNCT__,exprintf("!t || i<0 !! i>=3"));
1212 }
1213 if ( a!=(*t)(i)){
1214 throw ErrorException(__FUNCT__,exprintf("a!=(*t)(i)"));
1215 }
1216 int k=0;
1217 do { // turn around vertex in direct sens (trigo)
1218 k++;
1219 if (k>=20000){
1220 throw ErrorException(__FUNCT__,exprintf("k>=20000"));
1221 }
1222 // in no crack => ta.EdgeVertex(1) == a
1223 if ((rr=ta.EdgeVertex(0)) < r) r = rr;
1224 ta = ta.Adj();
1225 if ((rr=ta.EdgeVertex(1)) < r) r =rr;
1226 --ta;
1227 } while (t != (Triangle*) ta);
1228 return r;
1229 }
1230}
1231#endif
Note: See TracBrowser for help on using the repository browser.