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

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

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