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

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

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