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

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

passed MaxAngle into options and renamed BamgArgs BamgOpts

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