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

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

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