Context Navigation

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

Last change on this file since 2984 was 2984, checked in by Mathieu Morlighem, 15 years ago
minor commenting and formatting
File size: 41.6 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: