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

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

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