Context Navigation

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

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

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

Download in other formats: