Context Navigation

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

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

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

Download in other formats: