Context Navigation

Mesh.cpp@ 22085

Last change on this file since 22085 was 22085, checked in by Eric.Larour, 8 years ago
CHG: new bamg version.
File size: 148.7 KB

Rev	Line
[3913]	1	#include <cstdio>
	2	#include <cstring>
	3	#include <cmath>
	4	#include <ctime>
	5
[15066]	6	#include "../shared/shared.h"
[15065]	7	#include "./bamgobjects.h"
[14949]	8	#include "./det.h"
[3913]	9
	10	namespace bamg {
	11
	12	/Constructors/Destructors/
[18064]	13	Mesh::Mesh(BamgGeom* bamggeom,BamgMesh* bamgmesh, BamgOpts* bamgopts):Gh((new Geometry())),BTh(this){ /{{{/
[3913]	14
	15	/Initialize fields/
[5120]	16	Init(0);
[3913]	17
	18	/Read Geometry if provided/
	19	if(bamggeom->Edges) {
	20	Gh.ReadGeometry(bamggeom,bamgopts);
[5397]	21	Gh.PostRead();
[3913]	22	}
	23
	24	/Read background mesh/
	25	ReadMesh(bamgmesh,bamgopts);
	26
	27	/Build Geometry if not provided/
	28	if(bamggeom->Edges==NULL) {
	29	/Recreate geometry if needed/
[15104]	30	_printf_("WARNING: mesh present but no geometry found. Reconstructing...\n");
[3913]	31	BuildGeometryFromMesh(bamgopts);
[5397]	32	Gh.PostRead();
[3913]	33	}
	34
	35	/Set integer coordinates/
	36	SetIntCoor();
	37
	38	/Fill holes and generate mesh properties/
	39	ReconstructExistingMesh();
	40	}
[12365]	41	/}}}/
[18064]	42	Mesh::Mesh(int* index,double* x,double* y,int nods,int nels):Gh((new Geometry())),BTh(this){/{{{/
[3913]	43
[5120]	44	Init(0);
[3913]	45	ReadMesh(index,x,y,nods,nels);
	46	SetIntCoor();
	47	ReconstructExistingMesh();
	48	}
[12365]	49	/}}}/
[18064]	50	Mesh::Mesh(double* x,double* y,int nods):Gh((new Geometry())),BTh(this){/{{{/
[10205]	51	Triangulate(x,y,nods);
	52	}
[12365]	53	/}}}/
[18064]	54	Mesh::Mesh(const Mesh & Tho,const int flag ,const int bb,BamgOpts* bamgopts) : Gh((new Geometry())), BTh(this) {/{{{/
[3913]	55	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Triangles)/
	56
	57	int i,k,itadj;
	58	int kt=0;
	59	int * kk = new int [Tho.nbv];
	60	long * reft = new long[Tho.nbt];
	61	long nbInT = Tho.TriangleReferenceList(reft);
	62	long * refv = new long[Tho.nbv];
	63
	64	for (i=0;i<Tho.nbv;i++)
	65	kk[i]=-1;
	66	for (i=0;i<Tho.nbv;i++)
	67	refv[i]=0;
	68	int nbNewBedge =0;
	69	// int nbOldBedge =0;
	70	for (i=0;i<Tho.nbt;i++)
	71	if( reft[i] >=0 && flag[i])
	72	{
	73	const Triangle & t = Tho.triangles[i];
	74	kt++;
[5149]	75	kk[Tho.GetId(t[0])]=1;
	76	kk[Tho.GetId(t[1])]=1;
	77	kk[Tho.GetId(t[2])]=1;
	78	itadj=Tho.GetId(t.TriangleAdj(0));
[3913]	79	if ( reft[itadj] >=0 && !flag[itadj])
	80	{ nbNewBedge++;
[5149]	81	refv[Tho.GetId(t[VerticesOfTriangularEdge[0][0]])]=bb[i];
	82	refv[Tho.GetId(t[VerticesOfTriangularEdge[0][1]])]=bb[i];
[3913]	83	}
[5149]	84	itadj=Tho.GetId(t.TriangleAdj(1));
[3913]	85	if ( reft[itadj] >=0 && !flag[itadj])
	86	{ nbNewBedge++;
[5149]	87	refv[Tho.GetId(t[VerticesOfTriangularEdge[1][0]])]=bb[i];
	88	refv[Tho.GetId(t[VerticesOfTriangularEdge[1][1]])]=bb[i];}
	89	itadj=Tho.GetId(t.TriangleAdj(2));
[3913]	90	if ( reft[itadj] >=0 && !flag[itadj])
	91	{ nbNewBedge++;
[5149]	92	refv[Tho.GetId(t[VerticesOfTriangularEdge[2][0]])]=bb[i];
	93	refv[Tho.GetId(t[VerticesOfTriangularEdge[2][1]])]=bb[i];}
[3913]	94	}
	95	k=0;
	96	for (i=0;i<Tho.nbv;i++){
	97	if (kk[i]>=0) kk[i]=k++;
	98	}
[15100]	99	_printf_(" number of vertices " << k << ", remove = " << Tho.nbv - k << "\n");
	100	_printf_(" number of triangles " << kt << ", remove = " << nbInT-kt << "\n");
	101	_printf_(" number of New boundary edge " << nbNewBedge << "\n");
[5120]	102	long imaxnbv =k;
	103	Init(imaxnbv);
[3913]	104	for (i=0;i<Tho.nbv;i++)
	105	if (kk[i]>=0)
	106	{
	107	vertices[nbv] = Tho.vertices[i];
[5148]	108	if (!vertices[nbv].GetReferenceNumber())
[3913]	109	vertices[nbv].ReferenceNumber = refv[i];
	110	nbv++;
	111	}
[5120]	112	if (imaxnbv != nbv){
[9306]	113	delete [] kk;
	114	delete [] refv;
[13036]	115	_error_("imaxnbv != nbv");
[3913]	116	}
	117	for (i=0;i<Tho.nbt;i++)
[16152]	118	if(reft[i] >=0 && flag[i]){
[3913]	119	const Triangle & t = Tho.triangles[i];
[5149]	120	int i0 = Tho.GetId(t[0]);
	121	int i1 = Tho.GetId(t[1]);
	122	int i2 = Tho.GetId(t[2]);
[16152]	123	if(i0<0 \|\| i1<0 \|\| i2<0){
[9309]	124	delete [] refv;
[13036]	125	_error_("i0<0 \|\| i1<0 \|\| i2< 0");
[3913]	126	}
[16152]	127	if(i0>=Tho.nbv \|\| i1>=Tho.nbv \|\| i2>=Tho.nbv){
	128	delete [] refv;
[13036]	129	_error_("i0>=Tho.nbv \|\| i1>=Tho.nbv \|\| i2>=Tho.nbv");
[3913]	130	}
	131	triangles[nbt] = Triangle(this,kk[i0],kk[i1],kk[i2]);
[5148]	132	triangles[nbt].color = Tho.subdomains[reft[i]].ReferenceNumber;
[3913]	133	nbt++;
	134	}
[22085]	135	if (nbt==0 && nbv==0){
	136	delete [] refv;
[13036]	137	_error_("All triangles have been removed");
[3913]	138	}
	139	delete [] kk;
	140	delete [] reft;
	141	delete [] refv;
	142	BuildGeometryFromMesh(bamgopts);
[5397]	143	Gh.PostRead();
[3913]	144	SetIntCoor();
	145	ReconstructExistingMesh();
	146
[16237]	147	/Final checks/
	148	_assert_(kt==nbt);
	149	_assert_(nbsubdomains);
	150	_assert_(subdomains[0].head && subdomains[0].head->link);
[3913]	151	}
[12365]	152	/}}}/
[18064]	153	Mesh::Mesh(Mesh & Th,Geometry * pGh,Mesh * pBth,long maxnbv_in)/{{{/
[3913]	154	: Gh((pGh?pGh:&Th.Gh)), BTh((pBth?pBth:this)) {
	155	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Triangles)/
	156	Gh.NbRef++;
[5120]	157	maxnbv_in = Max(maxnbv_in,Th.nbv);
[3913]	158	long i;
	159	// do all the allocation to be sure all the pointer existe
	160
[5120]	161	Init(maxnbv_in);// to make the allocation
[3913]	162	// copy of triangles
	163	nbv = Th.nbv;
	164	nbt = Th.nbt;
	165	nbe = Th.nbe;
[5148]	166	nbsubdomains = Th.nbsubdomains;
[5146]	167	nbtout = Th.nbtout;
[3913]	168	NbVerticesOnGeomVertex = Th.NbVerticesOnGeomVertex;
	169	if(NbVerticesOnGeomVertex)
	170	VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];
	171	NbVerticesOnGeomEdge = Th.NbVerticesOnGeomEdge;
	172	if (NbVerticesOnGeomEdge)
	173	VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge] ;
	174	if (& BTh == & Th.BTh){ // same background
	175	BTh.NbRef++;
	176	NbVertexOnBThVertex = Th.NbVertexOnBThVertex;
	177	if(NbVertexOnBThVertex)
	178	VertexOnBThVertex = new VertexOnVertex[NbVertexOnBThVertex];
	179	NbVertexOnBThEdge = Th.NbVertexOnBThEdge;
	180	if(NbVertexOnBThEdge)
	181	VertexOnBThEdge = new VertexOnEdge[NbVertexOnBThEdge];
	182	}
	183	else { // no add on background mesh
	184	BTh.NbRef++;
	185	NbVertexOnBThVertex=0;
	186	VertexOnBThVertex=0;
	187	NbVertexOnBThEdge=0;
	188	VertexOnBThEdge=0;
	189	}
	190
	191	if(nbe)
	192	edges = new Edge[nbe];
[5148]	193	if(nbsubdomains)
	194	subdomains = new SubDomain[nbsubdomains];
[3913]	195	pmin = Th.pmin;
	196	pmax = Th.pmax;
	197	coefIcoor = Th.coefIcoor;
	198	for(i=0;i<nbt;i++)
	199	triangles[i].Set(Th.triangles[i],Th,*this);
	200	for(i=0;i<nbe;i++)
	201	edges[i].Set(Th,i,*this);
	202	for(i=0;i<nbv;i++)
	203	vertices[i].Set(Th.vertices[i],Th,*this);
[5148]	204	for(i=0;i<nbsubdomains;i++)
[3913]	205	subdomains[i].Set(Th,i,*this);
	206	for (i=0;i<NbVerticesOnGeomVertex;i++)
	207	VerticesOnGeomVertex[i].Set(Th.VerticesOnGeomVertex[i],Th,*this);
	208	for (i=0;i<NbVerticesOnGeomEdge;i++)
	209	VerticesOnGeomEdge[i].Set(Th.VerticesOnGeomEdge[i],Th,*this);
	210	quadtree=0;
	211
	212	}
[12365]	213	/}}}/
[18064]	214	Mesh::Mesh(long imaxnbv,Mesh & BT,BamgOpts* bamgopts,int keepBackVertices) :Gh(BT.Gh),BTh(BT) {/{{{/
[5460]	215	this->Init(imaxnbv);
	216	TriangulateFromGeom1(bamgopts,keepBackVertices);
[5130]	217	}
[12365]	218	/}}}/
[18064]	219	Mesh::Mesh(long imaxnbv,Geometry & G,BamgOpts* bamgopts):Gh(G),BTh(this){/{{{*/
[5460]	220	Init(imaxnbv);
	221	TriangulateFromGeom0(bamgopts);
[5130]	222	}
[12365]	223	/}}}/
[18064]	224	Mesh::~Mesh() {/{{{/
[3913]	225	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Triangles)/
	226
[18453]	227	if(vertices) delete [] vertices;
	228	if(edges) delete [] edges;
	229	if(triangles) delete [] triangles;
	230	if(quadtree) delete quadtree;
	231	if(orderedvertices) delete [] orderedvertices;
	232	if(subdomains) delete [] subdomains;
	233	if(VerticesOnGeomEdge) delete [] VerticesOnGeomEdge;
	234	if(VerticesOnGeomVertex) delete [] VerticesOnGeomVertex;
	235	if(VertexOnBThVertex) delete [] VertexOnBThVertex;
	236	if(VertexOnBThEdge) delete [] VertexOnBThEdge;
[3913]	237
[20621]	238	if (Gh.NbRef>0) Gh.NbRef--;
	239	else if (Gh.NbRef==0) delete &Gh;
	240	if(&BTh != this){
[3913]	241	if (BTh.NbRef>0) BTh.NbRef--;
	242	else if (BTh.NbRef==0) delete &BTh;
	243	}
[5120]	244	Init(0); // set all to zero
[3913]	245	}
[12365]	246	/}}}/
[3913]	247
	248	/IO/
[18064]	249	void Mesh::ReadMesh(int* index,double* x,double* y,int nods,int nels){/{{{/
[3913]	250
[5143]	251	long i1,i2,i3;
[13761]	252	long i;
[3913]	253	Metric M1(1);
	254	int verbose=0;
[9202]	255	bool* nodeflags=NULL;
[3913]	256
	257	nbv=nods;
[5120]	258	maxnbv=nbv;
[3913]	259	nbt=nels;
	260
	261	//Vertices
[15104]	262	if (verbose) _printf_("Reading vertices (" << nbv << ")\n");
[12456]	263	vertices=xNew<BamgVertex>(nbv);
	264	orderedvertices=xNew<BamgVertex*>(nbv);
[3913]	265	for (i=0;i<nbv;i++){
	266	vertices[i].r.x=x[i];
	267	vertices[i].r.y=y[i];
	268	vertices[i].ReferenceNumber=1;
	269	vertices[i].m=M1;
	270	vertices[i].color=0;
	271	}
[5146]	272	maxnbt=2*maxnbv-2; // for filling The Holes and quadrilaterals
[3913]	273
	274	//Triangles
[15104]	275	if (verbose) _printf_("Reading triangles (" << nbt << ")\n");
[5146]	276	triangles =new Triangle[maxnbt]; //we cannot allocate only nbt triangles since
[12456]	277	nodeflags=xNew<bool>(nbv);
[9202]	278	for(i=0;i<nbv;i++) nodeflags[i]=false;
[3913]	279	//other triangles will be added for each edge
	280	for (i=0;i<nbt;i++){
	281	Triangle & t = triangles[i];
	282	i1=(long)index[i*3+0]-1; //for C indexing
	283	i2=(long)index[i*3+1]-1; //for C indexing
	284	i3=(long)index[i*3+2]-1; //for C indexing
	285	t=Triangle(this,i1,i2,i3);
	286	t.color=1;
[9202]	287	nodeflags[i1]=nodeflags[i2]=nodeflags[i3]=true;
[3913]	288	}
	289
	290	/Recreate geometry: /
[15104]	291	if (verbose) _printf_("Building Geometry\n");
[3913]	292	BuildGeometryFromMesh();
[15104]	293	if (verbose) _printf_("Completing geometry\n");
[5397]	294	Gh.PostRead();
[9202]	295
	296	/Check that there is no orphan/
	297	bool isorphan=false;
	298	for(i=0;i<nbv;i++){
	299	if(!nodeflags[i]){
[15104]	300	_printf_("Vertex " << i+1 << " does not belong to any element\n");
[9202]	301	isorphan=true;
	302	}
	303	}
[13036]	304	if(isorphan) _error_("Orphan found in mesh, see ids above");
[9202]	305
	306	/Clean up/
[12456]	307	xDelete<bool>(nodeflags);
[3913]	308	}
[12365]	309	/}}}/
[18064]	310	void Mesh::ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts){/{{{/
[3913]	311
[5143]	312	int verbose;
	313	double Hmin = HUGE_VAL; // the infinie value
	314	long i1,i2,i3;
	315	long i,j;
[3913]	316	Metric M1(1);
	317
	318	verbose=bamgopts->verbose;
	319
	320	nbv=bamgmesh->VerticesSize[0];
[5120]	321	maxnbv=nbv;
[3913]	322	nbt=bamgmesh->TrianglesSize[0];
	323
	324	//Vertices
	325	if(bamgmesh->Vertices){
[15104]	326	if(verbose>5) _printf_(" processing Vertices\n");
[3913]	327
[12456]	328	vertices=xNew<BamgVertex>(nbv);
	329	orderedvertices=xNew<BamgVertex*>(nbv);
[3913]	330
	331	for (i=0;i<nbv;i++){
	332	vertices[i].r.x=bamgmesh->Vertices[i*3+0];
	333	vertices[i].r.y=bamgmesh->Vertices[i*3+1];
	334	vertices[i].ReferenceNumber=(long)bamgmesh->Vertices[i*3+2];
	335	vertices[i].m=M1;
	336	vertices[i].color=0;
	337	}
[5146]	338	maxnbt=2*maxnbv-2; // for filling The Holes and quadrilaterals
[3913]	339	}
	340	else{
[13036]	341	if(verbose>5) _error_("no Vertices found in the initial mesh");
[3913]	342	}
	343
	344	//Triangles
	345	if(bamgmesh->Triangles){
[15104]	346	if(verbose>5) _printf_(" processing Triangles\n");
[5146]	347	triangles =new Triangle[maxnbt]; //we cannot allocate only nbt triangles since
[3913]	348	//other triangles will be added for each edge
	349	for (i=0;i<nbt;i++){
[5143]	350	Triangle &t=triangles[i];
[3913]	351	i1=(long)bamgmesh->Triangles[i*4+0]-1; //for C indexing
	352	i2=(long)bamgmesh->Triangles[i*4+1]-1; //for C indexing
	353	i3=(long)bamgmesh->Triangles[i*4+2]-1; //for C indexing
	354	t=Triangle(this,i1,i2,i3);
	355	t.color=(long)bamgmesh->Triangles[i*4+3];
	356	}
	357	}
	358	else{
[13036]	359	if(verbose>5) _error_("no Triangles found in the initial mesh");
[3913]	360	}
	361
[5573]	362	//VerticesOnGeomEdge
	363	if(bamgmesh->VerticesOnGeomEdge){
[15104]	364	if(verbose>5) _printf_(" processing VerticesOnGeomEdge\n");
[5573]	365	NbVerticesOnGeomEdge=bamgmesh->VerticesOnGeomEdgeSize[0];
[3913]	366	VerticesOnGeomEdge= new VertexOnGeom[NbVerticesOnGeomEdge] ;
	367	for (i=0;i<NbVerticesOnGeomEdge;i++){
	368	long i1,i2;
	369	double s;
[5573]	370	i1=(long) bamgmesh->VerticesOnGeomEdge[i*3+0]-1; //for C indexing
	371	i2=(long) bamgmesh->VerticesOnGeomEdge[i*3+1]-1; //for C indexing
	372	s =(double)bamgmesh->VerticesOnGeomEdge[i*3+2];
[3913]	373	VerticesOnGeomEdge[i]=VertexOnGeom(vertices[i1],Gh.edges[i2],s);
	374	}
	375	}
	376
[5573]	377	//VerticesOnGeomVertex
	378	if(bamgmesh->VerticesOnGeomVertexSize[0]){
[15104]	379	if(verbose>5) _printf_(" processing VerticesOnGeomVertex\n");
[5573]	380	NbVerticesOnGeomVertex=bamgmesh->VerticesOnGeomVertexSize[0];
[3913]	381	VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex] ;
	382	for (i=0;i<NbVerticesOnGeomVertex;i++){
	383	long i1,i2;
[5573]	384	i1=(long)bamgmesh->VerticesOnGeomVertex[i*2+0]-1; //for C indexing
	385	i2=(long)bamgmesh->VerticesOnGeomVertex[i*2+1]-1; //for C indexing
[3913]	386	VerticesOnGeomVertex[i]=VertexOnGeom(vertices[i1],Gh.vertices[i2]);
	387	}
	388	}
	389
	390	//Edges
	391	if (bamgmesh->Edges){
	392	int i1,i2;
	393	double* len=NULL;
	394
[15104]	395	if(verbose>5) _printf_(" processing Edges\n");
[3913]	396	nbe=bamgmesh->EdgesSize[0];
	397	edges= new Edge[nbe];
	398	//initialize length of each edge (used to provided metric)
	399	len= new double[nbv];
	400	for(i=0;i<nbv;i++) len[i]=0;
	401
	402	for (i=0;i<nbe;i++){
	403	i1=(int)bamgmesh->Edges[i*3+0]-1; //-1 for C indexing
	404	i2=(int)bamgmesh->Edges[i*3+1]-1; //-1 for C indexing
[5148]	405	edges[i].ReferenceNumber=(long)bamgmesh->Edges[i*3+2];
[3913]	406	edges[i].v[0]= vertices +i1;
	407	edges[i].v[1]= vertices +i2;
	408	edges[i].adj[0]=NULL;
	409	edges[i].adj[1]=NULL;
	410	R2 x12=vertices[i2].r-vertices[i1].r;
	411	double l12=sqrt((x12,x12));
	412
	413	//prepare metric
	414	vertices[i1].color++;
	415	vertices[i2].color++;
	416	len[i1]+=l12;
	417	len[i2]+=l12;
	418	Hmin = Min(Hmin,l12);
	419	}
	420
	421	// definition the default of the given mesh size
	422	for (i=0;i<nbv;i++){
	423	if (vertices[i].color>0)
	424	vertices[i].m=Metric(len[i]/(double)vertices[i].color);
	425	else
	426	vertices[i].m=Metric(Hmin);
	427	}
	428	delete [] len;
	429
	430	// construction of edges[].adj
	431	for (i=0;i<nbv;i++){
	432	vertices[i].color=(vertices[i].color ==2) ?-1:-2;
	433	}
	434	for (i=0;i<nbe;i++){
	435	for (j=0;j<2;j++) {
[5120]	436	BamgVertex *v=edges[i].v[j];
[3913]	437	long i0=v->color,j0;
	438	if(i0==-1){
	439	v->color=i*2+j;
	440	}
	441	else if (i0>=0) {// i and i0 edge are adjacent by the vertex v
	442	j0 = i0%2;
	443	i0 = i0/2;
[6412]	444	_assert_(v==edges[i0 ].v[j0]);
[3913]	445	edges[i ].adj[j ] =edges +i0;
	446	edges[i0].adj[j0] =edges +i ;
	447	v->color = -3;
	448	}
	449	}
	450	}
	451	}
	452
[5573]	453	//EdgeOnGeomEdge
	454	if(bamgmesh->EdgesOnGeomEdge){
[15104]	455	if(verbose>5) _printf_(" processing EdgesOnGeomEdge\n");
[3913]	456	int i1,i2,i,j;
[5573]	457	i2=bamgmesh->EdgesOnGeomEdgeSize[0];
[3913]	458	for (i1=0;i1<i2;i1++) {
[5573]	459	i=(int)bamgmesh->EdgesOnGeomEdge[i1*2+0]-1; //C indexing
	460	j=(int)bamgmesh->EdgesOnGeomEdge[i1*2+1]-1; //C indexing
[3913]	461	//Check value
	462	if(!(i>=0 && j>=0 && i<nbe && j<Gh.nbe)) {
[13036]	463	_error_("ReadMesh error: EdgesOnGeomEdge edge provided (line " << i1+1 << ": [" << i+1 << " " << j+1 << "]) is incorrect (must be positive, [0<i<nbe=" << nbe << " 0<j<Gh.nbe=" << Gh.nbe << "]");
[3913]	464	}
[5573]	465	edges[i].GeomEdgeHook=Gh.edges+j;
[3913]	466	}
	467	}
	468
	469	//SubDomain
	470	if(bamgmesh->SubDomains){
[5148]	471	long i3,head,direction;
[15104]	472	if(verbose>5) _printf_(" processing SubDomains\n");
[5148]	473	nbsubdomains=bamgmesh->SubDomainsSize[0];
	474	subdomains = new SubDomain [ nbsubdomains ];
	475	for (i=0;i<nbsubdomains;i++) {
[21759]	476	i3 =(int)bamgmesh->SubDomains[i*4+0];
	477	head=(int)bamgmesh->SubDomains[i*4+1]-1;//C indexing
	478	direction=(int)bamgmesh->SubDomains[i*4+2];
[16166]	479	if (i3!=3) _error_("Bad Subdomain definition: first number should be 3");
[13036]	480	if (head<0 \|\| head>=nbt) _error_("Bad Subdomain definition: head should in [1 " << nbt << "] (triangle number)");
[3913]	481	subdomains[i].head = triangles+head;
[16166]	482	subdomains[i].direction = direction;
	483	subdomains[i].ReferenceNumber = i3;
[3913]	484	}
	485	}
	486
	487	}
[12365]	488	/}}}/
[18064]	489	void Mesh::WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts){/{{{/
[3913]	490
	491	/Intermediary/
	492	int i,j,k,num,i1,i2;
	493	long n;
	494	int* head_1=NULL;
	495	int* next_1=NULL;
	496	int* connectivitysize_1=NULL;
	497	int connectivitymax_1=0;
	498
	499	/Get options/
	500	int verbose=bamgopts->verbose;
	501
	502	/Build reft that holds the number the subdomain number of each triangle, and the real numbering of the elements/
	503	long* reft = new long[nbt];
	504	long* numt = new long[nbt];
	505	long nbInT = TriangleReferenceList(reft);
	506	TriangleIntNumbering(numt);
	507
	508	/Chaining algorithm used to generate connectivity tables and other outputs/
	509
	510	//Memory Allocation
[12456]	511	head_1=xNew<int>(nbv);
	512	next_1=xNew<int>(3*nbt);
	513	connectivitysize_1=xNew<int>(nbv);
[3913]	514
	515	//Initialization
	516	for (i=0;i<nbv;i++) head_1[i]=-1;
	517	for (i=0;i<nbv;i++) connectivitysize_1[i]=0;
	518	k=0;
	519	//Chains generation
	520	for (i=0;i<nbt;i++) {
	521	//Do not take into account outside triangles (reft<0)
	522	if (reft[i]>=0){
	523	for (j=0;j<3;j++){
[5149]	524	int v=GetId(triangles[i][j]); //jth vertex of the ith triangle
[13036]	525	if (k>3*nbt-1 \|\| k<0) _error_("k = " << k << ", nbt = " << nbt);
[3913]	526	next_1[k]=head_1[v];
[13036]	527	if (v>nbv-1 \|\| v<0) _error_("v = " << v << ", nbv = " << nbv);
[3913]	528	head_1[v]=k++;
	529	connectivitysize_1[v]+=1;
	530	}
	531	}
	532	}
	533	//Get maximum connectivity
	534	connectivitymax_1=0;
	535	for (i=0;i<nbv;i++){
	536	if (connectivitysize_1[i]>connectivitymax_1) connectivitymax_1=connectivitysize_1[i];
	537	}
	538
	539	/OK, now build outputs/
	540
	541	/Vertices/
[15104]	542	if(verbose>5) _printf_(" writing Vertices\n");
[3913]	543	bamgmesh->VerticesSize[0]=nbv;
	544	bamgmesh->VerticesSize[1]=3;
[18455]	545	if(nbv){
[12456]	546	bamgmesh->Vertices=xNew<double>(3*nbv);
[18455]	547	bamgmesh->PreviousNumbering=xNew<double>(nbv);
[3913]	548	for (i=0;i<nbv;i++){
	549	bamgmesh->Vertices[i*3+0]=vertices[i].r.x;
	550	bamgmesh->Vertices[i*3+1]=vertices[i].r.y;
[5148]	551	bamgmesh->Vertices[i*3+2]=vertices[i].GetReferenceNumber();
[18455]	552	bamgmesh->PreviousNumbering[i]=vertices[i].PreviousNumber;
[3913]	553	}
	554	}
	555
	556	/Edges/
[15104]	557	if(verbose>5) _printf_(" writing Edges\n");
[3913]	558	bamgmesh->EdgesSize[0]=nbe;
	559	bamgmesh->EdgesSize[1]=3;
[5091]	560	int NumIssmSegments=0;
[3913]	561	if (nbe){
[12456]	562	bamgmesh->Edges=xNew<double>(3*nbe);
[3913]	563	for (i=0;i<nbe;i++){
[5149]	564	bamgmesh->Edges[i*3+0]=GetId(edges[i][0])+1; //back to M indexing
	565	bamgmesh->Edges[i*3+1]=GetId(edges[i][1])+1; //back to M indexing
[5148]	566	bamgmesh->Edges[i*3+2]=edges[i].ReferenceNumber;
[5573]	567	if(edges[i].GeomEdgeHook){
[5091]	568	NumIssmSegments++;
[3913]	569	}
	570	}
	571	}
	572
	573	/Element edges/
[15104]	574	if(verbose>5) _printf_(" writing element edges\n");
[3913]	575	SetOfEdges4* edge4=new SetOfEdges4(nbt*3,nbv);
	576	double* elemedge=NULL;
[12456]	577	elemedge=xNew<double>(3*nbt);
[12574]	578	for (i=0;i<3*nbt;i++) elemedge[i]=-2.;//will become -1
[3913]	579	k=0;
	580	for (i=0;i<nbt;i++){
	581	//Do not take into account outside triangles (reft<0)
	582	if (reft[i]>=0){
	583	for (j=0;j<3;j++) {
[5149]	584	i1=GetId(triangles[i][VerticesOfTriangularEdge[j][0]]);
	585	i2=GetId(triangles[i][VerticesOfTriangularEdge[j][1]]);
[3913]	586	n =edge4->SortAndFind(i1,i2);
	587	if (n==-1){
	588	//first time
	589	n=edge4->SortAndAdd(i1,i2);
	590	elemedge[n*2+0]=double(k);
	591	}
	592	else{
	593	//second time
	594	elemedge[n*2+1]=double(k);
	595	}
	596	}
	597	k++;
	598	}
	599	}
[5091]	600	bamgmesh->IssmEdgesSize[0]=edge4->nb();
	601	bamgmesh->IssmEdgesSize[1]=4;
[12456]	602	bamgmesh->IssmEdges=xNew<double>(4*edge4->nb());
[3913]	603	for (i=0;i<edge4->nb();i++){
	604	/Invert first two vertices if necessary/
	605	bool found=false;
	606	for (j=0;j<3;j++){
	607	if (triangles[(int)elemedge[2*i+0]](j)==vertices+edge4->i(i)){
	608	if (triangles[(int)elemedge[2*i+0]]((j+1)%3)==vertices+edge4->j(i)){
	609	//trigonometric direction
[5091]	610	bamgmesh->IssmEdges[i*4+0]=edge4->i(i)+1;// back to M indexing
	611	bamgmesh->IssmEdges[i*4+1]=edge4->j(i)+1;// back to M indexing
[3913]	612	}
	613	else{
[5091]	614	bamgmesh->IssmEdges[i*4+0]=edge4->j(i)+1;// back to M indexing
	615	bamgmesh->IssmEdges[i*4+1]=edge4->i(i)+1;// back to M indexing
[3913]	616	}
	617	found=true;
	618	break;
	619	}
	620	}
[6412]	621	_assert_(found);
[5091]	622	bamgmesh->IssmEdges[i4+2]=elemedge[2i+0]+1; // back to M indexing
	623	bamgmesh->IssmEdges[i4+3]=elemedge[2i+1]+1; // back to M indexing
[3913]	624	}
	625	//clean up
	626	delete edge4;
[12456]	627	xDelete<double>(elemedge);
[3913]	628
[5091]	629	/IssmSegments/
[15104]	630	if(verbose>5) _printf_(" writing IssmSegments\n");
[5091]	631	bamgmesh->IssmSegmentsSize[0]=NumIssmSegments;
	632	bamgmesh->IssmSegmentsSize[1]=4;
[12456]	633	bamgmesh->IssmSegments=xNew<double>(4*NumIssmSegments);
[3913]	634	num=0;
	635	for (i=0;i<nbe;i++){
[5573]	636	if(edges[i].GeomEdgeHook){
[3913]	637	//build segment
[5149]	638	int i1=GetId(edges[i][0]);
	639	int i2=GetId(edges[i][1]);
[3913]	640	bool stop=false;
	641	for(j=head_1[i1];j!=-1;j=next_1[j]){
	642	for(k=0;k<3;k++){
[5149]	643	if (GetId(triangles[(int)j/3][k])==i1){
	644	if (GetId(triangles[(int)j/3][(int)((k+1)%3)])==i2){
	645	bamgmesh->IssmSegments[num*4+0]=GetId(edges[i][0])+1; //back to M indexing
	646	bamgmesh->IssmSegments[num*4+1]=GetId(edges[i][1])+1; //back to M indexing
[5091]	647	bamgmesh->IssmSegments[num*4+2]=(int)j/3+1; //back to M indexing
[5148]	648	bamgmesh->IssmSegments[num*4+3]=edges[i].ReferenceNumber;
[3913]	649	num+=1;
	650	stop=true;
	651	break;
	652	}
[5149]	653	if (GetId(triangles[(int)j/3][(int)((k+2)%3)])==i2){
	654	bamgmesh->IssmSegments[num*4+0]=GetId(edges[i][1])+1; //back to M indexing
	655	bamgmesh->IssmSegments[num*4+1]=GetId(edges[i][0])+1; //back to M indexing
[5091]	656	bamgmesh->IssmSegments[num*4+2]=(int)j/3+1; //back to M indexing
[5148]	657	bamgmesh->IssmSegments[num*4+3]=edges[i].ReferenceNumber;
[3913]	658	num+=1;
	659	stop=true;
	660	break;
	661	}
	662	}
	663	}
	664	if(stop) break;
	665	}
	666	if (!stop){
[13036]	667	_error_("Element holding segment [" << i1+1 << " " << i2+1 << "] not found...");
[3913]	668	}
	669	}
	670	}
	671
	672	/Triangles/
[15104]	673	if(verbose>5) _printf_(" writing Triangles\n");
[21759]	674	k=nbInT;
[3913]	675	num=0;
	676	bamgmesh->TrianglesSize[0]=k;
	677	bamgmesh->TrianglesSize[1]=4;
	678	if (k){
[12456]	679	bamgmesh->Triangles=xNew<double>(4*k);
[3913]	680	for (i=0;i<nbt;i++){
	681	Triangle &t=triangles[i];
	682	//reft[i]=-1 for outside triangle
	683	if (reft[i]>=0 && !( t.Hidden(0) \|\| t.Hidden(1) \|\| t.Hidden(2) )){
[5149]	684	bamgmesh->Triangles[num*4+0]=GetId(t[0])+1; //back to M indexing
	685	bamgmesh->Triangles[num*4+1]=GetId(t[1])+1; //back to M indexing
	686	bamgmesh->Triangles[num*4+2]=GetId(t[2])+1; //back to M indexing
[5148]	687	bamgmesh->Triangles[num*4+3]=subdomains[reft[i]].ReferenceNumber;
[3913]	688	num=num+1;
	689	}
	690	}
	691	}
	692
	693	/SubDomains/
[15104]	694	if(verbose>5) _printf_(" writing SubDomains\n");
[5148]	695	bamgmesh->SubDomainsSize[0]=nbsubdomains;
[3913]	696	bamgmesh->SubDomainsSize[1]=4;
[5148]	697	if (nbsubdomains){
[12456]	698	bamgmesh->SubDomains=xNew<double>(4*nbsubdomains);
[5148]	699	for (i=0;i<nbsubdomains;i++){
[3913]	700	bamgmesh->SubDomains[i*4+0]=3;
[22085]	701	bamgmesh->SubDomains[i*4+1]=reft[GetId(subdomains[i].head)]+1;//MATLAB indexing
[3913]	702	bamgmesh->SubDomains[i*4+2]=1;
[5148]	703	bamgmesh->SubDomains[i*4+3]=subdomains[i].ReferenceNumber;
[3913]	704	}
	705	}
	706
	707	/SubDomainsFromGeom/
[15104]	708	if(verbose>5) _printf_(" writing SubDomainsFromGeom\n");
[5148]	709	bamgmesh->SubDomainsFromGeomSize[0]=Gh.nbsubdomains;
[3913]	710	bamgmesh->SubDomainsFromGeomSize[1]=4;
[5148]	711	if (Gh.nbsubdomains){
[12456]	712	bamgmesh->SubDomainsFromGeom=xNew<double>(4*Gh.nbsubdomains);
[5148]	713	for (i=0;i<Gh.nbsubdomains;i++){
[3913]	714	bamgmesh->SubDomainsFromGeom[i*4+0]=2;
[5149]	715	bamgmesh->SubDomainsFromGeom[i*4+1]=GetId(subdomains[i].edge)+1; //back to Matlab indexing
[5148]	716	bamgmesh->SubDomainsFromGeom[i*4+2]=subdomains[i].direction;
	717	bamgmesh->SubDomainsFromGeom[i*4+3]=Gh.subdomains[i].ReferenceNumber;
[3913]	718	}
	719	}
	720
	721	/VerticesOnGeomVertex/
[15104]	722	if(verbose>5) _printf_(" writing VerticesOnGeomVertex\n");
[5573]	723	bamgmesh->VerticesOnGeomVertexSize[0]=NbVerticesOnGeomVertex;
	724	bamgmesh->VerticesOnGeomVertexSize[1]=2;
[3913]	725	if (NbVerticesOnGeomVertex){
[12456]	726	bamgmesh->VerticesOnGeomVertex=xNew<double>(2*NbVerticesOnGeomVertex);
[3913]	727	for (i=0;i<NbVerticesOnGeomVertex;i++){
	728	VertexOnGeom &v=VerticesOnGeomVertex[i];
[6412]	729	_assert_(v.OnGeomVertex());
[5573]	730	bamgmesh->VerticesOnGeomVertex[i2+0]=GetId((BamgVertex)v)+1; //back to Matlab indexing
	731	bamgmesh->VerticesOnGeomVertex[i2+1]=Gh.GetId((GeomVertex)v)+1; //back to Matlab indexing
[3913]	732	}
	733	}
	734
[5573]	735	/VertexOnGeomEdge/
[15104]	736	if(verbose>5) _printf_(" writing VerticesOnGeomEdge\n");
[5573]	737	bamgmesh->VerticesOnGeomEdgeSize[0]=NbVerticesOnGeomEdge;
	738	bamgmesh->VerticesOnGeomEdgeSize[1]=3;
[3913]	739	if (NbVerticesOnGeomEdge){
[12456]	740	bamgmesh->VerticesOnGeomEdge=xNew<double>(3*NbVerticesOnGeomEdge);
[3913]	741	for (i=0;i<NbVerticesOnGeomEdge;i++){
	742	const VertexOnGeom &v=VerticesOnGeomEdge[i];
	743	if (!v.OnGeomEdge()){
[13036]	744	_error_("A vertices supposed to be OnGeomEdge is actually not");
[3913]	745	}
[5573]	746	bamgmesh->VerticesOnGeomEdge[i3+0]=GetId((BamgVertex)v)+1; //back to Matlab indexing
	747	bamgmesh->VerticesOnGeomEdge[i3+1]=Gh.GetId((const GeomEdge)v)+1; //back to Matlab indexing
	748	bamgmesh->VerticesOnGeomEdge[i*3+2]=(double)v; //absisce
[3913]	749	}
	750	}
	751
[5573]	752	/EdgesOnGeomEdge/
[15104]	753	if(verbose>5) _printf_(" writing EdgesOnGeomEdge\n");
[3913]	754	k=0;
	755	for (i=0;i<nbe;i++){
[5573]	756	if (edges[i].GeomEdgeHook) k=k+1;
[3913]	757	}
[5573]	758	bamgmesh->EdgesOnGeomEdgeSize[0]=k;
	759	bamgmesh->EdgesOnGeomEdgeSize[1]=2;
[3913]	760	if (k){
[12456]	761	bamgmesh->EdgesOnGeomEdge=xNew<double>(2*(int)k);
[3913]	762	int count=0;
	763	for (i=0;i<nbe;i++){
[5573]	764	if (edges[i].GeomEdgeHook){
	765	bamgmesh->EdgesOnGeomEdge[count*2+0]=(double)i+1; //back to Matlab indexing
	766	bamgmesh->EdgesOnGeomEdge[count*2+1]=(double)Gh.GetId(edges[i].GeomEdgeHook)+1; //back to Matlab indexing
[3913]	767	count=count+1;
	768	}
	769	}
	770	}
	771
	772	/Element Connectivity/
[15104]	773	if(verbose>5) _printf_(" writing Element connectivity\n");
[5146]	774	bamgmesh->ElementConnectivitySize[0]=nbt-nbtout;
[3913]	775	bamgmesh->ElementConnectivitySize[1]=3;
[12456]	776	bamgmesh->ElementConnectivity=xNew<double>(3*(nbt-nbtout));
[5146]	777	for (i=0;i<3*(nbt-nbtout);i++) bamgmesh->ElementConnectivity[i]=NAN;
[3913]	778	num=0;
	779	for (i=0;i<nbt;i++){
	780	if (reft[i]>=0){
	781	for (j=0;j<3;j++){
[5149]	782	k=GetId(triangles[i].TriangleAdj(j));
[3913]	783	if (reft[k]>=0){
[6412]	784	_assert_(3num+j<3(nbt-nbtout));
[3913]	785	bamgmesh->ElementConnectivity[3*num+j]=k+1; // back to Matlab indexing
	786	}
	787	}
	788	num+=1;
	789	}
	790	}
	791
	792	/ElementNodal Connectivity/
[15104]	793	if(verbose>5) _printf_(" writing Nodal element connectivity\n");
[3913]	794	bamgmesh->NodalElementConnectivitySize[0]=nbv;
	795	bamgmesh->NodalElementConnectivitySize[1]=connectivitymax_1;
[12456]	796	bamgmesh->NodalElementConnectivity=xNew<double>(connectivitymax_1*nbv);
[3913]	797	for (i=0;i<connectivitymax_1*nbv;i++) bamgmesh->NodalElementConnectivity[i]=NAN;
	798	for (i=0;i<nbv;i++){
	799	k=0;
	800	for(j=head_1[i];j!=-1;j=next_1[j]){
[6412]	801	_assert_(connectivitymax_1i+k < connectivitymax_1nbv);
[11093]	802	bamgmesh->NodalElementConnectivity[connectivitymax_1*i+k]=floor((double)j/3)+1;
[3913]	803	k++;
	804	}
	805	}
	806
	807	/Nodal Connectivity/
[15104]	808	if(verbose>5) _printf_(" writing Nodal connectivity\n");
[3913]	809	//chaining algorithm (again...)
	810	int* head_2=NULL;
	811	int* next_2=NULL;
	812	int* connectivitysize_2=NULL;
	813	int connectivitymax_2=0;
[5091]	814	i1=bamgmesh->IssmEdgesSize[0];
	815	i2=bamgmesh->IssmEdgesSize[1];
[12456]	816	head_2=xNew<int>(nbv);
	817	next_2=xNew<int>(2*i1);
	818	connectivitysize_2=xNew<int>(nbv);
[3913]	819	//Initialization
	820	for (i=0;i<nbv;i++) head_2[i]=-1;
	821	for (i=0;i<nbv;i++) connectivitysize_2[i]=0;
	822	k=0;
	823	//Chains generation
	824	for (i=0;i<i1;i++) {
	825	for (j=0;j<2;j++){
[5091]	826	int v=(int)bamgmesh->IssmEdges[i*i2+j]-1; //back to C indexing
[13036]	827	if (k>2i1-1 \|\| k<0) _error_("Index exceed matrix dimensions (k=" << k << " not in [0 " << 2i1-1 << "]");
[3913]	828	next_2[k]=head_2[v];
[13036]	829	if (v>nbv-1 \|\| v<0) _error_("Index exceed matrix dimensions (v=" << v << " not in [0 " << nbv-1 << "])");
[3913]	830	head_2[v]=k++;
	831	connectivitysize_2[v]+=1;
	832	}
	833	}
	834	//Get maximum connectivity
	835	for (i=0;i<nbv;i++){
	836	if (connectivitysize_2[i]>connectivitymax_2) connectivitymax_2=connectivitysize_2[i];
	837	}
	838	//Build output
[14280]	839	connectivitymax_2++;//add last column for size
[3913]	840	bamgmesh->NodalConnectivitySize[0]=nbv;
	841	bamgmesh->NodalConnectivitySize[1]=connectivitymax_2;
[12456]	842	bamgmesh->NodalConnectivity=xNew<double>(connectivitymax_2*nbv);
[14280]	843	for (i=0;i<connectivitymax_2*nbv;i++) bamgmesh->NodalConnectivity[i]=0;
[3913]	844	for (i=0;i<nbv;i++){
	845	k=0;
	846	for(j=head_2[i];j!=-1;j=next_2[j]){
[6412]	847	_assert_(connectivitymax_2i+k < connectivitymax_2nbv);
[5091]	848	num=(int)bamgmesh->IssmEdges[int(j/2)*i2+0];
[3913]	849	if (i+1==num){ //carefull, ElementEdge is in M indexing
	850	//i is the first vertex of the edge, it is therefore connected to the second vertex
[5091]	851	bamgmesh->NodalConnectivity[connectivitymax_2i+k]=bamgmesh->IssmEdges[int(j/2)i2+1];
[3913]	852	}
	853	else{
	854	bamgmesh->NodalConnectivity[connectivitymax_2*i+k]=num;
	855	}
	856	k++;
	857	}
[14280]	858	bamgmesh->NodalConnectivity[connectivitymax_2*(i+1)-1]=k;
[3913]	859	}
	860
	861	/Cracked vertices/
[15104]	862	if(verbose>5) _printf_(" writing Cracked vertices\n");
[3913]	863	bamgmesh->CrackedVerticesSize[0]=NbCrackedVertices;
	864	bamgmesh->CrackedVerticesSize[1]=2;
	865	if (NbCrackedVertices){
[12456]	866	bamgmesh->CrackedVertices=xNew<double>(2*NbCrackedVertices);
[3913]	867	for (i=0;i<NbCrackedVertices;i++){
	868	bamgmesh->CrackedVertices[i2+0]=CrackedVertices[i2+0]+1; //M indexing
	869	bamgmesh->CrackedVertices[i2+1]=CrackedVertices[i2+1]+1; //M indexing
	870	}
	871	}
	872
	873	/Cracked vertices/
[15104]	874	if(verbose>5) _printf_(" writing Cracked vertices\n");
[3913]	875	bamgmesh->CrackedEdgesSize[0]=NbCrackedEdges;
	876	bamgmesh->CrackedEdgesSize[1]=4;
	877	if (NbCrackedEdges){
[12456]	878	bamgmesh->CrackedEdges=xNew<double>(2*NbCrackedEdges);
[3913]	879	for (i=0;i<NbCrackedEdges;i++){
	880	bamgmesh->CrackedEdges[i*2+0]=0;//CrackedEdges[i]->+1; //M indexing
	881	bamgmesh->CrackedEdges[i*2+1]=0;//CrackedEdges[i]-]->+1; //M indexing
	882	}
	883	}
	884
	885	//clean up
[12456]	886	xDelete<int>(connectivitysize_1);
	887	xDelete<int>(head_1);
	888	xDelete<int>(next_1);
	889	xDelete<int>(connectivitysize_2);
	890	xDelete<int>(head_2);
	891	xDelete<int>(next_2);
[3913]	892	delete [] reft;
	893	delete [] numt;
	894	}
[12365]	895	/}}}/
[18064]	896	void Mesh::ReadMetric(const BamgOpts* bamgopts) {/{{{/
[3913]	897
	898	/Intermediary/
	899	int i,j;
	900
[15104]	901	if(bamgopts->verbose>3) _printf_(" processing metric\n");
[3913]	902	double hmin = Max(bamgopts->hmin,MinimalHmin());
	903	double hmax = Min(bamgopts->hmax,MaximalHmax());
[5091]	904	double coef = bamgopts->coeff;
[3913]	905
	906	//for now we only use j==3
	907	j=3;
	908
	909	for (i=0;i<nbv;i++){
	910	double h;
	911	if (j == 1){
	912	h=bamgopts->metric[i];
	913	vertices[i].m=Metric(Max(hmin,Min(hmax, h*coef)));
	914	}
	915	else if (j==3){
	916	//do not erase metric computed by hVertices
	917	if (vertices[i].m.a11==1 && vertices[i].m.a21==0 && vertices[i].m.a22==1){
	918	double a,b,c;
	919	a=bamgopts->metric[i*3+0];
	920	b=bamgopts->metric[i*3+1];
	921	c=bamgopts->metric[i*3+2];
	922	Metric M(a,b,c);
[5400]	923	EigenMetric Vp(M/coef);
[3913]	924
	925	Vp.Maxh(hmax);
	926	Vp.Minh(hmin);
	927	vertices[i].m = Vp;
	928	}
	929	}
	930	}
	931	}
[12365]	932	/}}}/
[18064]	933	void Mesh::WriteMetric(BamgOpts* bamgopts) {/{{{/
[3913]	934	int i;
[12456]	935	xDelete<double>(bamgopts->metric);
	936	bamgopts->metric=xNew<double>(3*nbv);
[3913]	937	for (i=0;i<nbv;i++){
	938	bamgopts->metric[i*3+0]=vertices[i].m.a11;
	939	bamgopts->metric[i*3+1]=vertices[i].m.a21;
	940	bamgopts->metric[i*3+2]=vertices[i].m.a22;
	941	}
	942	}
[12365]	943	/}}}/
[18064]	944	void Mesh::WriteIndex(int** pindex,int* pnels){/{{{/
[3913]	945
[10205]	946	/Intermediary/
[13800]	947	int i,k;
[10205]	948
	949	/output/
	950	int* index=NULL;
[13800]	951	int num=0;
[10205]	952
	953	/Get number of triangles/
	954	k=0;
	955	for (i=0;i<nbt;i++){
	956	Triangle &t=triangles[i];
	957	if(t.det>0) k++;
	958	}
	959
	960	if (k){
[12456]	961	index=xNew<int>(3*k);
[10205]	962	for (i=0;i<nbt;i++){
	963	Triangle &t=triangles[i];
	964	if (t.det>0 && !(t.Hidden(0)\|\|t.Hidden(1) \|\| t.Hidden(2) )){
[16280]	965	/Remove triangles that have a bad aspect ratio/
[12408]	966	//if(t.Anisotropy()<2 & t.Length()<1.e+5){
[11584]	967	index[num*3+0]=GetId(t[0])+1; //back to M indexing
	968	index[num*3+1]=GetId(t[1])+1; //back to M indexing
	969	index[num*3+2]=GetId(t[2])+1; //back to M indexing
	970	num=num+1;
[12408]	971	//}
[10205]	972	}
	973	}
	974	}
	975
	976	/Assign output pointers/
	977	*pindex=index;
[11584]	978	*pnels=num;
[10205]	979	}
[12365]	980	/}}}/
[10205]	981
[3913]	982	/Methods/
[18064]	983	void Mesh::AddMetric(BamgOpts* bamgopts){/{{{/
[15219]	984	// Hessiantype = 0 => H is computed using double L2 projection
[3913]	985	// Hessiantype = 1 => H is computed with green formula
	986
	987	/Options/
	988	int Hessiantype=bamgopts->Hessiantype;
	989
	990	if (Hessiantype==0){
	991	BuildMetric0(bamgopts);
	992	}
	993	else if (Hessiantype==1){
	994	BuildMetric1(bamgopts);
	995	}
	996	else{
[15219]	997	_error_("Hessiantype " << Hessiantype << " not supported yet (1->use Green formula, 0-> double L2 projection)");
[3913]	998	}
	999	}
[12365]	1000	/}}}/
[22085]	1001	void Mesh::AddVertex( BamgVertex &s,Triangle* t, long long* det3){/{{{/
[3913]	1002	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Add)/
[5605]	1003	// -------------------------------
[3913]	1004	// s2
[5605]	1005	// !
	1006	// /\|\ !
	1007	// / \| \ !
	1008	// / \| \ !
	1009	// tt1 / \| \ tt0 !
	1010	// / \|s \ !
	1011	// / . \ !
	1012	// / . ` \ !
	1013	// / . ` \ !
	1014	// ---------------- !
[3913]	1015	// s0 tt2 s1
[5605]	1016	//-------------------------------
[21759]	1017
[5605]	1018	/Intermediaries/
	1019	Triangle* tt[3]; //the three triangles
[22085]	1020	long long det3local[3]; //three determinants (integer)
[5605]	1021	int nbzerodet =0; //number of zeros in det3
	1022	int izerodet=-1; //egde containing the vertex s
	1023	int iedge;
	1024
	1025	/three vertices of t/
[19293]	1026	BamgVertex* s0=t->GetVertex(0);
	1027	BamgVertex* s1=t->GetVertex(1);
	1028	BamgVertex* s2=t->GetVertex(2);
[5605]	1029
[3913]	1030	//determinant of t
[22085]	1031	long long detOld=t->det;
[3913]	1032
[5605]	1033	/* infvertexindex = index of the infinite vertex (NULL)
	1034	if no infinite vertex (NULL) infvertexindex=-1
	1035	else if v_i is infinite, infvertexindex=i*/
[19293]	1036	int infvertexindex = s0 ? ((s1? (s2?-1:2):1)):0;
[3913]	1037
	1038	//some checks
[16521]	1039	if(((infvertexindex <0 ) && (detOld <0)) \|\| ((infvertexindex >=0) && (detOld >0)) ){
[13036]	1040	_error_("inconsistent configuration (Contact ISSM developers)");
[3913]	1041	}
	1042
	1043	// if det3 does not exist, build it
[5605]	1044	if (!det3){
[3913]	1045	//allocate
	1046	det3 = det3local;
	1047	//if no infinite vertex
[5605]	1048	if (infvertexindex<0 ) {
[19293]	1049	det3[0]=bamg::det(s .GetIntegerCoordinates(),s1->GetIntegerCoordinates(),s2->GetIntegerCoordinates());
	1050	det3[1]=bamg::det(s0->GetIntegerCoordinates(),s .GetIntegerCoordinates(),s2->GetIntegerCoordinates());
	1051	det3[2]=bamg::det(s0->GetIntegerCoordinates(),s1->GetIntegerCoordinates(),s.GetIntegerCoordinates());}
[3913]	1052	else {
[5605]	1053	// one of &s1 &s2 &s0 is NULL
[19293]	1054	det3[0]= s0 ? -1 : bamg::det(s.GetIntegerCoordinates(),s1->GetIntegerCoordinates(),s2->GetIntegerCoordinates()) ;
	1055	det3[1]= s1 ? -1 : bamg::det(s0->GetIntegerCoordinates(),s.GetIntegerCoordinates(),s2->GetIntegerCoordinates()) ;
	1056	det3[2]= s2 ? -1 : bamg::det(s0->GetIntegerCoordinates(),s1->GetIntegerCoordinates(),s.GetIntegerCoordinates()) ;
[3913]	1057	}
	1058	}
	1059
	1060	if (!det3[0]) izerodet=0,nbzerodet++;
	1061	if (!det3[1]) izerodet=1,nbzerodet++;
	1062	if (!det3[2]) izerodet=2,nbzerodet++;
	1063
	1064	//if nbzerodet>0, point s is on an egde or on a vertex
[5605]	1065	if (nbzerodet>0){
	1066	/s is on an edge/
	1067	if (nbzerodet==1) {
[3913]	1068	iedge = OppositeEdge[izerodet];
[5143]	1069	AdjacentTriangle ta = t->Adj(iedge);
[3913]	1070
[5605]	1071	/*if the point is one the boundary
	1072	add the point in outside part */
	1073	if (t->det>=0){ // inside triangle
	1074	if (((Triangle*)ta)->det<0 ) {
[3913]	1075	// add in outside triangle
[5605]	1076	AddVertex(s,( Triangle *)ta);
[3913]	1077	return;
	1078	}
	1079	}
	1080	}
[5605]	1081	else{
[13036]	1082	_error_("Cannot add a vertex more than once. Check duplicates");
[3913]	1083	}
	1084	}
	1085
	1086	// remove de MarkUnSwap edge
[5605]	1087	t->SetUnMarkUnSwap(0);
	1088	t->SetUnMarkUnSwap(1);
[3913]	1089	t->SetUnMarkUnSwap(2);
	1090
	1091	tt[0]= t;
	1092	tt[1]= &triangles[nbt++];
	1093	tt[2]= &triangles[nbt++];
	1094
[13036]	1095	if (nbt>maxnbt) _error_("Not enough triangles");
[3913]	1096
[5623]	1097	tt[1]=tt[2]=*t;
[3913]	1098	tt[0]->link=tt[1];
	1099	tt[1]->link=tt[2];
	1100
[5623]	1101	(*tt[0])(OppositeVertex[0])=&s;
	1102	(*tt[1])(OppositeVertex[1])=&s;
	1103	(*tt[2])(OppositeVertex[2])=&s;
[3913]	1104
	1105	tt[0]->det=det3[0];
	1106	tt[1]->det=det3[1];
	1107	tt[2]->det=det3[2];
	1108
	1109	// update adj des triangles externe
	1110	tt[0]->SetAdjAdj(0);
	1111	tt[1]->SetAdjAdj(1);
	1112	tt[2]->SetAdjAdj(2);
	1113	// update des adj des 3 triangle interne
	1114	const int i0 = 0;
	1115	const int i1= NextEdge[i0];
	1116	const int i2 = PreviousEdge[i0];
	1117
	1118	tt[i0]->SetAdj2(i2,tt[i2],i0);
	1119	tt[i1]->SetAdj2(i0,tt[i0],i1);
	1120	tt[i2]->SetAdj2(i1,tt[i1],i2);
	1121
[5124]	1122	tt[0]->SetSingleVertexToTriangleConnectivity();
	1123	tt[1]->SetSingleVertexToTriangleConnectivity();
	1124	tt[2]->SetSingleVertexToTriangleConnectivity();
[3913]	1125
	1126	// swap if the point s is on a edge
	1127	if(izerodet>=0) {
[5623]	1128	int rswap=tt[izerodet]->swap(iedge);
[3913]	1129
	1130	if (!rswap) {
[13036]	1131	_error_("swap the point s is on a edge");
[3913]	1132	}
	1133	}
[21759]	1134
	1135	}/}}}/
[18064]	1136	void Mesh::BoundAnisotropy(double anisomax,double hminaniso) {/{{{/
[3913]	1137	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/BoundAnisotropy)/
	1138
[5091]	1139	long int verbose=0;
[3913]	1140	double lminaniso = 1/ (Max(hminaniso*hminaniso,1e-100));
	1141
	1142	//display info
[15100]	1143	if (verbose > 1) _printf_(" BoundAnisotropy by " << anisomax << "\n");
[3913]	1144
	1145	double h1=1.e30,h2=1e-30;
	1146	double coef = 1./(anisomax*anisomax);
	1147	double hn1=1.e30,hn2=1e-30,rnx =1.e-30,rx=0;
	1148
	1149	//loop over all vertices
	1150	for (int i=0;i<nbv;i++){
[5400]	1151	EigenMetric Vp(vertices[i]);
[3913]	1152	double lmax=Vp.lmax();
	1153	Vp*=Min(lminaniso,lmax)/lmax;
	1154	Vp.BoundAniso2(coef);
	1155	vertices[i].m = Vp;
	1156
	1157	//info to be displayed
[5091]	1158	if (verbose>2){
[3913]	1159	h1 =Min(h1,Vp.lmin());
	1160	h2 =Max(h2,Vp.lmax());
	1161	hn1=Min(hn1,Vp.lmin());
	1162	hn2=Max(hn2,Vp.lmax());
	1163	rx =Max(rx,Vp.Aniso2());
	1164	rnx= Max(rnx,Vp.Aniso2());
	1165	}
	1166	}
	1167
	1168	//display info
[5091]	1169	if (verbose>2){
[15100]	1170	_printf_(" input: Hmin = " << pow(h2,-0.5) << ", Hmax = " << pow(h1,-0.5) << ", factor of anisotropy max = " << pow(rx,0.5) << "\n");
	1171	_printf_(" output: Hmin = " << pow(hn2,-0.5) << ", Hmax = " << pow(hn1,-0.5)<< ", factor of anisotropy max = " <<pow(rnx,0.5) << "\n");
[3913]	1172	}
	1173	}
[12365]	1174	/}}}/
[18064]	1175	void Mesh::BuildGeometryFromMesh(BamgOpts* bamgopts){/{{{/
[3913]	1176	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/ConsGeometry)/
	1177
	1178	/Reconstruct Geometry from Mesh/
	1179
	1180	/Intermediary/
	1181	int i,j,k,kk,it,jt;
[5091]	1182	int verbose=0;
[3913]	1183
	1184	/Recover options/
[22085]	1185	if(bamgopts){
[5091]	1186	verbose=bamgopts->verbose;
[3913]	1187	}
	1188
	1189	//display info
[15104]	1190	if (verbose>1) _printf_(" construction of the geometry from the 2d mesh\n");
[3913]	1191
	1192	//check that the mesh is not empty
[22085]	1193	if(nbt<=0 \|\| nbv <=0 ) _error_("nbt or nbv is negative (Mesh empty?)");
[3913]	1194
	1195	/Construction of the edges/
	1196
	1197	//initialize st and edge4
	1198	SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv);
	1199	long* st = new long[nbt*3];
	1200
	1201	//initialize st as -1 (chaining algorithm)
	1202	for (i=0;i<nbt*3;i++) st[i]=-1;
	1203
	1204	//build edge4 (chain)
	1205	for (i=0;i<nbe;i++){
[5149]	1206	edge4->SortAndAdd(GetId(edges[i][0]),GetId(edges[i][1]));
[3913]	1207	}
	1208	//check that there is no double edge
	1209	if (nbe != edge4->nb()){
[9306]	1210	delete [] st;
[13036]	1211	_error_("Some Double edge in the mesh, the number is " << nbe << ", nbe4=" << edge4->nb());
[3913]	1212	}
	1213	//keep nbe in nbeold
	1214	long nbeold = nbe;
	1215
	1216	//Go through the triangles and ass the edges in edge4 if they are not there yet
	1217	for (i=0;i<nbt;i++){
	1218	//3 edges per triangle
	1219	for (j=0;j<3;j++) {
	1220	//Add Edge to edge4 (k=numberofedges in edge4)
[5149]	1221	long k =edge4->SortAndAdd(GetId(triangles[i][VerticesOfTriangularEdge[j][0]]), GetId(triangles[i][VerticesOfTriangularEdge[j][1]]));
[3913]	1222	long invisible = triangles[i].Hidden(j);
	1223
	1224	//if st[k] has not been changed yet, add 3*i+j (= vertex position in the index)
	1225	if(st[k]==-1) st[k]=3*i+j;
	1226
	1227	//else st[k]>=0 -> the edge already exist, check
	1228	else if(st[k]>=0) {
	1229	//check that it is not an edge on boundary (should not already exist)
	1230	if (triangles[i].TriangleAdj(j) \|\| triangles[st[k]/3].TriangleAdj((int) (st[k]%3))){
[13036]	1231	_error_("problem in Geometry reconstruction: an edge on boundary is duplicated (double element?)");
[3913]	1232	}
	1233	//OK, the element is not on boundary, is belongs to 2 triangles -> build Adjacent triangles list
	1234	triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3));
	1235	if (invisible) triangles[i].SetHidden(j);
	1236	// if k < nbe mark the edge as on Boundary (Locked)
	1237	if (k<nbe) {
	1238	triangles[i].SetLocked(j);
	1239	}
	1240	//set st[k] as negative so that the edge should not be called again
	1241	st[k]=-2-st[k];
	1242	}
	1243	//else (see 3 lines above), the edge has been called more than twice: return error
	1244	else {
[15104]	1245	_printf_("The edge (" << GetId(triangles[i][VerticesOfTriangularEdge[j][0]]) << "," << GetId(triangles[i][VerticesOfTriangularEdge[j][1]]) << ") belongs to more than 2 triangles (" << k << ")\n");
[15100]	1246	_printf_("Edge " << j << " of triangle " << i << "\n");
	1247	_printf_("Edge " << (-st[k]+2)%3 << " of triangle " << (-st[k]+2)/3 << "\n");
	1248	_printf_("Edge " << triangles[(-st[k]+2)/3].NuEdgeTriangleAdj((int)((-st[k]+2)%3)) << " of triangle " << GetId(triangles[(-st[k]+2)/3].TriangleAdj((int)((-st[k]+2)%3))) << "\n");
[13036]	1249	_error_("An edge belongs to more than 2 triangles");
[3913]	1250	}
	1251	}
	1252	}
	1253
	1254	//delete edge4
	1255	long nbedges = edge4->nb(); // the total number of edges
	1256	delete edge4; edge4=NULL;
	1257
	1258	//display info
[5091]	1259	if(verbose>5) {
[15104]	1260	_printf_(" info on Mesh:\n");
[15100]	1261	_printf_(" - number of vertices = " << nbv << "\n");
	1262	_printf_(" - number of triangles = " << nbt << "\n");
	1263	_printf_(" - number of given edges = " << nbe << "\n");
	1264	_printf_(" - number of all edges = " << nbedges << "\n");
	1265	_printf_(" - Euler number 1 - nb of holes = " << nbt-nbedges+nbv << "\n");
[3913]	1266	}
	1267
	1268	// check consistency of edge[].adj and geometrical required vertices
	1269	k=0; kk=0;
	1270	for (i=0;i<nbedges;i++){
	1271	//internal edge
	1272	if (st[i] <-1) {
	1273	//get triangle number back
	1274	it = (-2-st[i])/3;
	1275	//get edge position back
	1276	j = (int) ((-2-st[i])%3);
	1277	Triangle &tt=*triangles[it].TriangleAdj(j);
	1278	if (triangles[it].color != tt.color\|\| i < nbeold) k++;
	1279	}
	1280	//boundary edge (alone)
	1281	else if (st[i] >=0)
	1282	kk++;
	1283	}
	1284
	1285	/Constructions of edges/
	1286
	1287	k += kk;
	1288	kk=0;
	1289	if (k) {
	1290	nbe = k;
	1291	Edge* edgessave=edges;
	1292	edges = new Edge[nbe];
	1293	k =0;
	1294
	1295	//display info
[15100]	1296	if(verbose>4) _printf_(" Construction of the edges " << nbe << "\n");
[3913]	1297
	1298	for (i=0;i<nbedges;i++){
	1299	long add= -1;
	1300
	1301	//internal edge (belongs to two triangles)
	1302	if (st[i] <-1){
	1303	it = (-2-st[i])/3;
	1304	j = (int) ((-2-st[i])%3);
	1305	Triangle & tt = * triangles[it].TriangleAdj(j);
	1306	if (triangles[it].color != tt.color \|\| i < nbeold) add=k++;
	1307	}
	1308	//boundary edge
	1309	else if (st[i] >=0){
	1310	it = st[i]/3;
	1311	j = (int) (st[i]%3);
	1312	add=k++;
	1313	}
	1314	if (add>=0 && add < nbe){
	1315	edges[add].v[0] = &triangles[it][VerticesOfTriangularEdge[j][0]];
	1316	edges[add].v[1] = &triangles[it][VerticesOfTriangularEdge[j][1]];
[5573]	1317	edges[add].GeomEdgeHook=NULL;
[3913]	1318	//if already existed
	1319	if (i<nbeold){
[5148]	1320	edges[add].ReferenceNumber=edgessave[i].ReferenceNumber;
[5573]	1321	edges[add].GeomEdgeHook=edgessave[i].GeomEdgeHook; // HACK to get required edges
[15104]	1322	_printf_("oh no...\n");
[3913]	1323	}
	1324	else
[5148]	1325	edges[add].ReferenceNumber=Min(edges[add].v[0]->GetReferenceNumber(),edges[add].v[1]->GetReferenceNumber());
[3913]	1326	}
	1327	}
	1328
	1329	//check that we have been through all edges
	1330	if (k!=nbe){
[13036]	1331	_error_("problem in edge construction process: k!=nbe (should not happen)");
[3913]	1332	}
	1333	//delete edgessave
	1334	if (edgessave) delete [] edgessave;
	1335	}
	1336
	1337	/Color the vertices/
	1338
	1339	//initialize color of all vertices as 0
	1340	for (i=0;i<nbv;i++) vertices[i].color =0;
	1341
	1342	//go through the edges and add a color to corresponding vertices
	1343	//(A vertex in 4 edges will have a color 4)
	1344	for (i=0;i<nbe;i++){
	1345	for (j=0;j<2;j++) edges[i].v[j]->color++;
	1346	}
	1347
	1348	//change the color: if a vertex belongs to 2 edges -1, else -2
	1349	for (i=0;i<nbv;i++) {
	1350	vertices[i].color=(vertices[i].color ==2)? -1 : -2;
	1351	}
	1352
	1353	/Build edges[i].adj: adjacency of each edge (if on the same curve)/
	1354	for (i=0;i<nbe;i++){
	1355	for (j=0;j<2;j++){
	1356	//get current vertex
[5120]	1357	BamgVertex* v=edges[i].v[j];
[3913]	1358	//get vertex color (i0)
	1359	long i0=v->color;
	1360	long j0;
	1361
	1362	//if color<0 (first time), no adjacent edge
	1363	if(i0<0) edges[i].adj[j]=NULL;
	1364
	1365	//if color=-1 (corner),change the vertex color as 2*i+j (position of the vertex in edges)
	1366	if(i0==-1) v->color=i*2+j;
	1367
	1368	//if color>=0 (i and i0 edge are adjacent by the vertex v)
	1369	else if (i0>=0) {
	1370	//get position of v in edges back
	1371	j0 = i0%2; //column in edges
	1372	i0 = i0/2; //line in edges
	1373
	1374	//check that we have the correct vertex
	1375	if (v!=edges[i0 ].v[j0]){
[13036]	1376	_error_("v!=edges[i0 ].v[j0]: this should not happen as the vertex belongs to this edge");
[3913]	1377	}
	1378
	1379	//Add adjacence
	1380	edges[i ].adj[j ]=edges +i0;
	1381	edges[i0].adj[j0]=edges +i ;
	1382
	1383	//change color to -3
	1384	v->color = -3;
	1385	}
	1386	}
	1387	}
	1388
	1389	/Reconstruct subdomains info/
	1390
[5148]	1391	//check that nbsubdomains is empty
[22085]	1392	if(nbsubdomains) _error_("nbsubdomains should be 0");
[5148]	1393	nbsubdomains=0;
[3913]	1394
	1395	//color the subdomains
	1396	long* colorT= new long[nbt];
	1397	Triangle tt,t;
	1398
	1399	//initialize the color of each triangle as -1
	1400	for (it=0;it<nbt;it++) colorT[it]=-1;
	1401
	1402	//loop over the triangles
	1403	for (it=0;it<nbt;it++){
	1404
	1405	//if the triangle has not been colored yet:
	1406	if (colorT[it]<0){
	1407
	1408	//color = number of subdomains
[5148]	1409	colorT[it]=nbsubdomains;
[3913]	1410
	1411	//color all the adjacent triangles of T that share a non marked edge
	1412	int level =1;
	1413	int kolor=triangles[it].color;
	1414	st[0]=it; // stack
	1415	st[1]=0;
	1416	k=1;
	1417	while (level>0){
	1418	if( (j=st[level]++)<3 ){
	1419	t = &triangles[st[level-1]];
	1420	tt=t->TriangleAdj((int)j);
	1421
	1422	//color the adjacent triangle
[5149]	1423	if ( ! t->Locked(j) && tt && (colorT[jt = GetId(tt)] == -1) && ( tt->color==kolor)) {
[5148]	1424	colorT[jt]=nbsubdomains;
[3913]	1425	st[++level]=jt;
	1426	st[++level]=0;
	1427	k++;
	1428	}
	1429	}
	1430	else level-=2;
	1431	}
[5148]	1432	nbsubdomains++;
[3913]	1433	}
	1434	}
[15100]	1435	if (verbose> 3) _printf_(" The Number of sub domain = " << nbsubdomains << "\n");
[3913]	1436
	1437	//build subdomains
	1438	long isd;
[5148]	1439	subdomains = new SubDomain[nbsubdomains];
[3913]	1440
	1441	//initialize subdomains[isd].head as 0
[5148]	1442	for (isd=0;isd<nbsubdomains;isd++) subdomains[isd].head =0;
[13622]	1443
[3913]	1444	k=0;
	1445	for (it=0;it<nbt;it++){
	1446	for (int j=0;j<3;j++){
	1447	tt=triangles[it].TriangleAdj(j);
	1448	if ((!tt \|\| tt->color != triangles[it].color) && !subdomains[isd=colorT[it]].head){
	1449	subdomains[isd].head = triangles+it;
[5148]	1450	subdomains[isd].ReferenceNumber = triangles[it].color;
	1451	subdomains[isd].direction = j; // hack
[3913]	1452	subdomains[isd].edge = 0;
	1453	k++;
	1454	}
	1455	}
	1456	}
	1457	//check that we have been through all subdomains
[5148]	1458	if (k!= nbsubdomains){
[9309]	1459	delete [] colorT;
[13036]	1460	_error_("k!= nbsubdomains");
[3913]	1461	}
	1462	//delete colorT and st
	1463	delete [] colorT;
	1464	delete [] st;
	1465
	1466	/Reconstruct Geometry Gh/
	1467
	1468	//build colorV -1 for all vertex and 0 for the vertices belonging to edges
	1469	long* colorV = new long[nbv];
	1470	for (i=0;i<nbv;i++) colorV[i]=-1;
	1471	for (i=0;i<nbe;i++){
[5149]	1472	for ( j=0;j<2;j++) colorV[GetId(edges[i][j])]=0;
[3913]	1473	}
	1474	//number the vertices belonging to edges
	1475	k=0;
	1476	for (i=0;i<nbv;i++){
	1477	if(!colorV[i]) colorV[i]=k++;
	1478	}
	1479
	1480	//Build Gh
	1481	Gh.nbv=k;
	1482	Gh.nbe = nbe;
[5573]	1483	Gh.vertices = new GeomVertex[k];
	1484	Gh.edges = new GeomEdge[nbe];
[5148]	1485	Gh.nbsubdomains = nbsubdomains;
[5573]	1486	Gh.subdomains = new GeomSubDomain[nbsubdomains];
[15100]	1487	if (verbose>3) _printf_(" number of vertices = " << Gh.nbv << "\n number of edges = " << Gh.nbe << "\n");
[3913]	1488	NbVerticesOnGeomVertex = Gh.nbv;
	1489	VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];
	1490	NbVerticesOnGeomEdge =0;
	1491	VerticesOnGeomEdge =0;
	1492
	1493	//Build VertexOnGeom
	1494	for (i=0;i<nbv;i++){
	1495	if((j=colorV[i])>=0){
[5120]	1496	BamgVertex & v = Gh.vertices[j];
[3913]	1497	v = vertices[i];
	1498	v.color =0;
	1499	VerticesOnGeomVertex[j] = VertexOnGeom(vertices[i], Gh.vertices[j]);
	1500	}
	1501	}
	1502
	1503	//Buid pmin and pmax of Gh (extrema coordinates)
	1504	Gh.pmin = Gh.vertices[0].r;
	1505	Gh.pmax = Gh.vertices[0].r;
	1506	// recherche des extrema des vertices pmin,pmax
	1507	for (i=0;i<Gh.nbv;i++) {
	1508	Gh.pmin.x = Min(Gh.pmin.x,Gh.vertices[i].r.x);
	1509	Gh.pmin.y = Min(Gh.pmin.y,Gh.vertices[i].r.y);
	1510	Gh.pmax.x = Max(Gh.pmax.x,Gh.vertices[i].r.x);
	1511	Gh.pmax.y = Max(Gh.pmax.y,Gh.vertices[i].r.y);
	1512	}
	1513	R2 DD05 = (Gh.pmax-Gh.pmin)*0.05;
	1514	Gh.pmin -= DD05;
	1515	Gh.pmax += DD05;
	1516
	1517	//Build Gh.coefIcoor
[21759]	1518	int MaxICoord = 1073741823; //2^30 - 1 = =111...111 (29 times one)
	1519	Gh.coefIcoor= (MaxICoord)/(Max(Gh.pmax.x-Gh.pmin.x,Gh.pmax.y-Gh.pmin.y));
[3913]	1520	if (Gh.coefIcoor<=0){
[9306]	1521	delete [] colorV;
[13036]	1522	_error_("Gh.coefIcoor<=0 in infered Geometry (this should not happen)");
[3913]	1523	}
	1524
	1525	/Build Gh.edges/
	1526
	1527	//initialize len as 0
	1528	double * len = new double[Gh.nbv];
	1529	for(i=0;i<Gh.nbv;i++) len[i]=0;
	1530
	1531	//initialize edge4 again
	1532	edge4= new SetOfEdges4(nbe,nbv);
	1533	double hmin = HUGE_VAL;
	1534	int kreq=0;
	1535	for (i=0;i<nbe;i++){
	1536
[5149]	1537	long i0 = GetId(edges[i][0]);
	1538	long i1 = GetId(edges[i][1]);
[3913]	1539	long j0 = colorV[i0];
	1540	long j1 = colorV[i1];
	1541
	1542	Gh.edges[i].v[0] = Gh.vertices + j0;
	1543	Gh.edges[i].v[1] = Gh.vertices + j1;
	1544
[5146]	1545	Gh.edges[i].type = 0;
[3913]	1546
	1547	Gh.edges[i].tg[0]=R2();
	1548	Gh.edges[i].tg[1]=R2();
	1549
[5573]	1550	bool required= edges[i].GeomEdgeHook;
[3913]	1551	if(required) kreq++;
[5573]	1552	edges[i].GeomEdgeHook = Gh.edges + i;
[3913]	1553	if(required){
	1554	Gh.edges[i].v[0]->SetRequired();
	1555	Gh.edges[i].v[1]->SetRequired();
	1556	Gh.edges[i].SetRequired();
	1557	}
	1558
	1559	R2 x12 = Gh.vertices[j0].r-Gh.vertices[j1].r;
	1560	double l12=Norme2(x12);
	1561	hmin = Min(hmin,l12);
	1562
	1563	Gh.vertices[j1].color++;
	1564	Gh.vertices[j0].color++;
	1565
	1566	len[j0]+= l12;
	1567	len[j1] += l12;
	1568	hmin = Min(hmin,l12);
[5148]	1569	Gh.edges[i].ReferenceNumber = edges[i].ReferenceNumber;
[3913]	1570
	1571	k = edge4->SortAndAdd(i0,i1);
	1572	if (k != i){
[9306]	1573	delete [] len;
[9309]	1574	delete [] colorV;
[13036]	1575	_error_("problem in Edge4 construction: k != i");
[3913]	1576	}
	1577	}
	1578
	1579	//Build metric for all vertices of Gh
	1580	for (i=0;i<Gh.nbv;i++){
	1581	if (Gh.vertices[i].color > 0)
	1582	Gh.vertices[i].m= Metric(len[i] /(double) Gh.vertices[i].color);
	1583	else
	1584	Gh.vertices[i].m= Metric(hmin);
	1585	}
	1586	//delete len
	1587	delete [] len;
	1588
	1589	//Build Gh.subdomains
[5148]	1590	for (i=0;i<nbsubdomains;i++){
[5149]	1591	it = GetId(subdomains[i].head);
[5148]	1592	j = subdomains[i].direction;
[5149]	1593	long i0 = GetId(triangles[it][VerticesOfTriangularEdge[j][0]]);
	1594	long i1 = GetId(triangles[it][VerticesOfTriangularEdge[j][1]]);
[3913]	1595	k = edge4->SortAndFind(i0,i1);
[22085]	1596	if(k<0){
	1597	delete [] colorV;
	1598	_error_("This should not happen");
	1599	}
[16237]	1600	subdomains[i].direction = (vertices + i0 == edges[k].v[0]) ? 1 : -1;
	1601	subdomains[i].edge = edges+k;
	1602	Gh.subdomains[i].edge = Gh.edges + k;
	1603	Gh.subdomains[i].direction = subdomains[i].direction;
	1604	Gh.subdomains[i].ReferenceNumber = subdomains[i].ReferenceNumber;
	1605	}
[3913]	1606
	1607	delete edge4;
	1608	delete [] colorV;
	1609
	1610	//unset adj
	1611	for (i=0;i<nbt;i++){
[21759]	1612	for (j=0;j<3;j++){
[3913]	1613	triangles[i].SetAdj2(j,0,triangles[i].GetAllflag(j));
	1614	}
	1615	}
	1616
	1617	}
[12365]	1618	/}}}/
[18064]	1619	void Mesh::BuildMetric0(BamgOpts* bamgopts){/{{{/
[3913]	1620
	1621	/Options/
	1622	double* s=NULL;
	1623	long nbsol;
[5091]	1624	int verbose;
[3913]	1625
	1626	int i,j,k,iA,iB,iC;
	1627	int iv;
	1628
	1629	/Recover options/
[5091]	1630	verbose=bamgopts->verbose;
[3913]	1631
	1632	/Get and process fields/
	1633	s=bamgopts->field;
[5187]	1634	nbsol=bamgopts->fieldSize[1];
[3913]	1635
[5187]	1636	/Check size/
[13036]	1637	if (bamgopts->fieldSize[0] != nbv) _error_("'field' should have " << nbv << " rows");
[5187]	1638
[3913]	1639	//initialization of some variables
	1640	double* ss=(double*)s;
	1641	double sA,sB,sC;
	1642	double* detT = new double[nbt];
	1643	double* sumareas = new double[nbv];
	1644	double* alpha= new double[nbt*3];
	1645	double* beta = new double[nbt*3];
	1646	double* dx_elem = new double[nbt];
	1647	double* dy_elem = new double[nbt];
	1648	double* dx_vertex = new double[nbv];
	1649	double* dy_vertex = new double[nbv];
	1650	double* dxdx_elem = new double[nbt];
	1651	double* dxdy_elem = new double[nbt];
	1652	double* dydy_elem = new double[nbt];
	1653	double* dxdx_vertex= new double[nbv];
	1654	double* dxdy_vertex= new double[nbv];
	1655	double* dydy_vertex= new double[nbv];
	1656
	1657	//display infos
[5091]	1658	if(verbose>1) {
[15104]	1659	_printf_(" Construction of Metric: number of field: " << nbsol << " (nbt=" << nbt << ", nbv=" << nbv << ")\n");
[3913]	1660	}
	1661
	1662	//first, build the chains that will be used for the Hessian computation, as weel as the area of each element
	1663	int* head_s=NULL;
[12456]	1664	head_s=xNew<int>(nbv);
[3913]	1665	int* next_p=NULL;
[12456]	1666	next_p=xNew<int>(3*nbt);
[3913]	1667	int p=0;
	1668	//initialization
	1669	for(i=0;i<nbv;i++){
	1670	sumareas[i]=0;
	1671	head_s[i]=-1;
	1672	}
	1673	for(i=0;i<nbt;i++){
	1674
	1675	//lopp over the real triangles (no boundary elements)
	1676	if(triangles[i].link){
	1677
	1678	//get current triangle t
	1679	const Triangle &t=triangles[i];
	1680
	1681	// coor of 3 vertices
	1682	R2 A=t[0];
	1683	R2 B=t[1];
	1684	R2 C=t[2];
	1685
	1686	//compute triangle determinant (2*Area)
	1687	double dett = bamg::Area2(A,B,C);
	1688	detT[i]=dett;
	1689
	1690	/*The nodal functions are such that for a vertex A:
	1691	* N_A(x,y)=alphaA x + beta_A y +gamma_A
	1692	* N_A(A) = 1, N_A(B) = 0, N_A(C) = 0
	1693	* solving this system of equation (determinant = 2Area(T) != 0 if A,B and C are not inlined)
	1694	* leads to:
	1695	* N_A = (xB yC - xC yB + x(yB-yC) +y(xC-xB))/(2*Area(T))
	1696	* and this gives:
	1697	* alpha_A = (yB-yC)/(2Area(T))/
	1698	alpha[i*3+0]=(B.y-C.y)/dett;
	1699	alpha[i*3+1]=(C.y-A.y)/dett;
	1700	alpha[i*3+2]=(A.y-B.y)/dett;
	1701	beta[ i*3+0]=(C.x-B.x)/dett;
	1702	beta[ i*3+1]=(A.x-C.x)/dett;
	1703	beta[ i*3+2]=(B.x-A.x)/dett;
	1704
	1705	//compute chains
	1706	for(j=0;j<3;j++){
[5149]	1707	k=GetId(triangles[i][j]);
[3913]	1708	next_p[p]=head_s[k];
	1709	head_s[k]=p++;
	1710
	1711	//add area to sumareas
	1712	sumareas[k]+=dett;
	1713	}
	1714
	1715	}
	1716	}
	1717
	1718	//for all Solutions
	1719	for (int nusol=0;nusol<nbsol;nusol++) {
	1720	double smin=ss[nusol],smax=ss[nusol];
	1721
	1722	//get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy)
	1723	for ( iv=0,k=0; iv<nbv; iv++){
	1724	smin=Min(smin,ss[iv*nbsol+nusol]);
	1725	smax=Max(smax,ss[iv*nbsol+nusol]);
	1726	}
	1727	double sdelta=smax-smin;
	1728	double absmax=Max(Abs(smin),Abs(smax));
	1729
	1730	//display info
[15100]	1731	if(verbose>2) _printf_(" Solution " << nusol << ", Min = " << smin << ", Max = " << smax << ", Delta = " << sdelta << "\n");
[3913]	1732
	1733	//skip constant field
	1734	if (sdelta < 1.0e-10*Max(absmax,1e-20)){
[15104]	1735	_printf_(" Solution " << nusol << " is constant, skipping...\n");
[3913]	1736	continue;
	1737	}
	1738
	1739	//initialize the hessian and gradient matrices
	1740	for ( iv=0,k=0; iv<nbv; iv++) dxdx_vertex[iv]=dxdy_vertex[iv]=dydy_vertex[iv]=dx_vertex[iv]=dy_vertex[iv]=0;
	1741
	1742	//1: Compute gradient for each element (exact)
	1743	for (i=0;i<nbt;i++){
	1744	if(triangles[i].link){
	1745	// number of the 3 vertices
[5149]	1746	iA = GetId(triangles[i][0]);
	1747	iB = GetId(triangles[i][1]);
	1748	iC = GetId(triangles[i][2]);
[3913]	1749
	1750	// value of the P1 fonction on 3 vertices
	1751	sA = ss[iA*nbsol+nusol];
	1752	sB = ss[iB*nbsol+nusol];
	1753	sC = ss[iC*nbsol+nusol];
	1754
	1755	//gradient = (sum alpha_i s_i, sum_i beta_i s_i)
	1756	dx_elem[i]=sAalpha[3i+0]+sBalpha[3i+1]+sCalpha[3i+2];
	1757	dy_elem[i]=sAbeta[ 3i+0]+sBbeta[ 3i+1]+sCbeta[ 3i+2];
	1758	}
	1759	}
	1760
	1761	//2: then compute a gradient for each vertex using a P2 projection
	1762	for(i=0;i<nbv;i++){
	1763	for(p=head_s[i];p!=-1;p=next_p[p]){
	1764	//Get triangle number
	1765	k=(long)(p/3);
	1766	dx_vertex[i]+=dx_elem[k]*detT[k]/sumareas[i];
	1767	dy_vertex[i]+=dy_elem[k]*detT[k]/sumareas[i];
	1768	}
	1769	}
	1770
	1771	//3: compute Hessian matrix on each element
	1772	for (i=0;i<nbt;i++){
	1773	if(triangles[i].link){
	1774	// number of the 3 vertices
[5149]	1775	iA = GetId(triangles[i][0]);
	1776	iB = GetId(triangles[i][1]);
	1777	iC = GetId(triangles[i][2]);
[3913]	1778
	1779	//Hessian
	1780	dxdx_elem[i]=dx_vertex[iA]alpha[3i+0]+dx_vertex[iB]alpha[3i+1]+dx_vertex[iC]alpha[3i+2];
	1781	dxdy_elem[i]=dy_vertex[iA]alpha[3i+0]+dy_vertex[iB]alpha[3i+1]+dy_vertex[iC]alpha[3i+2];
	1782	dydy_elem[i]=dy_vertex[iA]beta[3i+0]+dy_vertex[iB]beta[3i+1]+dy_vertex[iC]beta[3i+2];
	1783	}
	1784	}
	1785
	1786	//4: finaly compute Hessian on each vertex using the second P2 projection
	1787	for(i=0;i<nbv;i++){
	1788	for(p=head_s[i];p!=-1;p=next_p[p]){
	1789	//Get triangle number
	1790	k=(long)(p/3);
	1791	dxdx_vertex[i]+=dxdx_elem[k]*detT[k]/sumareas[i];
	1792	dxdy_vertex[i]+=dxdy_elem[k]*detT[k]/sumareas[i];
	1793	dydy_vertex[i]+=dydy_elem[k]*detT[k]/sumareas[i];
	1794	}
	1795	}
	1796
	1797	/Compute Metric from Hessian/
	1798	for ( iv=0;iv<nbv;iv++){
	1799	vertices[iv].MetricFromHessian(dxdx_vertex[iv],dxdy_vertex[iv],dydy_vertex[iv],smin,smax,ss[iv*nbsol+nusol],bamgopts->err[nusol],bamgopts);
	1800	}
	1801
	1802	}//for all solutions
	1803
	1804	//clean up
[12456]	1805	xDelete<int>(head_s);
	1806	xDelete<int>(next_p);
[3913]	1807	delete [] detT;
	1808	delete [] alpha;
	1809	delete [] beta;
	1810	delete [] sumareas;
	1811	delete [] dx_elem;
	1812	delete [] dy_elem;
	1813	delete [] dx_vertex;
	1814	delete [] dy_vertex;
	1815	delete [] dxdx_elem;
	1816	delete [] dxdy_elem;
	1817	delete [] dydy_elem;
	1818	delete [] dxdx_vertex;
	1819	delete [] dxdy_vertex;
	1820	delete [] dydy_vertex;
	1821	}
[12365]	1822	/}}}/
[18064]	1823	void Mesh::BuildMetric1(BamgOpts* bamgopts){/{{{/
[3913]	1824	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/IntersectConsMetric)/
	1825
	1826	/Options/
	1827	double* s=NULL;
	1828	long nbsol;
	1829	int NbJacobi;
[5091]	1830	int verbose;
[3913]	1831
	1832	/Recover options/
[5091]	1833	verbose=bamgopts->verbose;
[3913]	1834	NbJacobi=bamgopts->nbjacobi;
	1835
	1836	/Get and process fields/
	1837	s=bamgopts->field;
[5187]	1838	nbsol=bamgopts->fieldSize[1];
[3913]	1839
[5187]	1840	/Check size/
[13036]	1841	if (bamgopts->fieldSize[0] != nbv) _error_("'field' should have " << nbv << " rows");
[5187]	1842
[3913]	1843	//initialization of some variables
	1844	long i,k,iA,iB,iC,iv;
	1845	R2 O(0,0);
	1846	double* ss=(double*)s;
	1847	double sA,sB,sC;
	1848	double* detT = new double[nbt];
	1849	double* Mmass= new double[nbv];
	1850	double* Mmassxx= new double[nbv];
	1851	double* dxdx= new double[nbv];
	1852	double* dxdy= new double[nbv];
	1853	double* dydy= new double[nbv];
	1854	double* workT= new double[nbt];
	1855	double* workV= new double[nbv];
	1856	int* OnBoundary = new int[nbv];
	1857
	1858	//display infos
[5091]	1859	if(verbose>1) {
[15104]	1860	_printf_(" Construction of Metric: number of field: " << nbsol << " (nbt=" << nbt << ", nbv=" << nbv << ")\n");
[3913]	1861	}
	1862
	1863	//initialize Mmass, OnBoundary and Massxx by zero
	1864	for (iv=0;iv<nbv;iv++){
	1865	Mmass[iv]=0;
	1866	OnBoundary[iv]=0;
	1867	Mmassxx[iv]=0;
	1868	}
	1869
	1870	//Build detT Mmas Mmassxx workT and OnBoundary
	1871	for (i=0;i<nbt;i++){
	1872
	1873	//lopp over the real triangles (no boundary elements)
	1874	if(triangles[i].link){
	1875
	1876	//get current triangle t
	1877	const Triangle &t=triangles[i];
	1878
	1879	// coor of 3 vertices
	1880	R2 A=t[0];
	1881	R2 B=t[1];
	1882	R2 C=t[2];
	1883
	1884	// number of the 3 vertices
[5149]	1885	iA = GetId(t[0]);
	1886	iB = GetId(t[1]);
	1887	iC = GetId(t[2]);
[3913]	1888
	1889	//compute triangle determinant (2*Area)
	1890	double dett = bamg::Area2(A,B,C);
	1891	detT[i]=dett;
	1892	dett /= 6;
	1893
	1894	// construction of OnBoundary (flag=1 if on boundary, else 0)
	1895	int nbb=0;
	1896	for(int j=0;j<3;j++){
	1897	//get adjacent triangle
	1898	Triangle *ta=t.Adj(j);
	1899	//if there is no adjacent triangle, the edge of the triangle t is on boundary
	1900	if ( !ta \|\| !ta->link){
	1901	//mark the two vertices of the edge as OnBoundary
[5149]	1902	OnBoundary[GetId(t[VerticesOfTriangularEdge[j][0]])]=1;
	1903	OnBoundary[GetId(t[VerticesOfTriangularEdge[j][1]])]=1;
[3913]	1904	nbb++;
	1905	}
	1906	}
	1907
	1908	//number of vertices on boundary for current triangle t
	1909	workT[i] = nbb;
	1910
	1911	//Build Mmass Mmass[i] = Mmass[i] + Area/3
	1912	Mmass[iA] += dett;
	1913	Mmass[iB] += dett;
	1914	Mmass[iC] += dett;
	1915
	1916	//Build Massxx = Mmass
	1917	Mmassxx[iA] += dett;
	1918	Mmassxx[iB] += dett;
	1919	Mmassxx[iC] += dett;
	1920	}
	1921
	1922	//else: the triangle is a boundary triangle -> workT=-1
	1923	else workT[i]=-1;
	1924	}
	1925
	1926	//for all Solution
	1927	for (int nusol=0;nusol<nbsol;nusol++) {
	1928
	1929	double smin=ss[nusol],smax=ss[nusol];
	1930	double h1=1.e30,h2=1e-30,rx=0;
	1931	double hn1=1.e30,hn2=1e-30,rnx =1.e-30;
	1932
	1933	//get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy)
	1934	for ( iv=0,k=0; iv<nbv; iv++ ){
	1935	dxdx[iv]=dxdy[iv]=dydy[iv]=0;
	1936	smin=Min(smin,ss[iv*nbsol+nusol]);
	1937	smax=Max(smax,ss[iv*nbsol+nusol]);
	1938	}
	1939	double sdelta=smax-smin;
	1940	double absmax=Max(Abs(smin),Abs(smax));
	1941
	1942	//display info
[15100]	1943	if(verbose>2) _printf_(" Solution " << nusol << ", Min = " << smin << ", Max = " << smax << ", Delta = " << sdelta << ", number of fields = " << nbsol << "\n");
[3913]	1944
	1945	//skip constant field
	1946	if (sdelta < 1.0e-10*Max(absmax,1e-20) ){
[15104]	1947	if (verbose>2) _printf_(" Solution " << nusol << " is constant, skipping...\n");
[3913]	1948	continue;
	1949	}
	1950
	1951	//pointer toward ss that is also a pointer toward s (solutions)
	1952	double* sf=ss;
	1953
	1954	//initialize the hessian matrix
	1955	for ( iv=0,k=0; iv<nbv; iv++) dxdx[iv]=dxdy[iv]=dydy[iv]=0;
	1956
	1957	//loop over the triangles
	1958	for (i=0;i<nbt;i++){
	1959
	1960	//for real all triangles
	1961	if(triangles[i].link){
	1962
	1963	// coor of 3 vertices
	1964	R2 A=triangles[i][0];
	1965	R2 B=triangles[i][1];
	1966	R2 C=triangles[i][2];
	1967
	1968	//warning: the normal is internal and the size is the length of the edge
	1969	R2 nAB = Orthogonal(B-A);
	1970	R2 nBC = Orthogonal(C-B);
	1971	R2 nCA = Orthogonal(A-C);
	1972	//note that : nAB + nBC + nCA == 0
	1973
	1974	// number of the 3 vertices
[5149]	1975	iA = GetId(triangles[i][0]);
	1976	iB = GetId(triangles[i][1]);
	1977	iC = GetId(triangles[i][2]);
[3913]	1978
	1979	// for the test of boundary edge
	1980	// the 3 adj triangles
	1981	Triangle *tBC = triangles[i].TriangleAdj(OppositeEdge[0]);
	1982	Triangle *tCA = triangles[i].TriangleAdj(OppositeEdge[1]);
	1983	Triangle *tAB = triangles[i].TriangleAdj(OppositeEdge[2]);
	1984
	1985	// value of the P1 fonction on 3 vertices
	1986	sA = ss[iA*nbsol+nusol];
	1987	sB = ss[iB*nbsol+nusol];
	1988	sC = ss[iC*nbsol+nusol];
	1989
	1990	/*The nodal functions are such that for a vertex A:
	1991	N_A(x,y)=alphaA x + beta_A y +gamma_A
	1992	N_A(A) = 1, N_A(B) = 0, N_A(C) = 0
	1993	solving this system of equation (determinant = 2Area(T) != 0 if A,B and C are not inlined)
	1994	leads to:
	1995	N_A = (xB yC - xC yB + x(yB-yC) +y(xC-xB))/(2*Area(T))
	1996	and this gives:
	1997	alpha_A = (yB-yC)/(2*Area(T))
	1998	beta_A = (xC-xB)/(2*Area(T))
	1999	and therefore:
	2000	grad N_A = nA / detT
	2001	for an interpolation of a solution s:
	2002	grad(s) = s * sum_{i=A,B,C} grad(N_i) */
	2003
	2004	R2 Grads=(nABsC+nBCsA+nCA*sB)/detT[i];
	2005
	2006	//Use Green to compute Hessian Matrix
	2007
	2008	// if edge on boundary no contribution => normal = 0
	2009	if ( !tBC \|\| !tBC->link ) nBC=O;
	2010	if ( !tCA \|\| !tCA->link ) nCA=O;
	2011	if ( !tAB \|\| !tAB->link ) nAB=O;
	2012
	2013	// remark we forgot a 1/2 because
	2014	// int_{edge} w_i = 1/2 if i is in edge
	2015	// 0 if not
	2016	// if we don't take the boundary
	2017	dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x;
	2018	dxdx[iB] += ( nAB.x + nBC.x ) *Grads.x;
	2019	dxdx[iC] += ( nBC.x + nCA.x ) *Grads.x;
	2020
	2021	//warning optimization (1) the division by 2 is done on the metric construction
	2022	dxdy[iA] += (( nCA.y + nAB.y ) Grads.x + ( nCA.x + nAB.x ) Grads.y) ;
	2023	dxdy[iB] += (( nAB.y + nBC.y ) Grads.x + ( nAB.x + nBC.x ) Grads.y) ;
	2024	dxdy[iC] += (( nBC.y + nCA.y ) Grads.x + ( nBC.x + nCA.x ) Grads.y) ;
	2025
	2026	dydy[iA] += ( nCA.y + nAB.y ) *Grads.y;
	2027	dydy[iB] += ( nAB.y + nBC.y ) *Grads.y;
	2028	dydy[iC] += ( nBC.y + nCA.y ) *Grads.y;
	2029
	2030	} // for real all triangles
	2031	}
	2032
	2033	long kk=0;
	2034	for ( iv=0,k=0 ; iv<nbv; iv++){
	2035	if(Mmassxx[iv]>0){
	2036	dxdx[iv] /= 2*Mmassxx[iv];
	2037	// warning optimization (1) on term dxdy[iv]*ci/2
	2038	dxdy[iv] /= 4*Mmassxx[iv];
	2039	dydy[iv] /= 2*Mmassxx[iv];
	2040	// Compute the matrix with abs(eigen value)
	2041	Metric M(dxdx[iv], dxdy[iv], dydy[iv]);
[5400]	2042	EigenMetric Vp(M);
[3913]	2043	Vp.Abs();
	2044	M = Vp;
	2045	dxdx[iv] = M.a11;
	2046	dxdy[iv] = M.a21;
	2047	dydy[iv] = M.a22;
	2048	}
	2049	else kk++;
	2050	}
	2051
	2052	// correction of second derivative
	2053	// by a laplacien
	2054	double* dd;
	2055	for (int xy = 0;xy<3;xy++) {
[16156]	2056	if (xy==0) dd=dxdx;
	2057	else if (xy==1) dd=dxdy;
	2058	else if (xy==2) dd=dydy;
[22085]	2059	else{
	2060	delete [] detT;
	2061	delete [] Mmass;
	2062	delete [] workT;
	2063	delete [] workV;
	2064	delete [] Mmassxx;
	2065	delete [] OnBoundary;
	2066	_error_("not supported yet");
	2067	}
[3913]	2068	// do leat 2 iteration for boundary problem
	2069	for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){
	2070	for (i=0;i<nbt;i++)
	2071	if(triangles[i].link){// the real triangles
	2072	// number of the 3 vertices
[5149]	2073	iA = GetId(triangles[i][0]);
	2074	iB = GetId(triangles[i][1]);
	2075	iC = GetId(triangles[i][2]);
[3913]	2076	double cc=3;
	2077	if(ijacobi==0)
	2078	cc = Max((double) ((Mmassxx[iA]>0)+(Mmassxx[iB]>0)+(Mmassxx[iC]>0)),1.);
	2079	workT[i] = (dd[iA]+dd[iB]+dd[iC])/cc;
	2080	}
	2081	for (iv=0;iv<nbv;iv++) workV[iv]=0;
	2082
	2083	for (i=0;i<nbt;i++){
	2084	if(triangles[i].link){ // the real triangles
	2085	// number of the 3 vertices
[5149]	2086	iA = GetId(triangles[i][0]);
	2087	iB = GetId(triangles[i][1]);
	2088	iC = GetId(triangles[i][2]);
[3913]	2089	double cc = workT[i]*detT[i];
	2090	workV[iA] += cc;
	2091	workV[iB] += cc;
	2092	workV[iC] += cc;
	2093	}
	2094	}
	2095
	2096	for (iv=0;iv<nbv;iv++){
	2097	if( ijacobi<NbJacobi \|\| OnBoundary[iv]){
	2098	dd[iv] = workV[iv]/(Mmass[iv]*6);
	2099	}
	2100	}
	2101	}
	2102	}
	2103
	2104	/Compute Metric from Hessian/
	2105	for ( iv=0;iv<nbv;iv++){
	2106	vertices[iv].MetricFromHessian(dxdx[iv],dxdy[iv],dydy[iv],smin,smax,ss[iv*nbsol+nusol],bamgopts->err[nusol],bamgopts);
	2107	}
	2108
	2109	}// end for all solution
	2110
	2111	delete [] detT;
	2112	delete [] Mmass;
	2113	delete [] dxdx;
	2114	delete [] dxdy;
	2115	delete [] dydy;
[22085]	2116	delete [] workT;
[3913]	2117	delete [] workV;
	2118	delete [] Mmassxx;
[22085]	2119	delete [] OnBoundary;
[3913]	2120
	2121	}
[12365]	2122	/}}}/
[18064]	2123	void Mesh::CrackMesh(BamgOpts* bamgopts) {/{{{/
[3913]	2124	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/CrackMesh)/
	2125
	2126	/Intermediary/
	2127	int i,j,k,num,count;
	2128	int i1,i2;
	2129	int j1,j2;
	2130
	2131	/Options/
	2132	int verbose=bamgopts->verbose;
	2133
	2134	// computed the number of cracked edge
	2135	for (k=i=0;i<nbe;i++){
[5573]	2136	if(edges[i].GeomEdgeHook->Cracked()) k++;
[3913]	2137	}
	2138
	2139	//Return if no edge is cracked
	2140	if(k==0) return;
[15100]	2141	if (verbose>4) _printf_(" number of Cracked Edges = " << k << "\n");
[3913]	2142
	2143	//Initialize Cracked edge
	2144	NbCrackedEdges=k;
	2145	CrackedEdges=new CrackedEdge[k];
	2146
	2147	//Compute number of Cracked Vertices
	2148	k=0;
	2149	NbCrackedVertices=0;
	2150
	2151	int* splitvertex=new int[nbv];
	2152	for (i=0;i<nbv;i++) splitvertex[i]=0;
	2153
	2154	for (i=0;i<nbe;i++){
[5573]	2155	if(edges[i].GeomEdgeHook->Cracked()){
[3913]	2156
	2157	//Fill edges fields of CrackedEdges
[5573]	2158	CrackedEdges[k ].E =edges[i].GeomEdgeHook;
[3913]	2159	CrackedEdges[k++].e1=&edges[i];
	2160
	2161	//Get number of the two vertices on the edge
[5149]	2162	i1=GetId(edges[i][0]);
	2163	i2=GetId(edges[i][1]);
[6412]	2164	_assert_(i1>=0 && i1<nbv && i2>=0 && i2<nbv);
[3913]	2165	splitvertex[i1]++;
	2166	splitvertex[i2]++;
	2167
	2168	//If the vertex has already been flagged once, it is a cracked vertex (tip otherwise)
	2169	if (splitvertex[i1]==2) NbCrackedVertices++;
	2170	if (splitvertex[i2]==2) NbCrackedVertices++;
	2171
	2172	//The vertex cannot be marked more than twice
	2173	if (splitvertex[i1]==3 \|\| splitvertex[i2]==3){
[9306]	2174	delete [] splitvertex;
[13036]	2175	_error_("Crossing rifts not supported yet");
[3913]	2176	}
	2177	}
	2178	}
[6412]	2179	_assert_(k==NbCrackedEdges);
[3913]	2180
	2181	//Add new vertices
[15100]	2182	if (verbose>4) _printf_(" number of Cracked Vertices = " << NbCrackedVertices << "\n");
[3913]	2183	if (NbCrackedVertices){
[12456]	2184	CrackedVertices=xNew<long>(2*NbCrackedVertices);
[3913]	2185	num=0;
	2186	for (i=0;i<nbv;i++){
	2187	if (splitvertex[i]==2){
	2188	CrackedVertices[num*2+0]=i; //index of first vertex
	2189	CrackedVertices[num*2+1]=nbv+num;//index of new vertex
	2190	num++;
	2191	}
	2192	}
[6412]	2193	_assert_(num==NbCrackedVertices);
[3913]	2194	}
	2195	delete [] splitvertex;
	2196
	2197	//Now, find the triangles that hold a cracked edge
[5124]	2198	CreateSingleVertexToTriangleConnectivity();
[3913]	2199
	2200	long* Edgeflags=new long[NbCrackedEdges];
	2201	for(i=0;i<NbCrackedEdges;i++) Edgeflags[i]=0;
	2202
	2203	for(i=0;i<NbCrackedEdges;i++){
	2204	//Get the numbers of the 2 vertices of the crren cracked edge
[5149]	2205	i1=GetId((*CrackedEdges[i].e1)[0]);
	2206	i2=GetId((*CrackedEdges[i].e1)[1]);
[3913]	2207
	2208	//find a triangle holding the vertex i1 (first vertex of the ith cracked edge)
	2209	Triangle* tbegin=vertices[i1].t;
[5460]	2210	k=vertices[i1].IndexInTriangle;//local number of i in triangle tbegin
[6412]	2211	_assert_(GetId((*tbegin)[k])==GetId(vertices[i1]));
[3913]	2212
	2213	//Now, we are going to go through the adjacent triangle that hold i1 till
	2214	//we find one that has the cracked edge
[5143]	2215	AdjacentTriangle ta(tbegin,EdgesVertexTriangle[k][0]);
[3913]	2216	count=0;
	2217	do {
	2218	for(j=0;j<3;j++){
	2219	//Find the position of i1 in the triangle index
[5149]	2220	if (GetId((*ta.t)[j])==i1){
[3913]	2221	j1=j;
	2222	break;
	2223	}
	2224	}
	2225	for(j=0;j<3;j++){
	2226	//Check wether i2 is also in the triangle index
[5149]	2227	if (GetId((*ta.t)[j])==i2){
[3913]	2228	j2=j;
	2229	//Invert j1 and j2 if necessary
	2230	if ((j1+1)%3==j2){
	2231	int j3=j1;
	2232	j1=j2;
	2233	j2=j3;
	2234	}
	2235	if (Edgeflags[i]==0){
	2236	//first element
	2237	CrackedEdges[i].a=ta.t;
	2238	CrackedEdges[i].length=Norme2((ta.t)[j1].r-(ta.t)[j2].r);
	2239	CrackedEdges[i].normal=Orthogonal((ta.t)[j1].r-(ta.t)[j2].r);
	2240	}
	2241	else{
	2242	//Second element -> to renumber
	2243	CrackedEdges[i].b=ta.t;
	2244	CrackedEdges[i].length=Norme2((ta.t)[j1].r-(ta.t)[j2].r);
	2245	CrackedEdges[i].normal=Orthogonal((ta.t)[j1].r-(ta.t)[j2].r);
	2246	}
	2247	Edgeflags[i]++;
	2248	break;
	2249	}
	2250	}
[15104]	2251	//_printf_(element_renu[GetId(ta.t)] << " -> " << GetId((ta.t)[0])+1 << " " << GetId((ta.t)[1])+1 << " " << GetId((*ta.t)[2])+1 << ", edge [" << i1 << "->" << j1 << " " << i2 << "->" << j2 << "]\n");
[3913]	2252	ta = Next(ta).Adj();
[13036]	2253	if (count++>50) _error_("Maximum number of iteration exceeded");
[3913]	2254	}while ((tbegin != ta));
	2255	}
	2256
	2257	//Check EdgeFlag
	2258	for(i=0;i<NbCrackedEdges;i++){
	2259	if (Edgeflags[i]!=2){
[13036]	2260	_error_("A problem occured: at least one crack edge (number " << i+1 << ") does not belong to 2 elements");
[3913]	2261	}
	2262	}
	2263	delete [] Edgeflags;
	2264
[5120]	2265	//Reset BamgVertex to On
[3913]	2266	SetVertexFieldOn();
	2267
	2268	}
[12365]	2269	/}}}/
[18064]	2270	void Mesh::Echo(void) {/{{{/
[5208]	2271
	2272	int i;
	2273
[15104]	2274	_printf_("Mesh Echo:\n");
[15100]	2275	_printf_(" nbv = " << nbv << "\n");
	2276	_printf_(" nbt = " << nbt << "\n");
	2277	_printf_(" nbe = " << nbe << "\n");
[15104]	2278	_printf_(" index:\n");
[5208]	2279	for (i=0;i<nbt;i++){
[15100]	2280	_printf_(" " << setw(4) << i+1 << ": ["
[21759]	2281	<< setw(4) << (((BamgVertex*)triangles[i](0))?GetId(triangles[i][0])+1:0) << " "
	2282	<< setw(4) << (((BamgVertex*)triangles[i](1))?GetId(triangles[i][1])+1:0) << " "
	2283	<< setw(4) << (((BamgVertex*)triangles[i](2))?GetId(triangles[i][2])+1:0) << "]\n");
[5208]	2284	}
[15104]	2285	_printf_(" coordinates:\n");
[5208]	2286	for (i=0;i<nbv;i++){
[21759]	2287	_printf_(" " << setw(4) << i+1 << ": [" << vertices[i].r.x << " " << vertices[i].r.y << "] and [" << vertices[i].i.x << " " << vertices[i].i.y << "]\n");
[5208]	2288	}
	2289
	2290	}
[12365]	2291	/}}}/
[18064]	2292	void Mesh::ForceBoundary() {/{{{/
	2293	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ForceBoundary)/
[3913]	2294
[18064]	2295	long int verbose=2;
	2296	int k=0;
	2297	int nbfe=0,nbswp=0,Nbswap=0;
[3913]	2298
[18064]	2299	//display
	2300	if (verbose > 2) _printf_(" ForceBoundary nb of edge: " << nbe << "\n");
[3913]	2301
[18064]	2302	//check that there is no triangle with 0 determinant
	2303	for (int t = 0; t < nbt; t++){
	2304	if (!triangles[t].det) k++;
	2305	}
	2306	if (k!=0) {
	2307	_error_("there is " << k << " triangles of mes = 0");
	2308	}
[3913]	2309
[18064]	2310	//Force Edges
	2311	AdjacentTriangle ta(0,0);
	2312	for (int i = 0; i < nbe; i++){
[3913]	2313
[18064]	2314	//Force edge i
	2315	nbswp = ForceEdge(edges[i][0],edges[i][1],ta);
	2316	if (nbswp<0) k++;
	2317	else Nbswap += nbswp;
[3913]	2318
[18064]	2319	if (nbswp) nbfe++;
	2320	if ( nbswp < 0 && k < 5){
	2321	_error_("Missing Edge " << i << ", v0=" << GetId(edges[i][0]) << ",v1=" << GetId(edges[i][1]));
[3913]	2322	}
[18064]	2323	}
[3913]	2324
[18064]	2325	if (k!=0) {
	2326	_error_("There are " << k << " lost edges, the boundary might be crossing");
[3913]	2327	}
[18064]	2328	for (int j=0;j<nbv;j++){
	2329	Nbswap += vertices[j].Optim(1,0);
	2330	}
	2331	if (verbose > 3) _printf_(" number of inforced edge = " << nbfe << ", number of swap= " << Nbswap << "\n");
	2332	}
[12365]	2333	/}}}/
[18064]	2334	void Mesh::FindSubDomain(int OutSide) {/{{{/
[3913]	2335	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FindSubDomain)/
	2336
[5091]	2337	long int verbose=0;
[3913]	2338
[5091]	2339	if (verbose >2){
[15104]	2340	if (OutSide) _printf_(" Find all external sub-domain\n");
	2341	else _printf_(" Find all internal sub-domain\n");
[3913]	2342	}
	2343	short * HeapArete = new short[nbt];
	2344	Triangle ** HeapTriangle = new Triangle* [nbt];
	2345	Triangle t,t1;
	2346	long k,it;
	2347
	2348	for (int itt=0;itt<nbt;itt++)
	2349	triangles[itt].link=0; // par defaut pas de couleur
	2350
	2351	long NbSubDomTot =0;
	2352	for ( it=0;it<nbt;it++) {
	2353	if ( ! triangles[it].link ) {
	2354	t = triangles + it;
	2355	NbSubDomTot++;; // new composante connexe
	2356	long i = 0; // niveau de la pile
	2357	t->link = t ; // sd forme d'un triangle cicular link
	2358
	2359	HeapTriangle[i] =t ;
	2360	HeapArete[i] = 3;
	2361
	2362	while (i >= 0) // boucle sur la pile
	2363	{ while ( HeapArete[i]--) // boucle sur les 3 aretes
	2364	{
	2365	int na = HeapArete[i];
	2366	Triangle * tc = HeapTriangle[i]; // triangle courant
	2367	if( ! tc->Locked(na)) // arete non frontiere
	2368	{
	2369	Triangle * ta = tc->TriangleAdj(na) ; // næ triangle adjacent
	2370	if (ta->link == 0 ) // non deja chainer => on enpile
	2371	{
	2372	i++;
	2373	ta->link = t->link ; // on chaine les triangles
	2374	t->link = ta ; // d'un meme sous domaine
	2375	HeapArete[i] = 3; // pour les 3 triangles adjacents
	2376	HeapTriangle[i] = ta;
	2377	}}
	2378	} // deplie fin de boucle sur les 3 adjacences
	2379	i--;
	2380	}
	2381	}
	2382	}
	2383
	2384	// supression de tous les sous domaine infini <=> contient le sommet NULL
	2385	it =0;
[5146]	2386	nbtout = 0;
[3913]	2387	while (it<nbt) {
	2388	if (triangles[it].link)
	2389	{
	2390	if (!( triangles[it](0) && triangles[it](1) && triangles[it](2) ))
	2391	{
	2392	// infini triangle
	2393	NbSubDomTot --;
	2394	t=&triangles[it];
[5146]	2395	nbtout--; // on fait un coup de trop.
[3913]	2396	while (t){
[5146]	2397	nbtout++;
[3913]	2398	t1=t;
	2399	t=t->link;
	2400	t1->link=0;
	2401	}
	2402	}
	2403	}
	2404	it++;} // end while (it<nbt)
[5146]	2405	if (nbt == nbtout \|\| !NbSubDomTot) {
[9306]	2406	delete [] HeapArete;
[22085]	2407	delete [] HeapTriangle;
[13036]	2408	_error_("The boundary is not close: all triangles are outside");
[3913]	2409	}
	2410
	2411	delete [] HeapArete;
	2412	delete [] HeapTriangle;
	2413
[5148]	2414	if (OutSide\|\| !Gh.subdomains \|\| !Gh.nbsubdomains )
[3913]	2415	{ // No geom sub domain
	2416	long i;
	2417	if (subdomains) delete [] subdomains;
	2418	subdomains = new SubDomain[ NbSubDomTot];
[5148]	2419	nbsubdomains= NbSubDomTot;
	2420	for ( i=0;i<nbsubdomains;i++) {
[3913]	2421	subdomains[i].head=NULL;
[5148]	2422	subdomains[i].ReferenceNumber=i+1;
[3913]	2423	}
	2424	long * mark = new long[nbt];
	2425	for (it=0;it<nbt;it++)
	2426	mark[it]=triangles[it].link ? -1 : -2;
	2427
	2428	it =0;
	2429	k = 0;
	2430	while (it<nbt) {
	2431	if (mark[it] == -1) {
	2432	t1 = & triangles[it];
	2433	t = t1->link;
	2434	mark[it]=k;
	2435	subdomains[k].head = t1;
	2436	do {
[5149]	2437	mark[GetId(t)]=k;
[3913]	2438	t=t->link;
	2439	} while (t!=t1);
	2440	mark[it]=k++;}
	2441	// else if(mark[it] == -2 ) triangles[it].Draw(999);
	2442	it++;} // end white (it<nbt)
[5148]	2443	if (k!=nbsubdomains){
[9306]	2444	delete [] mark;
[13036]	2445	_error_("k!=nbsubdomains");
[3913]	2446	}
	2447	if(OutSide)
	2448	{
	2449	// to remove all the sub domain by parity adjacents
	2450	// because in this case we have only the true boundary edge
	2451	// so teh boundary is manifold
[5148]	2452	long nbk = nbsubdomains;
[3913]	2453	while (nbk)
	2454	for (it=0;it<nbt && nbk ;it++)
	2455	for (int na=0;na<3 && nbk ;na++)
	2456	{
	2457	Triangle *ta = triangles[it].TriangleAdj(na);
[5149]	2458	long kl = ta ? mark[GetId(ta)] : -2;
[3913]	2459	long kr = mark[it];
	2460	if(kr !=kl) {
[5148]	2461	if (kl >=0 && subdomains[kl].ReferenceNumber <0 && kr >=0 && subdomains[kr].ReferenceNumber>=0)
	2462	nbk--,subdomains[kr].ReferenceNumber=subdomains[kl].ReferenceNumber-1;
	2463	if (kr >=0 && subdomains[kr].ReferenceNumber <0 && kl >=0 && subdomains[kl].ReferenceNumber>=0)
	2464	nbk--,subdomains[kl].ReferenceNumber=subdomains[kr].ReferenceNumber-1;
	2465	if(kr<0 && kl >=0 && subdomains[kl].ReferenceNumber>=0)
	2466	nbk--,subdomains[kl].ReferenceNumber=-1;
	2467	if(kl<0 && kr >=0 && subdomains[kr].ReferenceNumber>=0)
	2468	nbk--,subdomains[kr].ReferenceNumber=-1;
[3913]	2469	}
	2470	}
	2471	long j=0;
[5148]	2472	for ( i=0;i<nbsubdomains;i++)
	2473	if((-subdomains[i].ReferenceNumber) %2) { // good
[3913]	2474	if(i != j)
	2475	Exchange(subdomains[i],subdomains[j]);
	2476	j++;}
	2477	else{
	2478	t= subdomains[i].head;
	2479	while (t){
[5146]	2480	nbtout++;
[3913]	2481	t1=t;
	2482	t=t->link;
	2483	t1->link=0;
	2484	}//while (t)
	2485	}
[15100]	2486	if(verbose>4) _printf_(" Number of removes subdomains (OutSideMesh) = " << nbsubdomains-j << "\n");
[5148]	2487	nbsubdomains=j;
[3913]	2488	}
	2489
	2490	delete [] mark;
	2491
	2492	}
[16166]	2493	else{ // find the head for all subdomains
[5148]	2494	if (Gh.nbsubdomains != nbsubdomains && subdomains)
[3913]	2495	delete [] subdomains, subdomains=0;
	2496	if (! subdomains )
[5148]	2497	subdomains = new SubDomain[ Gh.nbsubdomains];
	2498	nbsubdomains =Gh.nbsubdomains;
[5124]	2499	CreateSingleVertexToTriangleConnectivity();
[3913]	2500	long * mark = new long[nbt];
[5573]	2501	Edge **GeomEdgetoEdge = MakeGeomEdgeToEdge();
[3913]	2502
	2503	for (it=0;it<nbt;it++)
	2504	mark[it]=triangles[it].link ? -1 : -2;
	2505	long inew =0;
[5148]	2506	for (int i=0;i<nbsubdomains;i++) {
[5573]	2507	GeomEdge &eg = *Gh.subdomains[i].edge;
[5148]	2508	subdomains[i].ReferenceNumber = Gh.subdomains[i].ReferenceNumber;
[5573]	2509	// by carefull is not easy to find a edge create from a GeomEdge
	2510	// see routine MakeGeomEdgeToEdge
	2511	Edge &e = *GeomEdgetoEdge[Gh.GetId(eg)];
[6412]	2512	_assert_(&e);
[5120]	2513	BamgVertex * v0 = e(0),*v1 = e(1);
[3913]	2514	Triangle *t = v0->t;
[5148]	2515	int direction = Gh.subdomains[i].direction;
	2516	// test if ge and e is in the same direction
	2517	if (((eg[0].r-eg[1].r),(e[0].r-e[1].r))<0) direction = -direction ;
	2518	subdomains[i].direction = direction;
[3913]	2519	subdomains[i].edge = &e;
[6412]	2520	_assert_(t && direction);
[3913]	2521
[5460]	2522	AdjacentTriangle ta(t,EdgesVertexTriangle[v0->IndexInTriangle][0]);// previous edges
[3913]	2523
	2524	while (1) {
[6412]	2525	_assert_(v0==ta.EdgeVertex(1));
[3913]	2526	if (ta.EdgeVertex(0) == v1) { // ok we find the edge
[5148]	2527	if (direction>0)
[3913]	2528	subdomains[i].head=t=Adj(ta);
	2529	else
	2530	subdomains[i].head=t=ta;
	2531	if(t<triangles \|\| t >= triangles+nbt \|\| t->det < 0 \|\| t->link == 0) {
[13036]	2532	_error_("bad definition of SubSomain " << i);
[3913]	2533	}
[5149]	2534	long it = GetId(t);
[3913]	2535	if (mark[it] >=0) {
	2536	break;
	2537	}
	2538	if(i != inew)
	2539	Exchange(subdomains[i],subdomains[inew]);
	2540	inew++;
	2541	Triangle *tt=t;
	2542	long kkk=0;
	2543	do
	2544	{
	2545	kkk++;
[5149]	2546	if (mark[GetId(tt)]>=0){
[13036]	2547	_error_("mark[GetId(tt)]>=0");
[3913]	2548	}
[5149]	2549	mark[GetId(tt)]=i;
[3913]	2550	tt=tt->link;
	2551	} while (tt!=t);
	2552	break;
	2553	}
	2554	ta = Previous(Adj(ta));
	2555	if(t == (Triangle *) ta) {
[13036]	2556	_error_("bad definition of SubSomain " << i);
[3913]	2557	}
	2558	}
	2559	}
	2560
[5148]	2561	if (inew < nbsubdomains) {
[15104]	2562	if (verbose>5) _printf_("WARNING: " << nbsubdomains-inew << " SubDomains are being removed\n");
[5148]	2563	nbsubdomains=inew;}
[3913]	2564
	2565	for (it=0;it<nbt;it++)
	2566	if ( mark[it] ==-1 )
[5146]	2567	nbtout++,triangles[it].link =0;
[5573]	2568	delete [] GeomEdgetoEdge;
[3913]	2569	delete [] mark;
	2570
	2571	}
[5146]	2572	nbtout=0;
[3913]	2573	for (it=0;it<nbt;it++)
[5146]	2574	if(!triangles[it].link) nbtout++;
[3913]	2575	}
[12365]	2576	/}}}/
[18064]	2577	long Mesh::GetId(const Triangle & t) const { /{{{/
[5150]	2578	return &t - triangles;
	2579	}
[12365]	2580	/}}}/
[18064]	2581	long Mesh::GetId(const Triangle * t) const { /{{{/
[5150]	2582	return t - triangles;
	2583	}
[12365]	2584	/}}}/
[18064]	2585	long Mesh::GetId(const BamgVertex & t) const { /{{{/
[5150]	2586	return &t - vertices;
	2587	}
[12365]	2588	/}}}/
[18064]	2589	long Mesh::GetId(const BamgVertex * t) const { /{{{/
[5150]	2590	return t - vertices;
	2591	}
[12365]	2592	/}}}/
[18064]	2593	long Mesh::GetId(const Edge & t) const { /{{{/
[5150]	2594	return &t - edges;
	2595	}
[12365]	2596	/}}}/
[18064]	2597	long Mesh::GetId(const Edge * t) const { /{{{/
[5150]	2598	return t - edges;
	2599	}
[12365]	2600	/}}}/
[18064]	2601	void Mesh::Init(long maxnbv_in) {/{{{/
[5120]	2602	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/PreInit)/
[3913]	2603
[5447]	2604	/Initialize fields/
[22085]	2605	this->NbRef = 0;
	2606	this->quadtree = NULL;
	2607	this->nbv = 0;
	2608	this->nbt = 0;
	2609	this->nbe = 0;
	2610	this->edges = NULL;
	2611	this->nbsubdomains = 0;
	2612	this->subdomains = NULL;
	2613	this->maxnbv = maxnbv_in;
	2614	this->maxnbt = 2 *maxnbv_in-2;
	2615	this->NbVertexOnBThVertex = 0;
	2616	this->VertexOnBThVertex = NULL;
	2617	this->NbVertexOnBThEdge = 0;
	2618	this->VertexOnBThEdge = NULL;
	2619	this->NbCrackedVertices = 0;
	2620	this->CrackedVertices = NULL;
	2621	this->NbCrackedEdges = 0;
	2622	this->CrackedEdges = NULL;
	2623	this->NbVerticesOnGeomVertex = 0;
	2624	this->VerticesOnGeomVertex = NULL;
	2625	this->NbVerticesOnGeomEdge = 0;
	2626	this->VerticesOnGeomEdge = NULL;
[3913]	2627
[22085]	2628	/Initialize random seed/
	2629	this->randomseed = 1;
	2630
[5447]	2631	/Allocate if maxnbv_in>0/
[18453]	2632	if(maxnbv_in){
[22085]	2633	this->vertices=new BamgVertex[this->maxnbv];
	2634	this->orderedvertices=new BamgVertex* [this->maxnbv];
	2635	this->triangles=new Triangle[this->maxnbt];
	2636	_assert_(this->vertices);
	2637	_assert_(this->orderedvertices);
	2638	_assert_(this->triangles);
[3913]	2639	}
[18453]	2640	else{
[22085]	2641	this->vertices = NULL;
	2642	this->orderedvertices = NULL;
	2643	this->triangles = NULL;
	2644	this->maxnbt = 0;
[5120]	2645	}
[3913]	2646	}
[12365]	2647	/}}}/
[22085]	2648	void Mesh::Insert(){/{{{/
[3913]	2649	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Insert)/
	2650
	2651	/Insert points in the existing Geometry/
	2652
	2653	//Intermediary
	2654	int i;
	2655
	2656	/Get options/
[5091]	2657	long int verbose=2;
[3913]	2658
	2659	//Display info
[15104]	2660	if (verbose>2) _printf_(" Insert initial " << nbv << " vertices\n");
[3913]	2661
[5581]	2662	//Compute integer coordinates for the existing vertices
[3913]	2663	SetIntCoor();
	2664
[5583]	2665	/*Now we want to build a list (orderedvertices) of the vertices in a random
[3913]	2666	* order. To do so, we use the following method:
	2667	*
	2668	* From an initial k0 in [0 nbv[ random (vertex number)
	2669	* the next k (vertex number) is computed using a big
	2670	* prime number (PN>>nbv) following:
	2671	*
	2672	* k_{i+1} = k_i + PN [nbv]
	2673	*
	2674	* let's show that:
	2675	*
[21759]	2676	* for all j in [0 nbv[, there is a unique i in [0 nbv[ such that k_i=j
[3913]	2677	*
	2678	* Let's assume that there are 2 distinct j1 and j2 such that
	2679	* k_j1 = k_j2
	2680	*
	2681	* This means that
	2682	*
	2683	* k0+j1PN = k0+j2PN [nbv]
	2684	* (j1-j2)*PN =0 [nbv]
	2685	* since PN is a prime number larger than nbv, and nbv!=1
	2686	* j1-j2=0 [nbv]
	2687	* BUT
	2688	* j1 and j2 are in [0 nbv[ which is impossible.
	2689	*
	2690	* We hence have built a random list of nbv elements of
	2691	* [0 nbv[ all distincts*/
[5581]	2692
	2693	//Get Prime number
	2694	const long PrimeNumber= BigPrimeNumber(nbv);
[22085]	2695	int k0=this->RandomNumber(nbv);
[5581]	2696
[5583]	2697	//Build orderedvertices
[5581]	2698	for (i=0; i<nbv; i++){
[5583]	2699	orderedvertices[i]=&vertices[k0=(k0+PrimeNumber)%nbv];
[3913]	2700	}
	2701
[5583]	2702	/Modify orderedvertices such that the first 3 vertices form a triangle/
[3913]	2703
	2704	//get first vertex i such that [0,1,i] are not aligned
[5583]	2705	for (i=2; det(orderedvertices[0]->i,orderedvertices[1]->i,orderedvertices[i]->i)==0;){
[3913]	2706	//if i is higher than nbv, it means that all the determinants are 0,
	2707	//all vertices are aligned!
[13036]	2708	if (++i>=nbv) _error_("all the vertices are aligned");
[3913]	2709	}
[5583]	2710	// exchange i et 2 in "orderedvertices" so that
[3913]	2711	// the first 3 vertices are not aligned (real triangle)
[5583]	2712	Exchange(orderedvertices[2], orderedvertices[i]);
[3913]	2713
	2714	/*Take the first edge formed by the first two vertices and build
	2715	* the initial simple mesh from this edge and 2 boundary triangles*/
	2716
[5605]	2717	BamgVertex v0=orderedvertices[0], v1=orderedvertices[1];
[3913]	2718
	2719	nbt = 2;
[5605]	2720	triangles[0](0) = NULL;//infinite vertex
[3913]	2721	triangles[0](1) = v0;
	2722	triangles[0](2) = v1;
	2723	triangles[1](0) = NULL;//infinite vertex
	2724	triangles[1](2) = v0;
	2725	triangles[1](1) = v1;
	2726
	2727	//Build adjacence
	2728	const int e0 = OppositeEdge[0];
	2729	const int e1 = NextEdge[e0];
	2730	const int e2 = PreviousEdge[e0];
	2731	triangles[0].SetAdj2(e0, &triangles[1] ,e0);
	2732	triangles[0].SetAdj2(e1, &triangles[1] ,e2);
	2733	triangles[0].SetAdj2(e2, &triangles[1] ,e1);
	2734
	2735	triangles[0].det = -1; //boundary triangle: det = -1
	2736	triangles[1].det = -1; //boundary triangle: det = -1
	2737
[5124]	2738	triangles[0].SetSingleVertexToTriangleConnectivity();
	2739	triangles[1].SetSingleVertexToTriangleConnectivity();
[3913]	2740
	2741	triangles[0].link=&triangles[1];
	2742	triangles[1].link=&triangles[0];
	2743
	2744	//build quadtree
[12218]	2745	if (!quadtree) quadtree = new BamgQuadtree(this,0);
[3913]	2746	quadtree->Add(*v0);
	2747	quadtree->Add(*v1);
	2748
	2749	/Now, add the vertices One by One/
	2750	long NbSwap=0;
[15104]	2751	if (verbose>3) _printf_(" Begining of insertion process...\n");
[3913]	2752
	2753	for (int icount=2; icount<nbv; icount++) {
	2754
	2755	//Get new vertex
[5583]	2756	BamgVertex *newvertex=orderedvertices[icount];
[3913]	2757
	2758	//Find the triangle in which newvertex is located
[22085]	2759	long long det3[3];
[5605]	2760	Triangle* tcvi = TriangleFindFromCoord(newvertex->i,det3); //(newvertex->i = integer coordinates)
[3913]	2761
	2762	//Add newvertex to the quadtree
	2763	quadtree->Add(*newvertex);
	2764
	2765	//Add newvertex to the existing mesh
[5605]	2766	AddVertex(*newvertex,tcvi,det3);
[3913]	2767
	2768	//Make the mesh Delaunay around newvertex by swaping the edges
	2769	NbSwap += newvertex->Optim(1,0);
	2770	}
	2771
	2772	//Display info
[5091]	2773	if (verbose>3) {
[15100]	2774	_printf_(" NbSwap of insertion: " << NbSwap << "\n");
	2775	_printf_(" NbSwap/nbv: " << NbSwap/nbv << "\n");
[3913]	2776	}
	2777	}
[12365]	2778	/}}}/
[22085]	2779	long Mesh::InsertNewPoints(long nbvold,long & NbTSwap) {/{{{/
[3913]	2780	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/InsertNewPoints)/
	2781
[5091]	2782	long int verbose=0;
[3913]	2783	double seuil= 1.414/2 ;// for two close point
	2784	long i;
	2785	long NbSwap=0;
[22085]	2786	long long det3[3];
[3913]	2787
	2788	//number of new points
	2789	const long nbvnew=nbv-nbvold;
	2790
	2791	//display info if required
[15104]	2792	if (verbose>5) _printf_(" Try to Insert " << nbvnew << " new points\n");
[3913]	2793
	2794	//return if no new points
	2795	if (!nbvnew) return 0;
	2796
	2797	/construction of a random order/
	2798	const long PrimeNumber= BigPrimeNumber(nbv) ;
[22085]	2799	long k3 = long(this->RandomNumber(nbvnew));
[3913]	2800	//loop over the new points
	2801	for (int is3=0; is3<nbvnew; is3++){
[17507]	2802	long j=nbvold +(k3 = (k3+PrimeNumber)%nbvnew);
	2803	long i=nbvold+is3;
[5583]	2804	orderedvertices[i]= vertices + j;
	2805	orderedvertices[i]->ReferenceNumber=i;
[3913]	2806	}
	2807
	2808	// for all the new point
	2809	long iv=nbvold;
[18822]	2810	for(i=nbvold;i<nbv;i++){
[5583]	2811	BamgVertex &vi=*orderedvertices[i];
[5180]	2812	vi.i=R2ToI2(vi.r);
	2813	vi.r=I2ToR2(vi.i);
[3913]	2814	double hx,hy;
	2815	vi.m.Box(hx,hy);
[21759]	2816	int hi=(int) (hxcoefIcoor),hj=(int) (hycoefIcoor);
	2817	if(!quadtree->TooClose(&vi,seuil,hi,hj)){
[3913]	2818	// a good new point
[5120]	2819	BamgVertex &vj = vertices[iv];
[3913]	2820	long j=vj.ReferenceNumber;
[5583]	2821	if (&vj!=orderedvertices[j]){
[13036]	2822	_error_("&vj!= orderedvertices[j]");
[3913]	2823	}
	2824	if(i!=j){
	2825	Exchange(vi,vj);
[5583]	2826	Exchange(orderedvertices[j],orderedvertices[i]);
[3913]	2827	}
	2828	vj.ReferenceNumber=0;
[5605]	2829	Triangle *tcvj=TriangleFindFromCoord(vj.i,det3);
[3913]	2830	if (tcvj && !tcvj->link){
[18488]	2831	_printf_("While trying to add the following point:\n");
	2832	vj.Echo();
	2833	_printf_("BAMG determined that it was inside the following triangle, which probably lies outside of the geometric domain\n");
[3913]	2834	tcvj->Echo();
[13036]	2835	_error_("problem inserting point in InsertNewPoints (tcvj=" << tcvj << " and tcvj->link=" << tcvj->link << ")");
[3913]	2836	}
	2837	quadtree->Add(vj);
[5605]	2838	AddVertex(vj,tcvj,det3);
[3913]	2839	NbSwap += vj.Optim(1);
	2840	iv++;
	2841	}
[18822]	2842	else{
	2843	vi.PreviousNumber = 0;
	2844	}
[3913]	2845	}
[5091]	2846	if (verbose>3) {
[15100]	2847	_printf_(" number of new points: " << iv << "\n");
	2848	_printf_(" number of to close (?) points: " << nbv-iv << "\n");
	2849	_printf_(" number of swap after: " << NbSwap << "\n");
[3913]	2850	}
	2851	nbv = iv;
	2852
	2853	for (i=nbvold;i<nbv;i++) NbSwap += vertices[i].Optim(1);
[15100]	2854	if (verbose>3) _printf_(" NbSwap=" << NbSwap << "\n");
[3913]	2855
	2856	NbTSwap += NbSwap ;
	2857	return nbv-nbvold;
	2858	}
[12365]	2859	/}}}/
[18064]	2860	Edge** Mesh::MakeGeomEdgeToEdge() {/{{{/
[5573]	2861	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeGeomEdgeToEdge)/
[3913]	2862
	2863	if (!Gh.nbe){
[13036]	2864	_error_("!Gh.nbe");
[3913]	2865	}
[16521]	2866	Edge *e= new Edge [Gh.nbe];
[3913]	2867
	2868	long i;
	2869	for ( i=0;i<Gh.nbe ; i++)
	2870	e[i]=NULL;
	2871	for ( i=0;i<nbe ; i++)
	2872	{
	2873	Edge * ei = edges+i;
[5573]	2874	GeomEdge *GeomEdgeHook = ei->GeomEdgeHook;
	2875	e[Gh.GetId(GeomEdgeHook)] = ei;
[3913]	2876	}
	2877	for ( i=0;i<nbe ; i++)
	2878	for (int ii=0;ii<2;ii++) {
	2879	Edge * ei = edges+i;
[5573]	2880	GeomEdge *GeomEdgeHook = ei->GeomEdgeHook;
[3913]	2881	int j= ii;
[5573]	2882	while (!(*GeomEdgeHook)[j].Required()) {
	2883	Adj(GeomEdgeHook,j); // next geom edge
[3913]	2884	j=1-j;
[5573]	2885	if (e[Gh.GetId(GeomEdgeHook)]) break; // optimisation
	2886	e[Gh.GetId(GeomEdgeHook)] = ei;
[3913]	2887	}
	2888	}
	2889
	2890	int kk=0;
	2891	for ( i=0;i<Gh.nbe ; i++){
	2892	if (!e[i]){
	2893	kk++;
[15104]	2894	if(kk<10) _printf_("BUG: the geometrical edge " << i << " is on no edge curve\n");
[3913]	2895	}
	2896	}
[22085]	2897	if(kk){
	2898	delete [] e;
	2899	_error_("See above");
	2900	}
[3913]	2901
	2902	return e;
	2903	}
[12365]	2904	/}}}/
[18064]	2905	void Mesh::MakeBamgQuadtree() { /{{{/
[12218]	2906	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeBamgQuadtree)/
[16166]	2907	if(!quadtree) quadtree = new BamgQuadtree(this);
[3913]	2908	}
[12365]	2909	/}}}/
[18064]	2910	double Mesh::MaximalHmax() {/{{{/
[5150]	2911	return Max(pmax.x-pmin.x,pmax.y-pmin.y);
	2912	}
[12365]	2913	/}}}/
[18064]	2914	void Mesh::MaxSubDivision(double maxsubdiv) {/{{{/
[3913]	2915	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/MaxSubDivision)/
	2916
[5091]	2917	long int verbose=0;
[3913]	2918
	2919	const double maxsubdiv2 = maxsubdiv*maxsubdiv;
[15100]	2920	if(verbose>1) _printf_(" Limit the subdivision of a edges in the new mesh by " << maxsubdiv << "\n");
[3913]	2921	// for all the edges
	2922	// if the len of the edge is to long
	2923	long it,nbchange=0;
	2924	double lmax=0;
	2925	for (it=0;it<nbt;it++){
	2926	Triangle &t=triangles[it];
	2927	for (int j=0;j<3;j++){
[20621]	2928	Triangle *ptt=t.TriangleAdj(j);
	2929	Triangle &tt = *ptt;
	2930	if ( (!ptt \|\| it < GetId(tt)) && ( tt.link \|\| t.link)){
[5120]	2931	BamgVertex &v0 = t[VerticesOfTriangularEdge[j][0]];
	2932	BamgVertex &v1 = t[VerticesOfTriangularEdge[j][1]];
[3913]	2933	R2 AB= (R2) v1-(R2) v0;
	2934	Metric M = v0;
	2935	double l = M(AB,AB);
	2936	lmax = Max(lmax,l);
	2937	if(l> maxsubdiv2){
[21759]	2938	R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M
[3913]	2939	double lc = M(AC,AC);
	2940	D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;
	2941	D2xD2 Rt1(Rt.inv());
	2942	D2xD2 D(maxsubdiv2,0,0,lc);
	2943	D2xD2 MM = Rt1DRt1.t();
	2944	v0.m = M = Metric(MM.x.x,MM.y.x,MM.y.y);
	2945	nbchange++;
	2946	}
	2947	M = v1;
	2948	l = M(AB,AB);
	2949	lmax = Max(lmax,l);
	2950	if(l> maxsubdiv2){
[21759]	2951	R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M
[3913]	2952	double lc = M(AC,AC);
	2953	D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;
	2954	D2xD2 Rt1(Rt.inv());
	2955	D2xD2 D(maxsubdiv2,0,0,lc);
	2956	D2xD2 MM = Rt1DRt1.t();
	2957	v1.m = M = Metric(MM.x.x,MM.y.x,MM.y.y);
	2958	nbchange++;
	2959	}
	2960	}
	2961	}
	2962	}
[5091]	2963	if(verbose>3){
[15100]	2964	_printf_(" number of metric changes = " << nbchange << ", maximum number of subdivision of a edges before change = " << pow(lmax,0.5) << "\n");
[3913]	2965	}
	2966	}
[12365]	2967	/}}}/
[22085]	2968	Metric Mesh::MetricAt(const R2 & A){ /{{{/
[3913]	2969	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MetricAt)/
	2970
[5180]	2971	I2 a = R2ToI2(A);
[22085]	2972	long long deta[3];
[5124]	2973	Triangle * t =TriangleFindFromCoord(a,deta);
[3913]	2974	if (t->det <0) { // outside
	2975	double ba,bb;
[5143]	2976	AdjacentTriangle edge= CloseBoundaryEdge(a,t,ba,bb) ;
[3913]	2977	return Metric(ba,edge.EdgeVertex(0),bb,edge.EdgeVertex(1));}
	2978	else { // inside
	2979	double aa[3];
	2980	double s = deta[0]+deta[1]+deta[2];
	2981	aa[0]=deta[0]/s;
	2982	aa[1]=deta[1]/s;
	2983	aa[2]=deta[2]/s;
	2984	return Metric(aa,(t)[0],(t)[1],(*t)[2]);
	2985	}
	2986	}
[12365]	2987	/}}}/
[18064]	2988	double Mesh::MinimalHmin() {/{{{/
[5150]	2989	return 2.0/coefIcoor;
	2990	}
[12365]	2991	/}}}/
[21759]	2992	BamgVertex* Mesh::NearestVertex(int i,int j) {/{{{/
[18064]	2993	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NearestVertex)/
	2994	return quadtree->NearestVertex(i,j);
	2995	}
	2996	/}}}/
	2997	void Mesh::NewPoints(Mesh & Bh,BamgOpts* bamgopts,int KeepVertices){/{{{/
[3913]	2998	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewPoints)/
	2999
	3000	int i,j,k;
	3001	long NbTSwap=0;
	3002	long nbtold=nbt;
	3003	long nbvold=nbv;
	3004	long Headt=0;
	3005	long next_t;
[5146]	3006	long* first_np_or_next_t=new long[maxnbt];
[3913]	3007	Triangle* t=NULL;
	3008
	3009	/Recover options/
[5091]	3010	int verbose=bamgopts->verbose;
[3913]	3011
	3012	/First, insert old points if requested/
[18486]	3013	if(KeepVertices && (&Bh != this) && (nbv+Bh.nbv< maxnbv)){
[15104]	3014	if (verbose>5) _printf_(" Inserting initial mesh points\n");
[18486]	3015	bool pointsoutside = false;
[18488]	3016	for(i=0;i<Bh.nbv;i++){
[5120]	3017	BamgVertex &bv=Bh[i];
[18486]	3018	/Do not insert if the point is outside/
	3019	long long det3[3];
	3020	Triangle* tcvj=TriangleFindFromCoord(bv.i,det3);
[18488]	3021	if(tcvj->det<0 \|\| !tcvj->link){
[18486]	3022	pointsoutside = true;
	3023	continue;
	3024	}
[18488]	3025	IssmDouble area_1=((bv.r.x -(tcvj)(2)->r.x)((tcvj)(1)->r.y-(tcvj)(2)->r.y)
	3026	- (bv.r.y -(tcvj)(2)->r.y)((tcvj)(1)->r.x-(tcvj)(2)->r.x));
	3027	IssmDouble area_2=(((tcvj)(0)->r.x -(tcvj)(2)->r.x)(bv.r.y -(tcvj)(2)->r.y)
	3028	- ((tcvj)(0)->r.y -(tcvj)(2)->r.y)(bv.r.x -(tcvj)(2)->r.x));
	3029	IssmDouble area_3 =((bv.r.x -(tcvj)(1)->r.x)((tcvj)(0)->r.y-(tcvj)(1)->r.y)
	3030	- (bv.r.y -(tcvj)(1)->r.y)((tcvj)(0)->r.x-(tcvj)(1)->r.x));
	3031	if(area_1<0 \|\| area_2<0 \|\| area_3<0){
	3032	pointsoutside = true;
	3033	continue;
	3034	}
	3035	if(!bv.GeomEdgeHook){
[18455]	3036	vertices[nbv].r = bv.r;
	3037	vertices[nbv].PreviousNumber = i+1;
[3913]	3038	vertices[nbv++].m = bv.m;
	3039	}
	3040	}
[22085]	3041	//if(pointsoutside) _printf_("WARNING: One or more points of the initial mesh fall outside of the geometric boundary\n");
[5124]	3042	Bh.CreateSingleVertexToTriangleConnectivity();
[22085]	3043	InsertNewPoints(nbvold,NbTSwap);
[18488]	3044	}
[5124]	3045	else Bh.CreateSingleVertexToTriangleConnectivity();
[3913]	3046
	3047	// generation of the list of next Triangle
	3048	for(i=0;i<nbt;i++) first_np_or_next_t[i]=-(i+1);
	3049	// the next traingle of i is -first_np_or_next_t[i]
	3050
	3051	// Big loop (most time consuming)
	3052	int iter=0;
[15104]	3053	if (verbose>5) _printf_(" Big loop\n");
[3913]	3054	do {
	3055	/Update variables/
	3056	iter++;
	3057	nbtold=nbt;
	3058	nbvold=nbv;
	3059
	3060	/We test all triangles/
	3061	i=Headt;
	3062	next_t=-first_np_or_next_t[i];
	3063	for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i]){
	3064
	3065	//check i
	3066	if (i<0 \|\| i>=nbt ){
[13036]	3067	_error_("Index problem in NewPoints (i=" << i << " not in [0 " << nbt-1 << "])");
[3913]	3068	}
	3069	//change first_np_or_next_t[i]
	3070	first_np_or_next_t[i] = iter;
	3071
	3072	//Loop over the edges of t
	3073	for(j=0;j<3;j++){
[5143]	3074	AdjacentTriangle tj(t,j);
[5120]	3075	BamgVertex &vA = *tj.EdgeVertex(0);
	3076	BamgVertex &vB = *tj.EdgeVertex(1);
[3913]	3077
	3078	//if t is a boundary triangle, or tj locked, continue
	3079	if (!t->link) continue;
	3080	if (t->det <0) continue;
	3081	if (t->Locked(j)) continue;
	3082
[5143]	3083	AdjacentTriangle tadjj = t->Adj(j);
[3913]	3084	Triangle* ta=tadjj;
	3085
	3086	//if the adjacent triangle is a boundary triangle, continur
	3087	if (ta->det<0) continue;
	3088
	3089	R2 A=vA;
	3090	R2 B=vB;
[5149]	3091	k=GetId(ta);
[3913]	3092
	3093	//if this edge has already been done, go to next edge of triangle
	3094	if(first_np_or_next_t[k]==iter) continue;
	3095
	3096	lIntTria.SplitEdge(Bh,A,B);
[5120]	3097	lIntTria.NewPoints(vertices,nbv,maxnbv);
[18064]	3098	} // end loop for each edge
	3099	}// for triangle
[3913]	3100
[22085]	3101	if (!InsertNewPoints(nbvold,NbTSwap)) break;
[3913]	3102	for (i=nbtold;i<nbt;i++) first_np_or_next_t[i]=iter;
	3103	Headt = nbt; // empty list
	3104
	3105	// for all the triangle containing the vertex i
	3106	for (i=nbvold;i<nbv;i++){
[5120]	3107	BamgVertex* s = vertices + i;
[5460]	3108	AdjacentTriangle ta(s->t, EdgesVertexTriangle[s->IndexInTriangle][1]);
[3913]	3109	Triangle* tbegin= (Triangle*) ta;
	3110	long kt;
	3111	do {
[5149]	3112	kt = GetId((Triangle*) ta);
[3913]	3113	if (first_np_or_next_t[kt]>0){
	3114	first_np_or_next_t[kt]=-Headt;
	3115	Headt=kt;
	3116	}
	3117	if (ta.EdgeVertex(0)!=s){
[13036]	3118	_error_("ta.EdgeVertex(0)!=s");
[3913]	3119	}
	3120	ta = Next(Adj(ta));
	3121	} while ( (tbegin != (Triangle*) ta));
	3122	}
	3123
[16166]	3124	}while(nbv!=nbvold);
	3125	delete [] first_np_or_next_t;
[3913]	3126
[16166]	3127	long NbSwapf =0;
	3128	for(i=0;i<nbv;i++) NbSwapf += vertices[i].Optim(0);
	3129	}/}}}/
[18064]	3130	GeomEdge* Mesh::ProjectOnCurve( Edge & BhAB, BamgVertex & vA, BamgVertex & vB,/{{{/
[5120]	3131	double theta,BamgVertex & R,VertexOnEdge & BR,VertexOnGeom & GR) {
[3913]	3132	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/ProjectOnCurve)/
	3133
	3134	void pA=0,pB=0;
	3135	double tA=0,tB=0;
	3136	R2 A=vA,B=vB;
[5120]	3137	BamgVertex * pvA=&vA, * pvB=&vB;
[5460]	3138	if (vA.IndexInTriangle == IsVertexOnVertex){
[5148]	3139	pA=vA.BackgroundVertexHook;
[3913]	3140	}
[5460]	3141	else if (vA.IndexInTriangle == IsVertexOnEdge){
[5148]	3142	pA=vA.BackgroundEdgeHook->be;
	3143	tA=vA.BackgroundEdgeHook->abcisse;
[3913]	3144	}
	3145	else {
[13036]	3146	_error_("ProjectOnCurve On BamgVertex " << BTh.GetId(vA) << " forget call to SetVertexFieldOnBTh");
[3913]	3147	}
	3148
[5460]	3149	if (vB.IndexInTriangle == IsVertexOnVertex){
[5148]	3150	pB=vB.BackgroundVertexHook;
[3913]	3151	}
[5460]	3152	else if(vB.IndexInTriangle == IsVertexOnEdge){
[5148]	3153	pB=vB.BackgroundEdgeHook->be;
	3154	tB=vB.BackgroundEdgeHook->abcisse;
[3913]	3155	}
	3156	else {
[13036]	3157	_error_("ProjectOnCurve On BamgVertex " << BTh.GetId(vB) << " forget call to SetVertexFieldOnBTh");
[3913]	3158	}
	3159	Edge * e = &BhAB;
	3160	if (!pA \|\| !pB \|\| !e){
[13036]	3161	_error_("!pA \|\| !pB \|\| !e");
[3913]	3162	}
	3163	// be carefull the back ground edge e is on same geom edge
	3164	// of the initiale edge def by the 2 vertex A B;
	3165	//check Is a background Mesh;
	3166	if (e<BTh.edges \|\| e>=BTh.edges+BTh.nbe){
[13036]	3167	_error_("e<BTh.edges \|\| e>=BTh.edges+BTh.nbe");
[3913]	3168	}
	3169	// walk on BTh edge
	3170	//not finish ProjectOnCurve with BackGround Mesh);
	3171	// 1 first find a back ground edge contening the vertex A
	3172	// 2 walk n back gound boundary to find the final vertex B
	3173
[5460]	3174	if( vA.IndexInTriangle == IsVertexOnEdge)
[3913]	3175	{ // find the start edge
[5148]	3176	e = vA.BackgroundEdgeHook->be;
[3913]	3177
	3178	}
[5460]	3179	else if (vB.IndexInTriangle == IsVertexOnEdge)
[3913]	3180	{
	3181	theta = 1-theta;
	3182	Exchange(tA,tB);
	3183	Exchange(pA,pB);
	3184	Exchange(pvA,pvB);
	3185	Exchange(A,B);
[5148]	3186	e = vB.BackgroundEdgeHook->be;
[3913]	3187
	3188	}
	3189	else{ // do the search by walking
[13036]	3190	_error_("case not supported yet");
[18064]	3191	}
[3913]	3192
[5148]	3193	// find the direction of walking with direction of edge and pA,PB;
[3913]	3194	R2 AB=B-A;
	3195
	3196	double cosE01AB = (( (R2) (e)[1] - (R2) (e)[0] ) , AB);
	3197	int kkk=0;
[5148]	3198	int direction = (cosE01AB>0) ? 1 : 0;
[3913]	3199
	3200	// double l=0; // length of the edge AB
	3201	double abscisse = -1;
	3202
[5124]	3203	for (int step=0;step<2;step++){
[3913]	3204	// 2 times algo:
	3205	// 1 for computing the length l
	3206	// 2 for find the vertex
	3207	int iii;
[5120]	3208	BamgVertex v0=pvA,v1;
[3913]	3209	Edge neee,eee;
	3210	double lg =0; // length of the curve
	3211	double te0;
	3212	// we suppose take the curve's abcisse
[5148]	3213	for ( eee=e,iii=direction,te0=tA;
[3913]	3214	eee && ((( void) eee) != pB) && (( void) (v1=&((*eee)[iii]))) != pB ;
	3215	neee = eee->adj[iii],iii = 1-neee->Intersection(*eee),eee = neee,v0=v1,te0=1-iii ) {
	3216
	3217	kkk=kkk+1;
[6412]	3218	_assert_(kkk<100);
	3219	_assert_(eee);
[3913]	3220	double lg0 = lg;
	3221	double dp = LengthInterpole(v0->m,v1->m,(R2) v1 - (R2) v0);
	3222	lg += dp;
[5124]	3223	if (step && abscisse <= lg) { // ok we find the geom edge
[3913]	3224	double sss = (abscisse-lg0)/dp;
	3225	double thetab = te0(1-sss)+ sssiii;
[6412]	3226	_assert_(thetab>=0 && thetab<=1);
[3913]	3227	BR = VertexOnEdge(&R,eee,thetab);
	3228	return Gh.ProjectOnCurve(*eee,thetab,R,GR);
[18064]	3229	}
	3230	}
[3913]	3231	// we find the end
	3232	if (v1 != pvB){
	3233	if (( void*) v1 == pB)
	3234	tB = iii;
	3235
	3236	double lg0 = lg;
[6412]	3237	_assert_(eee);
[3913]	3238	v1 = pvB;
	3239	double dp = LengthInterpole(v0->m,v1->m,(R2) v1 - (R2) v0);
	3240	lg += dp;
	3241	abscisse = lg*theta;
	3242	if (abscisse <= lg && abscisse >= lg0 ) // small optimisation we know the lenght because end
	3243	{ // ok we find the geom edge
	3244	double sss = (abscisse-lg0)/dp;
	3245	double thetab = te0(1-sss)+ ssstB;
[6412]	3246	_assert_(thetab>=0 && thetab<=1);
[3913]	3247	BR = VertexOnEdge(&R,eee,thetab);
	3248	return Gh.ProjectOnCurve(*eee,thetab,R,GR);
	3249	}
[18064]	3250	}
[3913]	3251	abscisse = lg*theta;
	3252
[18064]	3253	}
[13036]	3254	_error_("Big bug...");
[3913]	3255	return 0; // just for the compiler
	3256	}
[12365]	3257	/}}}/
[18064]	3258	void Mesh::ReconstructExistingMesh(){/{{{/
	3259	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FillHoleInMesh)/
[3913]	3260
[18064]	3261	/*This routine reconstruct an existing mesh to make it CONVEX:
	3262	* -all the holes are filled
	3263	* -concave boundaries are filled
	3264	* A convex mesh is required for a lot of operations. This is why every mesh
	3265	* goes through this process.
	3266	* This routine also generates mesh properties such as adjencies,...
	3267	*/
[3913]	3268
[18064]	3269	/Intermediary/
	3270	int verbose=0;
[3913]	3271
[18064]	3272	// generation of the integer coordinate
[3913]	3273
[18064]	3274	// find extrema coordinates of vertices pmin,pmax
	3275	long i;
	3276	if(verbose>2) _printf_(" Reconstruct mesh of " << nbv << " vertices\n");
[3913]	3277
[18064]	3278	//initialize orderedvertices
	3279	_assert_(orderedvertices);
	3280	for (i=0;i<nbv;i++) orderedvertices[i]=0;
[3913]	3281
[18064]	3282	//Initialize nbsubdomains
	3283	nbsubdomains =0;
[3913]	3284
[18064]	3285	/* generation of triangles adjacency*/
[3913]	3286
[18064]	3287	//First add existing edges
	3288	long kk =0;
	3289	SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv);
	3290	for (i=0;i<nbe;i++){
	3291	kk=kk+(i==edge4->SortAndAdd(GetId(edges[i][0]),GetId(edges[i][1])));
	3292	}
	3293	if (kk != nbe){
	3294	_error_("There are " << kk-nbe << " double edges in the mesh");
	3295	}
[3913]	3296
[18064]	3297	//Add edges of all triangles in existing mesh
	3298	long* st = new long[nbt*3];
	3299	for (i=0;i<nbt*3;i++) st[i]=-1;
	3300	for (i=0;i<nbt;i++){
	3301	for (int j=0;j<3;j++){
[3913]	3302
[18064]	3303	//Add current triangle edge to edge4
	3304	long k =edge4->SortAndAdd(GetId(triangles[i][VerticesOfTriangularEdge[j][0]]),GetId(triangles[i][VerticesOfTriangularEdge[j][1]]));
[3913]	3305
[18064]	3306	long invisible=triangles[i].Hidden(j);
[3913]	3307
[18064]	3308	//If the edge has not been added to st, add it
	3309	if(st[k]==-1) st[k]=3*i+j;
[3913]	3310
[18064]	3311	//If the edge already exists, add adjacency
	3312	else if(st[k]>=0) {
	3313	_assert_(!triangles[i].TriangleAdj(j));
	3314	_assert_(!triangles[st[k]/3].TriangleAdj((int) (st[k]%3)));
[3913]	3315
[18064]	3316	triangles[i].SetAdj2(j,triangles+st[k]/3,(int)(st[k]%3));
	3317	if (invisible) triangles[i].SetHidden(j);
	3318	if (k<nbe) triangles[i].SetLocked(j);
[3913]	3319
[18064]	3320	//Make st[k] negative so that it will throw an error message if it is found again
	3321	st[k]=-2-st[k];
	3322	}
[3913]	3323
[18064]	3324	//An edge belongs to 2 triangles
	3325	else {
	3326	_error_("The edge (" << GetId(triangles[i][VerticesOfTriangularEdge[j][0]]) << " , " << GetId(triangles[i][VerticesOfTriangularEdge[j][1]]) << ") belongs to more than 2 triangles");
	3327	}
[3913]	3328	}
	3329	}
	3330
[18064]	3331	//Display info if required
	3332	if(verbose>5) {
	3333	_printf_(" info of Mesh:\n");
	3334	_printf_(" - number of vertices = " << nbv << " \n");
	3335	_printf_(" - number of triangles = " << nbt << " \n");
	3336	_printf_(" - number of given edges = " << nbe << " \n");
	3337	_printf_(" - number of all edges = " << edge4->nb() << "\n");
	3338	_printf_(" - Euler number 1 - nb of holes = " << nbt-edge4->nb()+nbv << "\n");
	3339	}
[3913]	3340
[18064]	3341	//check the consistency of edge[].adj and the geometrical required vertex
	3342	long k=0;
	3343	for (i=0;i<edge4->nb();i++){
	3344	if (st[i]>=0){ // edge alone
	3345	if (i<nbe){
	3346	long i0=edge4->i(i);
	3347	orderedvertices[i0] = vertices+i0;
	3348	long i1=edge4->j(i);
	3349	orderedvertices[i1] = vertices+i1;
[3913]	3350	}
[18064]	3351	else {
	3352	k=k+1;
	3353	if (k<10) {
	3354	//print only 10 edges
	3355	_printf_("Lost boundary edges " << i << " : " << edge4->i(i) << " " << edge4->j(i) << "\n");
	3356	}
	3357	else if (k==10){
	3358	_printf_("Other lost boundary edges not shown...\n");
	3359	}
[3913]	3360	}
	3361	}
	3362	}
[18064]	3363	if(k) {
	3364	_error_(k << " boundary edges (from the geometry) are not defined as mesh edges");
	3365	}
[3913]	3366
[18064]	3367	/* mesh generation with boundary points*/
	3368	long nbvb=0;
	3369	for (i=0;i<nbv;i++){
	3370	vertices[i].t=0;
	3371	vertices[i].IndexInTriangle=0;
	3372	if (orderedvertices[i]) orderedvertices[nbvb++]=orderedvertices[i];
	3373	}
[3913]	3374
[18064]	3375	Triangle* savetriangles=triangles;
	3376	long savenbt=nbt;
	3377	long savemaxnbt=maxnbt;
	3378	SubDomain* savesubdomains=subdomains;
	3379	subdomains=0;
[3913]	3380
[18064]	3381	long Nbtriafillhole=2*nbvb;
	3382	Triangle* triafillhole=new Triangle[Nbtriafillhole];
	3383	triangles = triafillhole;
[3913]	3384
[18064]	3385	nbt=2;
	3386	maxnbt= Nbtriafillhole;
[3913]	3387
[18064]	3388	//Find a vertex that is not aligned with vertices 0 and 1
	3389	for (i=2;det(orderedvertices[0]->i,orderedvertices[1]->i,orderedvertices[i]->i)==0;)
	3390	if (++i>=nbvb) {
	3391	_error_("ReconstructExistingMesh: All the vertices are aligned");
	3392	}
	3393	//Move this vertex (i) to the 2d position in orderedvertices
	3394	Exchange(orderedvertices[2], orderedvertices[i]);
[3913]	3395
[18064]	3396	/Reconstruct mesh beginning with 2 triangles/
	3397	BamgVertex * v0=orderedvertices[0], *v1=orderedvertices[1];
[3913]	3398
[18064]	3399	triangles[0](0) = NULL; // Infinite vertex
	3400	triangles[0](1) = v0;
	3401	triangles[0](2) = v1;
[3913]	3402
[18064]	3403	triangles[1](0) = NULL;// Infinite vertex
	3404	triangles[1](2) = v0;
	3405	triangles[1](1) = v1;
	3406	const int e0 = OppositeEdge[0];
	3407	const int e1 = NextEdge[e0];
	3408	const int e2 = PreviousEdge[e0];
	3409	triangles[0].SetAdj2(e0, &triangles[1] ,e0);
	3410	triangles[0].SetAdj2(e1, &triangles[1] ,e2);
	3411	triangles[0].SetAdj2(e2, &triangles[1] ,e1);
[3913]	3412
[18064]	3413	triangles[0].det = -1; // boundary triangles
	3414	triangles[1].det = -1; // boundary triangles
[3913]	3415
[18064]	3416	triangles[0].SetSingleVertexToTriangleConnectivity();
	3417	triangles[1].SetSingleVertexToTriangleConnectivity();
[3913]	3418
[18064]	3419	triangles[0].link=&triangles[1];
	3420	triangles[1].link=&triangles[0];
[3913]	3421
[18064]	3422	if (!quadtree) delete quadtree; //ReInitialise;
	3423	quadtree = new BamgQuadtree(this,0);
	3424	quadtree->Add(*v0);
	3425	quadtree->Add(*v1);
[3913]	3426
[18064]	3427	// vertices are added one by one
	3428	long NbSwap=0;
	3429	for (int icount=2; icount<nbvb; icount++) {
	3430	BamgVertex *vi = orderedvertices[icount];
[22085]	3431	long long det3[3];
[18064]	3432	Triangle *tcvi = TriangleFindFromCoord(vi->i,det3);
	3433	quadtree->Add(*vi);
	3434	AddVertex(*vi,tcvi,det3);
	3435	NbSwap += vi->Optim(1,1);
	3436	}
[3913]	3437
[18064]	3438	//enforce the boundary
	3439	AdjacentTriangle ta(0,0);
	3440	long nbloss = 0,knbe=0;
	3441	for ( i = 0; i < nbe; i++){
	3442	if (st[i] >=0){ //edge alone => on border
	3443	BamgVertex &a=edges[i][0], &b=edges[i][1];
	3444	if (a.t && b.t){
	3445	knbe++;
	3446	if (ForceEdge(a,b,ta)<0) nbloss++;
	3447	}
[3913]	3448	}
	3449	}
[18064]	3450	if(nbloss) {
	3451	_error_("we lost " << nbloss << " existing edges other " << knbe);
	3452	}
[3913]	3453
[18064]	3454	FindSubDomain(1);
	3455	// remove all the hole
	3456	// remove all the good sub domain
	3457	long krm =0;
	3458	for (i=0;i<nbt;i++){
	3459	if (triangles[i].link){ // remove triangles
	3460	krm++;
	3461	for (int j=0;j<3;j++){
	3462	AdjacentTriangle ta = triangles[i].Adj(j);
	3463	Triangle &tta = (Triangle)ta;
	3464	//if edge between remove and not remove
	3465	if(! tta.link){
	3466	// change the link of ta;
	3467	int ja = ta;
	3468	BamgVertex *v0= ta.EdgeVertex(0);
	3469	BamgVertex *v1= ta.EdgeVertex(1);
	3470	long k =edge4->SortAndAdd(v0?GetId(v0):nbv,v1? GetId(v1):nbv);
[3913]	3471
[18064]	3472	_assert_(st[k]>=0);
	3473	tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3));
	3474	ta.SetLock();
	3475	st[k]=-2-st[k];
	3476	}
[3913]	3477	}
	3478	}
	3479	}
[18064]	3480	long NbTfillHoll =0;
	3481	for (i=0;i<nbt;i++){
	3482	if (triangles[i].link) {
	3483	triangles[i]=Triangle((BamgVertex ) NULL,(BamgVertex ) NULL,(BamgVertex *) NULL);
	3484	triangles[i].color=-1;
	3485	}
	3486	else{
	3487	triangles[i].color= savenbt+ NbTfillHoll++;
	3488	}
[3913]	3489	}
[18064]	3490	_assert_(savenbt+NbTfillHoll<=savemaxnbt);
[3913]	3491
[18064]	3492	// copy of the outside triangles in saveMesh
	3493	for (i=0;i<nbt;i++){
	3494	if(triangles[i].color>=0) {
	3495	savetriangles[savenbt]=triangles[i];
	3496	savetriangles[savenbt].link=0;
	3497	savenbt++;
	3498	}
[3913]	3499	}
[18064]	3500	// gestion of the adj
	3501	k =0;
	3502	Triangle * tmax = triangles + nbt;
	3503	for (i=0;i<savenbt;i++) {
	3504	Triangle & ti = savetriangles[i];
	3505	for (int j=0;j<3;j++){
	3506	Triangle * ta = ti.TriangleAdj(j);
	3507	int aa = ti.NuEdgeTriangleAdj(j);
	3508	int lck = ti.Locked(j);
	3509	if (!ta) k++; // bug
	3510	else if ( ta >= triangles && ta < tmax){
	3511	ta= savetriangles + ta->color;
	3512	ti.SetAdj2(j,ta,aa);
	3513	if(lck) ti.SetLocked(j);
	3514	}
[3913]	3515	}
	3516	}
	3517
[18064]	3518	// restore triangles;
	3519	nbt=savenbt;
	3520	maxnbt=savemaxnbt;
	3521	delete [] triangles;
	3522	delete [] subdomains;
	3523	triangles = savetriangles;
	3524	subdomains = savesubdomains;
	3525	if (k) {
	3526	_error_("number of triangles edges alone = " << k);
	3527	}
	3528	FindSubDomain();
[3913]	3529
[18064]	3530	delete edge4;
	3531	delete [] st;
	3532	for (i=0;i<nbv;i++) quadtree->Add(vertices[i]);
[3913]	3533
[18064]	3534	SetVertexFieldOn();
[3913]	3535
[18064]	3536	/Check requirements consistency/
	3537	for (i=0;i<nbe;i++){
	3538	/If the current mesh edge is on Geometry/
	3539	if(edges[i].GeomEdgeHook){
	3540	for(int j=0;j<2;j++){
	3541	/Go through the edges adjacent to current edge (if on the same curve)/
	3542	if (!edges[i].adj[j]){
	3543	/The edge is on Geometry and does not have 2 adjacent edges... (not on a closed curve)/
	3544	/Check that the 2 vertices are on geometry AND required/
	3545	if(!edges[i][j].GeomEdgeHook->IsRequiredVertex()){
	3546	_printf_("ReconstructExistingMesh error message: problem with the edge number " << i+1 << ": [" << GetId(edges[i][0])+1 << " " << GetId(edges[i][1])+1 << "]\n");
	3547	_printf_("This edge is on geometrical edge number " << Gh.GetId(edges[i].GeomEdgeHook)+1 << "\n");
	3548	if (edges[i][j].GeomEdgeHook->OnGeomVertex())
	3549	_printf_("the vertex number " << GetId(edges[i][j])+1 << " of this edge is a geometric BamgVertex number " << Gh.GetId(edges[i][j].GeomEdgeHook->gv)+1 << "\n");
	3550	else if (edges[i][j].GeomEdgeHook->OnGeomEdge())
	3551	_printf_("the vertex number " << GetId(edges[i][j])+1 << " of this edge is a geometric Edge number " << Gh.GetId(edges[i][j].GeomEdgeHook->ge)+1 << "\n");
	3552	else
	3553	_printf_("Its pointer is " << edges[i][j].GeomEdgeHook << "\n");
[3913]	3554
[18064]	3555	_printf_("This edge is on geometry and has no adjacent edge (open curve) and one of the tip is not required\n");
	3556	_error_("See above (might be cryptic...)");
	3557	}
[3913]	3558	}
	3559	}
	3560	}
	3561	}
	3562	}
[18064]	3563	/}}}/
	3564	void Mesh::TrianglesRenumberBySubDomain(bool justcompress){/{{{/
[3913]	3565	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingTheTriangleBySubDomain)/
	3566
	3567	long *renu= new long[nbt];
[17507]	3568	Triangle t0,t,*te=triangles+nbt;
	3569	long k=0,it,i,j;
[3913]	3570
	3571	for ( it=0;it<nbt;it++)
	3572	renu[it]=-1; // outside triangle
[5148]	3573	for ( i=0;i<nbsubdomains;i++)
[3913]	3574	{
	3575	t=t0=subdomains[i].head;
	3576	if (!t0){ // not empty sub domain
[13036]	3577	_error_("!t0");
[3913]	3578	}
	3579	do {
[5149]	3580	long kt = GetId(t);
[3913]	3581	if (kt<0 \|\| kt >= nbt ){
[13036]	3582	_error_("kt<0 \|\| kt >= nbt");
[3913]	3583	}
	3584	if (renu[kt]!=-1){
[13036]	3585	_error_("renu[kt]!=-1");
[3913]	3586	}
	3587	renu[kt]=k++;
	3588	}
	3589	while (t0 != (t=t->link));
	3590	}
	3591	// take is same numbering if possible
	3592	if(justcompress)
	3593	for ( k=0,it=0;it<nbt;it++)
	3594	if(renu[it] >=0 )
	3595	renu[it]=k++;
	3596
	3597	// put the outside triangles at the end
	3598	for ( it=0;it<nbt;it++){
	3599	if (renu[it]==-1) renu[it]=k++;
	3600	}
	3601	if (k != nbt){
[13036]	3602	_error_("k != nbt");
[3913]	3603	}
	3604	// do the change on all the pointeur
	3605	for ( it=0;it<nbt;it++)
[5124]	3606	triangles[it].Renumbering(triangles,te,renu);
[3913]	3607
[5148]	3608	for ( i=0;i<nbsubdomains;i++)
[5149]	3609	subdomains[i].head=triangles+renu[GetId(subdomains[i].head)];
[3913]	3610
	3611	// move the Triangles without a copy of the array
	3612	// be carefull not trivial code
	3613	for ( it=0;it<nbt;it++) // for all sub cycles of the permutation renu
	3614	if (renu[it] >= 0) // a new sub cycle
	3615	{
	3616	i=it;
	3617	Triangle ti=triangles[i],tj;
	3618	while ( (j=renu[i]) >= 0)
	3619	{ // i is old, and j is new
	3620	renu[i] = -1; // mark
	3621	tj = triangles[j]; // save new
	3622	triangles[j]= ti; // new <- old
	3623	i=j; // next
	3624	ti = tj;
	3625	}
	3626	}
	3627	delete [] renu;
	3628
	3629	}
[12365]	3630	/}}}/
[18064]	3631	void Mesh::SetIntCoor(const char * strfrom) {/{{{/
	3632	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SetIntCoor)/
[3913]	3633
[18064]	3634	/Set integer coordinate for existing vertices/
[3913]	3635
[18064]	3636	//Get extrema coordinates of the existing vertices
	3637	pmin = vertices[0].r;
	3638	pmax = vertices[0].r;
	3639	long i;
	3640	for (i=0;i<nbv;i++) {
	3641	pmin.x = Min(pmin.x,vertices[i].r.x);
	3642	pmin.y = Min(pmin.y,vertices[i].r.y);
	3643	pmax.x = Max(pmax.x,vertices[i].r.x);
	3644	pmax.y = Max(pmax.y,vertices[i].r.y);
	3645	}
	3646	R2 DD = (pmax-pmin)*0.05;
	3647	pmin = pmin-DD;
	3648	pmax = pmax+DD;
[3913]	3649
[18064]	3650	//Compute coefIcoor
[21759]	3651	int MaxICoord = 1073741823; //2^30 - 1 = =111...111 (29 times one)
	3652	coefIcoor= (MaxICoord)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
[18064]	3653	if (coefIcoor<=0){
	3654	_error_("coefIcoor should be positive, a problem in the geometry is likely");
	3655	}
[3913]	3656
[18064]	3657	// generation of integer coord
	3658	for (i=0;i<nbv;i++) {
	3659	vertices[i].i = R2ToI2(vertices[i].r);
	3660	}
[3913]	3661
[18064]	3662	// computation of the det
	3663	int number_of_errors=0;
	3664	for (i=0;i<nbt;i++) {
[19293]	3665	BamgVertex* v0 = triangles[i](0);
	3666	BamgVertex* v1 = triangles[i](1);
	3667	BamgVertex* v2 = triangles[i](2);
[3913]	3668
[18064]	3669	//If this is not a boundary triangle
[19293]	3670	if (v0 && v1 && v2){
[5010]	3671
[18064]	3672	/Compute determinant/
[19293]	3673	triangles[i].det= det(v0->GetIntegerCoordinates(),v1->GetIntegerCoordinates(),v2->GetIntegerCoordinates());
[5010]	3674
[18064]	3675	/Check that determinant is positive/
	3676	if (triangles[i].det <=0){
[5010]	3677
[18064]	3678	/increase number_of_errors and print error only for the first 20 triangles/
	3679	number_of_errors++;
	3680	if (number_of_errors<20){
	3681	_printf_("Area of Triangle " << i+1 << " < 0 (det=" << triangles[i].det << ")\n");
	3682	}
[5010]	3683	}
[3913]	3684	}
[18064]	3685
	3686	//else, set as -1
	3687	else triangles[i].det=-1;
[3913]	3688	}
	3689
[18064]	3690	if (number_of_errors) _error_("Fatal error: some triangles have negative areas, see above");
[3913]	3691	}
[18064]	3692	/}}}/
	3693	void Mesh::SmoothingVertex(int nbiter,double omega ) { /{{{/
	3694	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SmoothingVertex)/
[5010]	3695
[18064]	3696	long int verbose=0;
	3697	// if quatree exist remove it end reconstruct
	3698	if (quadtree) delete quadtree;
	3699	quadtree=0;
	3700	CreateSingleVertexToTriangleConnectivity();
	3701	Triangle vide; // a triangle to mark the boundary vertex
	3702	Triangle ** tstart= new Triangle* [nbv];
	3703	long i,j,k;
	3704	// attention si Background == Triangle alors on ne peut pas utiliser la rechech rapide
	3705	if ( this == & BTh)
[3913]	3706	for ( i=0;i<nbv;i++)
[18064]	3707	tstart[i]=vertices[i].t;
	3708	else
	3709	for ( i=0;i<nbv;i++)
	3710	tstart[i]=0;
	3711	for ( j=0;j<NbVerticesOnGeomVertex;j++ )
	3712	tstart[ GetId(VerticesOnGeomVertex[j].meshvertex)]=&vide;
	3713	for ( k=0;k<NbVerticesOnGeomEdge;k++ )
	3714	tstart[ GetId(VerticesOnGeomEdge[k].meshvertex)]=&vide;
	3715	if(verbose>2) _printf_(" SmoothingVertex: nb Iteration = " << nbiter << ", Omega=" << omega << "\n");
	3716	for (k=0;k<nbiter;k++)
	3717	{
	3718	long i,NbSwap =0;
	3719	double delta =0;
	3720	for ( i=0;i<nbv;i++)
	3721	if (tstart[i] != &vide) // not a boundary vertex
	3722	delta=Max(delta,vertices[i].Smoothing(*this,BTh,tstart[i],omega));
[21759]	3723	for ( i=0;i<nbv;i++)
	3724	if (tstart[i] != &vide) // not a boundary vertex
	3725	NbSwap += vertices[i].Optim(1);
[18064]	3726	if (verbose>3) _printf_(" move max = " << pow(delta,0.5) << ", iteration = " << k << ", nb of swap = " << NbSwap << "\n");
	3727	}
[3913]	3728
[18064]	3729	delete [] tstart;
	3730	if (quadtree) quadtree= new BamgQuadtree(this);
	3731	}
	3732	/}}}/
	3733	void Mesh::SmoothMetric(double raisonmax) { /{{{/
	3734	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/SmoothMetric)/
[3913]	3735
[18064]	3736	long int verbose=0;
[3913]	3737
[18064]	3738	if(raisonmax<1.1) return;
	3739	if(verbose > 1) _printf_(" Mesh::SmoothMetric raisonmax = " << raisonmax << "\n");
	3740	CreateSingleVertexToTriangleConnectivity();
	3741	long i,j,kch,kk,ip;
	3742	long *first_np_or_next_t0 = new long[nbv];
	3743	long *first_np_or_next_t1 = new long[nbv];
	3744	long Head0 =0,Head1=-1;
	3745	double logseuil= log(raisonmax);
[3913]	3746
[18064]	3747	for(i=0;i<nbv-1;i++)
	3748	first_np_or_next_t0[i]=i+1;
	3749	first_np_or_next_t0[nbv-1]=-1;// end;
	3750	for(i=0;i<nbv;i++)
	3751	first_np_or_next_t1[i]=-1;
	3752	kk=0;
	3753	while(Head0>=0&& kk++<100){
	3754	kch=0;
	3755	for(i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1) {
	3756	// pour tous les triangles autour du sommet s
	3757	Triangle* t= vertices[i].t;
	3758	if (!t){
	3759	_error_("!t");
[3913]	3760	}
[18064]	3761	BamgVertex & vi = vertices[i];
	3762	AdjacentTriangle ta(t,EdgesVertexTriangle[vertices[i].IndexInTriangle][0]);
	3763	BamgVertex *pvj0 = ta.EdgeVertex(0);
	3764	while (1) {
	3765	ta=Previous(Adj(ta));
	3766	if (vertices+i != ta.EdgeVertex(1)){
	3767	_error_("vertices+i != ta.EdgeVertex(1)");
[3913]	3768	}
[20621]	3769	BamgVertex *pvj = (ta.EdgeVertex(0));
	3770	BamgVertex & vj = *pvj;
	3771	if(pvj){
[18064]	3772	j= &vj-vertices;
	3773	if (j<0 \|\| j >= nbv){
	3774	_error_("j<0 \|\| j >= nbv");
	3775	}
	3776	R2 Aij = (R2) vj - (R2) vi;
	3777	double ll = Norme2(Aij);
	3778	if (0) {
[21759]	3779	double hi = ll/vi.m.Length(Aij.x,Aij.y);
	3780	double hj = ll/vj.m.Length(Aij.x,Aij.y);
[18064]	3781	if(hi < hj)
	3782	{
	3783	double dh=(hj-hi)/ll;
	3784	if (dh>logseuil) {
	3785	vj.m.IntersectWith(vi.m/(1 +logseuil*ll/hi));
	3786	if(first_np_or_next_t1[j]<0)
	3787	kch++,first_np_or_next_t1[j]=Head1,Head1=j;
	3788	}
	3789	}
	3790	}
	3791	else
[3913]	3792	{
[21759]	3793	double li = vi.m.Length(Aij.x,Aij.y);
[18064]	3794	if( vj.m.IntersectWith(vi.m/(1 +logseuil*li)) )
	3795	if(first_np_or_next_t1[j]<0) // if the metrix change
	3796	kch++,first_np_or_next_t1[j]=Head1,Head1=j;
[3913]	3797	}
[18064]	3798	}
	3799	if ( &vj == pvj0 ) break;
[3913]	3800	}
	3801	}
[18064]	3802	Head0 = Head1;
	3803	Head1 = -1;
	3804	Exchange(first_np_or_next_t0,first_np_or_next_t1);
[3913]	3805	}
[18064]	3806	if(verbose>2) _printf_(" number of iterations = " << kch << "\n");
	3807	delete [] first_np_or_next_t0;
	3808	delete [] first_np_or_next_t1;
[3913]	3809	}
[18064]	3810	/}}}/
	3811	long Mesh::SplitInternalEdgeWithBorderVertices(){/{{{/
	3812	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SplitInternalEdgeWithBorderVertices)/
[3913]	3813
[18064]	3814	long NbSplitEdge=0;
	3815	SetVertexFieldOn();
	3816	long it;
	3817	long nbvold=nbv;
	3818	long int verbose=2;
	3819	for (it=0;it<nbt;it++){
	3820	Triangle &t=triangles[it];
	3821	if (t.link)
	3822	for (int j=0;j<3;j++)
	3823	if(!t.Locked(j) && !t.Hidden(j)){
[20621]	3824
	3825	Triangle *ptt = t.TriangleAdj(j);
	3826	Triangle &tt = *ptt;
	3827
	3828	if( ptt && tt.link && it < GetId(tt))
[18064]	3829	{ // an internal edge
	3830	BamgVertex &v0 = t[VerticesOfTriangularEdge[j][0]];
	3831	BamgVertex &v1 = t[VerticesOfTriangularEdge[j][1]];
	3832	if (v0.GeomEdgeHook && v1.GeomEdgeHook){
	3833	R2 P= ((R2) v0 + (R2) v1)*0.5;
	3834	if ( nbv<maxnbv) {
	3835	vertices[nbv].r = P;
	3836	vertices[nbv++].m = Metric(0.5,v0.m,0.5,v1.m);
	3837	vertices[nbv].ReferenceNumber=0;
	3838	}
	3839	NbSplitEdge++;
[3913]	3840	}
[18064]	3841	}
	3842	}
	3843	}
	3844	CreateSingleVertexToTriangleConnectivity();
	3845	if (nbvold!=nbv){
	3846	long iv = nbvold;
	3847	long NbSwap = 0;
[22085]	3848	long long det3[3];
[18064]	3849	for (int i=nbvold;i<nbv;i++) {// for all the new point
	3850	BamgVertex & vi = vertices[i];
	3851	vi.i = R2ToI2(vi.r);
	3852	vi.r = I2ToR2(vi.i);
[3913]	3853
[18064]	3854	// a good new point
	3855	vi.ReferenceNumber=0;
	3856	Triangle *tcvi = TriangleFindFromCoord(vi.i,det3);
	3857	if (tcvi && !tcvi->link) {
	3858	_printf_("problem inserting point in SplitInternalEdgeWithBorderVertices (tcvj && !tcvj->link)\n");
	3859	}
[3913]	3860
[18064]	3861	quadtree->Add(vi);
	3862	if (!tcvi \|\| tcvi->det<0){// internal
	3863	_error_("!tcvi \|\| tcvi->det < 0");
	3864	}
	3865	AddVertex(vi,tcvi,det3);
	3866	NbSwap += vi.Optim(1);
	3867	iv++;
[3913]	3868	}
[18064]	3869	if (verbose>3) {
	3870	_printf_(" number of points: " << iv << "\n");
	3871	_printf_(" number of swap to split internal edges with border vertices: " << NbSwap << "\n");
	3872	nbv = iv;
	3873	}
[12876]	3874	}
[18064]	3875	if (NbSplitEdge>nbv-nbvold) _printf_("WARNING: not enough vertices to split all internal edges, we lost " << NbSplitEdge - ( nbv-nbvold) << " edges...\n");
	3876	if (verbose>2) _printf_("SplitInternalEdgeWithBorderVertices: Number of splited edge " << NbSplitEdge << "\n");
	3877
	3878	return NbSplitEdge;
[3913]	3879	}
[18064]	3880	/}}}/
	3881	I2 Mesh::R2ToI2(const R2 & P) const {/{{{/
[21759]	3882	return I2( (int) (coefIcoor(P.x-pmin.x)),(int) (coefIcoor(P.y-pmin.y)) );
[18064]	3883	}
	3884	/}}}/
	3885	R2 Mesh::I2ToR2(const I2 & P) const {/{{{/
	3886	return R2( (double) P.x/coefIcoor+pmin.x, (double) P.y/coefIcoor+pmin.y);
	3887	}
	3888	/}}}/
[22085]	3889	Triangle * Mesh::TriangleFindFromCoord(const I2 & B,long long det3[3], Triangle tstart){/{{{*/
[18064]	3890	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FindTriangleContening)/
[3913]	3891
[18064]	3892	Triangle * t=0;
	3893	int j,jp,jn,jj;
	3894	int counter;
[5124]	3895
[18064]	3896	/Get starting triangle. Take tsart if provided/
	3897	if (tstart) t=tstart;
[5124]	3898
[18064]	3899	/Else find the closest Triangle using the quadtree/
	3900	else {
[5124]	3901
[18064]	3902	/Check that the quadtree does exist/
	3903	if (!quadtree) _error_("no starting triangle provided and no quadtree available");
[5124]	3904
[18064]	3905	/Call NearestVertex/
	3906	BamgVertex *a = quadtree->NearestVertex(B.x,B.y) ;
[5124]	3907
[18064]	3908	/Check output (Vertex a)/
	3909	if (!a) _error_("problem while trying to find nearest vertex from a given point. No output found");
	3910	if (!a->t) _error_("no triangle is associated to vertex number " << GetId(a)+1 << " (orphan?)");
	3911	_assert_(a>=vertices && a<vertices+nbv);
[5124]	3912
[18064]	3913	/Get starting triangle/
	3914	t = a->t;
	3915	_assert_(t>=triangles && t<triangles+nbt);
	3916	}
[5124]	3917
[22085]	3918	long long detop ;
[5124]	3919
[18064]	3920	/initialize number of test triangle/
	3921	counter=0;
[5124]	3922
[18064]	3923	/The initial triangle might be outside/
	3924	while (t->det < 0){
[5124]	3925
[18064]	3926	/Get a real vertex from this triangle (k0)/
	3927	int k0=(t)(0)?(((t)(1)?((*t)(2)?-1:2):1)):0;
	3928	_assert_(k0>=0);// k0 the NULL vertex
	3929	int k1=NextVertex[k0],k2=PreviousVertex[k0];
	3930	det3[k0]=det(B,(t)[k1],(t)[k2]);
	3931	det3[k1]=det3[k2]=-1;
	3932	if (det3[k0] > 0) // outside B
	3933	return t;
	3934	t = t->TriangleAdj(OppositeEdge[k0]);
	3935	counter++;
	3936	_assert_(counter<2);
	3937	}
[5124]	3938
[18064]	3939	jj=0;
	3940	detop = det((t)(VerticesOfTriangularEdge[jj][0]),(t)(VerticesOfTriangularEdge[jj][1]),B);
[5124]	3941
[18064]	3942	while(t->det>0){
[5124]	3943
[18064]	3944	/Increase counter/
	3945	if (++counter>=10000) _error_("Maximum number of iteration reached (threshold = " << counter << ").");
[5124]	3946
[18064]	3947	j= OppositeVertex[jj];
	3948	det3[j] = detop; //det(b,s1,*s2);
	3949	jn = NextVertex[j];
	3950	jp = PreviousVertex[j];
	3951	det3[jp]= det((t)(j),(t)(jn),B);
	3952	det3[jn] = t->det-det3[j] -det3[jp];
[5124]	3953
[18064]	3954	// count the number k of det3 <0
	3955	int k=0,ii[3];
	3956	if (det3[0] < 0 ) ii[k++]=0;
	3957	if (det3[1] < 0 ) ii[k++]=1;
	3958	if (det3[2] < 0 ) ii[k++]=2;
	3959	// 0 => ok
	3960	// 1 => go in way 1
	3961	// 2 => two way go in way 1 or 2 randomly
[5124]	3962
[18064]	3963	if (k==0) break;
[22085]	3964	if (k==2 && this->RandomNumber(1)) Exchange(ii[0],ii[1]);
[18064]	3965	_assert_(k<3);
	3966	AdjacentTriangle t1 = t->Adj(jj=ii[0]);
	3967	if ((t1.det() < 0 ) && (k == 2))
	3968	t1 = t->Adj(jj=ii[1]);
	3969	t=t1;
	3970	j=t1;// for optimisation we now the -det[OppositeVertex[j]];
	3971	detop = -det3[OppositeVertex[jj]];
	3972	jj = j;
	3973	}
[5124]	3974
[18064]	3975	if (t->det<0) // outside triangle
	3976	det3[0]=det3[1]=det3[2]=-1,det3[OppositeVertex[jj]]=detop;
	3977	return t;
[5124]	3978	}
[18064]	3979	/}}}/
	3980	void Mesh::TriangleIntNumbering(long* renumbering){/{{{/
[5124]	3981
[18064]	3982	long num=0;
	3983	for (int i=0;i<nbt;i++){
	3984	if (triangles[i].det>0) renumbering[i]=num++;
	3985	else renumbering[i]=-1;
	3986	}
	3987	return;
[5124]	3988	}
[18064]	3989	/}}}/
	3990	long Mesh::TriangleReferenceList(long* reft) const {/{{{/
	3991	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ConsRefTriangle)/
[5124]	3992
[18064]	3993	Triangle t0,t;
	3994	long k=0, num;
[5124]	3995
[18064]	3996	//initialize all triangles as -1 (outside)
	3997	for (int it=0;it<nbt;it++) reft[it]=-1;
[5124]	3998
[18064]	3999	//loop over all subdomains
	4000	for (int i=0;i<nbsubdomains;i++){
[5124]	4001
[18064]	4002	//first triangle of the subdomain i
	4003	t=t0=subdomains[i].head;
[5124]	4004
[18064]	4005	//check that the subdomain is not empty
	4006	if (!t0){ _error_("At least one subdomain is empty");}
[5124]	4007
[18064]	4008	//loop
	4009	do{
	4010	k++;
[5124]	4011
[18064]	4012	//get current triangle number
	4013	num = GetId(t);
[5124]	4014
[18064]	4015	//check that num is in [0 nbt[
	4016	_assert_(num>=0 && num<nbt);
[5124]	4017
[18064]	4018	//reft of this triangle is the subdomain number
	4019	reft[num]=i;
[5124]	4020
[18064]	4021	} while (t0 != (t=t->link));
	4022	//stop when all triangles of subdomains have been tagged
[5124]	4023
[18064]	4024	}
	4025	return k;
[5124]	4026	}
[18064]	4027	/}}}/
	4028	void Mesh::Triangulate(double* x,double* y,int nods){/{{{/
[10205]	4029
[18064]	4030	int verbose=0;
	4031	int i;
	4032	Metric M1(1);
[10205]	4033
[18064]	4034	/Initialize mesh/
	4035	Init(nods);//this resets nbv to 0
	4036	nbv=nods;
[10205]	4037
[18064]	4038	//Vertices
	4039	if(verbose) _printf_("Reading vertices (" << nbv << ")\n");
	4040	for(i=0;i<nbv;i++){
	4041	vertices[i].r.x=x[i];
	4042	vertices[i].r.y=y[i];
	4043	vertices[i].ReferenceNumber=1;
	4044	vertices[i].m=M1;
	4045	vertices[i].color=0;
	4046	}
	4047	maxnbt=2*maxnbv-2; // for filling The Holes and quadrilaterals
	4048
	4049	/Insert Vertices/
[22085]	4050	Insert();
[10205]	4051	}
[18064]	4052	/}}}/
	4053	void Mesh::TriangulateFromGeom0(BamgOpts* bamgopts){/{{{/
[5120]	4054	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/GeomToTriangles0)/
	4055	/Generate mesh from geometry/
	4056
[5448]	4057	/Intermediaries/
[13761]	4058	int i,k;
[5448]	4059	int nbcurves = 0;
	4060	int NbNewPoints,NbEdgeCurve;
	4061	double lcurve,lstep,s;
	4062	const int MaxSubEdge = 10;
[5120]	4063
[10205]	4064	R2 AB;
[5573]	4065	GeomVertex a, b;
[10205]	4066	BamgVertex va,vb;
[5573]	4067	GeomEdge *e;
[5120]	4068
[5448]	4069	// add a ref to GH to make sure that it is not destroyed by mistake
	4070	Gh.NbRef++;
[5120]	4071
	4072	/Get options/
	4073	int verbose=bamgopts->verbose;
	4074
	4075	//build background mesh flag (1 if background, else 0)
[5448]	4076	bool background=(&BTh != this);
[5120]	4077
[5448]	4078	/Build VerticesOnGeomVertex/
	4079
	4080	//Compute the number of geometrical vertices that we are going to use to mesh
[5120]	4081	for (i=0;i<Gh.nbv;i++){
[5124]	4082	if (Gh[i].Required()) NbVerticesOnGeomVertex++;
[5120]	4083	}
[5448]	4084	//allocate
[5120]	4085	VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];
[13036]	4086	if(NbVerticesOnGeomVertex >= maxnbv) _error_("too many vertices on geometry: " << NbVerticesOnGeomVertex << " >= " << maxnbv);
[6412]	4087	_assert_(nbv==0);
[5448]	4088	//Build VerticesOnGeomVertex
[5120]	4089	for (i=0;i<Gh.nbv;i++){
[5124]	4090	/* Add vertex only if required*/
	4091	if (Gh[i].Required()) {//Gh vertices Required
[5120]	4092
[5448]	4093	//Add the vertex
[6412]	4094	_assert_(nbv<maxnbv);
[5448]	4095	vertices[nbv]=Gh[i];
[13622]	4096
[5120]	4097	//Add pointer from geometry (Gh) to vertex from mesh (Th)
[5460]	4098	Gh[i].MeshVertexHook=vertices+nbv;
[5120]	4099
	4100	//Build VerticesOnGeomVertex for current point
[5460]	4101	VerticesOnGeomVertex[nbv]=VertexOnGeom(*Gh[i].MeshVertexHook,Gh[i]);
[5120]	4102
	4103	//nbv increment
	4104	nbv++;
	4105	}
	4106	}
	4107
[5448]	4108	/Build VerticesOnGeomEdge/
[5120]	4109
[5448]	4110	//check that edges is still empty (Init)
[6412]	4111	_assert_(!edges);
[5448]	4112
[5120]	4113	/* Now we are going to create the first edges corresponding
	4114	* to the one present in the geometry provided.
	4115	* We proceed in 2 steps
[5531]	4116	* -step 0: we count all the edges
	4117	* we allocate the number of edges at the end of step 0
	4118	* -step 1: the edges are created */
[5120]	4119	for (int step=0;step<2;step++){
	4120
	4121	//initialize number of edges and number of edges max
	4122	long nbex=0;
	4123	nbe=0;
	4124	long NbVerticesOnGeomEdge0=NbVerticesOnGeomEdge;
	4125	Gh.UnMarkEdges();
[5148]	4126	nbcurves=0;
[5120]	4127
	4128	//go through the edges of the geometry
[5448]	4129	for (i=0;i<Gh.nbe;i++){
[5120]	4130
	4131	//ei = current Geometrical edge
[5573]	4132	GeomEdge &ei=Gh.edges[i];
[5120]	4133
[5448]	4134	//loop over the two vertices of the edge ei
[5120]	4135	for(int j=0;j<2;j++) {
	4136
[5531]	4137	/Take only required vertices (corner->beginning of a new curve)/
[5120]	4138	if (!ei.Mark() && ei[j].Required()){
	4139
	4140	long nbvend=0;
	4141	Edge* PreviousNewEdge=NULL;
	4142	lstep = -1;
	4143
[5448]	4144	/If Edge is required (do that only once for the 2 vertices)/
[5531]	4145	if(ei.Required()){
	4146	if (j==0){
	4147	//do not create internal points if required (take it as is)
	4148	if(step==0) nbe++;
	4149	else{
	4150	e=&ei;
	4151	a=ei(0);
	4152	b=ei(1);
[5120]	4153
[5531]	4154	//check that edges has been allocated
[6412]	4155	_assert_(edges);
[5531]	4156	edges[nbe].v[0]=a->MeshVertexHook;
	4157	edges[nbe].v[1]=b->MeshVertexHook;;
	4158	edges[nbe].ReferenceNumber = e->ReferenceNumber;
[5573]	4159	edges[nbe].GeomEdgeHook = e;
[5531]	4160	edges[nbe].adj[0] = 0;
	4161	edges[nbe].adj[1] = 0;
	4162	nbe++;
	4163	}
[5120]	4164	}
	4165	}
	4166
[5448]	4167	/If Edge is not required: we are on a curve/
[5120]	4168	else {
	4169	for (int kstep=0;kstep<=step;kstep++){
[5531]	4170	//kstep=0, compute number of edges (discretize curve)
	4171	//kstep=1 create the points and edge
[5120]	4172	PreviousNewEdge=0;
	4173	NbNewPoints=0;
	4174	NbEdgeCurve=0;
[13036]	4175	if (nbvend>=maxnbv) _error_("maximum number of vertices too low! Check the domain outline or increase maxnbv");
[5120]	4176	lcurve =0;
[5531]	4177	s = lstep; //-1 initially, then length of each sub edge
[5120]	4178
[5448]	4179	/reminder: i = edge number, j=[0;1] vertex index in edge/
	4180	k=j; // k = vertex index in edge (0 or 1)
[5120]	4181	e=&ei; // e = reference of current edge
[5124]	4182	a=ei(k); // a = pointer toward the kth vertex of the current edge
[5460]	4183	va = a->MeshVertexHook; // va = pointer toward mesh vertex associated
[5120]	4184	e->SetMark(); // Mark edge
	4185
[5531]	4186	/Loop until we reach the end of the curve/
[5120]	4187	for(;;){
[5448]	4188	k = 1-k; // other vertx index of the curve
	4189	b = (*e)(k); // b = pointer toward the other vertex of the current edge
[5120]	4190	AB= b->r - a->r; // AB = vector of the current edge
	4191	Metric MA = background ? BTh.MetricAt(a->r) :a->m ; //Get metric associated to A
[5531]	4192	Metric MB = background ? BTh.MetricAt(b->r) :b->m ; //Get metric associated to B
[21759]	4193	double ledge = (MA.Length(AB.x,AB.y) + MB.Length(AB.x,AB.y))/2; //Get edge length in metric
[5120]	4194
[5448]	4195	/* We are now creating the mesh edges from the geometrical edge selected above.
	4196	* The edge will be divided according to the metric previously computed and cannot
	4197	* be divided more than 10 times (MaxSubEdge). */
[5120]	4198
	4199	//By default, there is only one subedge that is the geometrical edge itself
	4200	int NbSubEdge = 1;
	4201
	4202	//initialize lSubEdge, holding the length of each subedge (cannot be higher than 10)
	4203	double lSubEdge[MaxSubEdge];
	4204
	4205	//Build Subedges according to the edge length
[5448]	4206	if (ledge < 1.5){
	4207	//if ledge < 1.5 (between one and 2), take the edge as is
	4208	lSubEdge[0] = ledge;
	4209	}
[5120]	4210	//else, divide the edge
	4211	else {
[5531]	4212	//compute number of subedges (division of the edge), Maximum is 10
[5120]	4213	NbSubEdge = Min( MaxSubEdge, (int) (ledge +0.5));
[5531]	4214	/*Now, we are going to divide the edge according to the metric.
	4215	* Get segment by sement along the edge.
	4216	* Build lSubEdge, which holds the distance between the first vertex
	4217	* of the edge and the next point on the edge according to the
	4218	* discretization (each SubEdge is AB)*/
[5120]	4219	R2 A,B;
	4220	A=a->r;
	4221	Metric MAs=MA,MBs;
	4222	ledge=0;
	4223	double x =0, xstep= 1./NbSubEdge;
	4224	for (int kk=0; kk < NbSubEdge; kk++,A=B,MAs=MBs ) {
	4225	x += xstep;
	4226	B = e->F(k ? x : 1-x);
[5489]	4227	MBs= background ? BTh.MetricAt(B) : Metric(1-x,MA,x,MB);
[5120]	4228	AB = A-B;
[21759]	4229	lSubEdge[kk]=(ledge+=(MAs.Length(AB.x,AB.y)+MBs.Length(AB.x,AB.y))/2);
[5120]	4230	}
	4231	}
	4232
[5489]	4233	double lcurveb = lcurve+ledge;
[5120]	4234
[5489]	4235	/Now, create corresponding points/
[5531]	4236	while(s>=lcurve && s<=lcurveb && nbv<nbvend){
[5120]	4237
[5581]	4238	/*Schematic of current curve
	4239	*
	4240	* a vb b // vertex
	4241	* 0 ll0 ll1 ledge // length from a
	4242	* + --- + - ... - + --S-- + --- + - ... - + // where is S
	4243	* 0 kk0 kk1 NbSubEdge // Sub edge index
	4244	*
	4245	*/
	4246
[5120]	4247	double ss = s-lcurve;
	4248
	4249	/Find the SubEdge containing ss using Dichotomy/
	4250	int kk0=-1,kk1=NbSubEdge-1,kkk;
	4251	double ll0=0,ll1=ledge,llk;
	4252	while (kk1-kk0>1){
	4253	if (ss < (llk=lSubEdge[kkk=(kk0+kk1)/2]))
	4254	kk1=kkk,ll1=llk;
	4255	else
	4256	kk0=kkk,ll0=llk;
	4257	}
[6412]	4258	_assert_(kk1!=kk0);
[5120]	4259
[5581]	4260	/Curvilinear coordinate in [0 1] of ss in current edge/
	4261	// WARNING: This is what we would do
	4262	// ssa = (ss-ll0)/(ll1-ll0);
	4263	// aa = (kk0+ssa)/NbSubEdge
	4264	// This is what Bamg does:
	4265	double sbb = (ss-ll0)/(ll1-ll0);
	4266	/Curvilinear coordinate in [0 1] of ss in current curve/
	4267	double bb = (kk1+sbb)/NbSubEdge;
	4268	double aa = 1-bb;
[5120]	4269
	4270	// new vertex on edge
	4271	vb = &vertices[nbv++];
	4272	vb->m = Metric(aa,a->m,bb,b->m);
[5148]	4273	vb->ReferenceNumber = e->ReferenceNumber;
[5120]	4274	double abcisse = k ? bb : aa;
[5581]	4275	vb->r = e->F(abcisse);
[5120]	4276	VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(vb,e,abcisse);
	4277
	4278	// to take into account the direction of the edge
	4279	s += lstep;
	4280	edges[nbe].v[0]=va;
	4281	edges[nbe].v[1]=vb;
[5148]	4282	edges[nbe].ReferenceNumber =e->ReferenceNumber;
[5573]	4283	edges[nbe].GeomEdgeHook = e;
[5120]	4284	edges[nbe].adj[0] = PreviousNewEdge;
	4285	if(PreviousNewEdge) PreviousNewEdge->adj[1]=&edges[nbe];
	4286	PreviousNewEdge=edges+nbe;
	4287	nbe++;
	4288	va = vb;
	4289	}
[13622]	4290
[5531]	4291	/We just added one edge to the curve: Go to the next one/
[5120]	4292	lcurve = lcurveb;
	4293	e->SetMark();
	4294	a=b;
[5531]	4295
	4296	/If b is required, we are on a new curve->break/
	4297	if (b->Required()) break;
[5120]	4298	int kprev=k;
[5401]	4299	k = e->AdjVertexIndex[kprev];// next vertices
[5120]	4300	e = e->Adj[kprev];
[6412]	4301	_assert_(e);
[5120]	4302	}// for(;;)
[5460]	4303	vb = b->MeshVertexHook;
[5531]	4304
	4305	/Number of edges in the last disretized curve/
[5120]	4306	NbEdgeCurve = Max((long) (lcurve +0.5), (long) 1);
[5531]	4307	/Number of internal vertices in the last disretized curve/
[5120]	4308	NbNewPoints = NbEdgeCurve-1;
	4309	if(!kstep){
	4310	NbVerticesOnGeomEdge0 += NbNewPoints;
[5148]	4311	nbcurves++;
[5120]	4312	}
	4313	nbvend=nbv+NbNewPoints;
[5531]	4314	lstep = lcurve / NbEdgeCurve; //approximately one
[5120]	4315	}// end of curve --
	4316	if (edges) { // last edges of the curves
	4317	edges[nbe].v[0]=va;
	4318	edges[nbe].v[1]=vb;
[5148]	4319	edges[nbe].ReferenceNumber = e->ReferenceNumber;
[5573]	4320	edges[nbe].GeomEdgeHook = e;
[5120]	4321	edges[nbe].adj[0] = PreviousNewEdge;
	4322	edges[nbe].adj[1] = 0;
	4323	if(PreviousNewEdge) PreviousNewEdge->adj[1] = & edges[nbe];
	4324	nbe++;
	4325	}
	4326	else nbe += NbEdgeCurve;
	4327	} // end on curve ---
	4328	}
	4329	}
	4330	} // for (i=0;i<nbe;i++)
	4331	if(!step) {
[6412]	4332	_assert_(!edges);
	4333	_assert_(!VerticesOnGeomEdge);
[5448]	4334
[5120]	4335	edges = new Edge[nbex=nbe];
[5448]	4336	if(NbVerticesOnGeomEdge0) VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge0];
	4337
[5120]	4338	// do the vertex on a geometrical vertex
[6412]	4339	_assert_(VerticesOnGeomEdge \|\| NbVerticesOnGeomEdge0==0);
[5120]	4340	NbVerticesOnGeomEdge0 = NbVerticesOnGeomEdge;
	4341	}
[5448]	4342	else{
[6412]	4343	_assert_(NbVerticesOnGeomEdge==NbVerticesOnGeomEdge0);
[5120]	4344	}
	4345	}
	4346
	4347	//Insert points inside existing triangles
[15100]	4348	if (verbose>4) _printf_(" -- current number of vertices = " << nbv << "\n");
[15104]	4349	if (verbose>3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n");
	4350	if (verbose>3) _printf_(" Inserting boundary points\n");
[22085]	4351	Insert();
[5120]	4352
	4353	//Force the boundary
[15104]	4354	if (verbose>3) _printf_(" Forcing boundaries\n");
[5120]	4355	ForceBoundary();
	4356
	4357	//Extract SubDomains
[15104]	4358	if (verbose>3) _printf_(" Extracting subdomains\n");
[5120]	4359	FindSubDomain();
	4360
[15104]	4361	if (verbose>3) _printf_(" Inserting internal points\n");
[5120]	4362	NewPoints(*this,bamgopts,0) ;
[15100]	4363	if (verbose>4) _printf_(" -- current number of vertices = " << nbv << "\n");
[5120]	4364	}
[12365]	4365	/}}}/
[18064]	4366	void Mesh::TriangulateFromGeom1(BamgOpts* bamgopts,int KeepVertices){ /{{{/
[5120]	4367	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/GeomToTriangles1)/
	4368
	4369	/Get options/
	4370	int verbose=bamgopts->verbose;
	4371
[18493]	4372	Gh.NbRef++;// add a ref to Gh
[5120]	4373
[18493]	4374	/*************************************************************************
	4375	* method in 2 steps
	4376	* 1 - compute the number of new edges to allocate
	4377	* 2 - construct the edges
	4378	* remark:
	4379	* in this part we suppose to have a background mesh with the same geometry
	4380	*
	4381	* To construct the discretization of the new mesh we have to
	4382	* rediscretize the boundary of background Mesh
	4383	* because we have only the pointeur from the background mesh to the geometry.
	4384	* We need the abcisse of the background mesh vertices on geometry
	4385	* so a vertex is
	4386	* 0 on GeomVertex ;
	4387	* 1 on GeomEdge + abcisse
	4388	* 2 internal
	4389	*************************************************************************/
[5120]	4390
[18493]	4391	//Check that background mesh and current mesh do have the same geometry
	4392	_assert_(&BTh.Gh==&Gh);
	4393	BTh.NbRef++; // add a ref to BackGround Mesh
[5120]	4394
[18493]	4395	//Initialize new mesh
	4396	BTh.SetVertexFieldOn();
	4397	int* bcurve = new int[Gh.nbcurves]; //
[5120]	4398
[18493]	4399	/* There are 2 ways to make the loop
	4400	* 1) on the geometry
	4401	* 2) on the background mesh
	4402	* if you do the loop on geometry, we don't have the pointeur on background,
	4403	* and if you do the loop in background we have the pointeur on geometry
	4404	* so do the walk on background */
[5120]	4405
[18493]	4406	NbVerticesOnGeomVertex=0;
	4407	NbVerticesOnGeomEdge=0;
[5120]	4408
[18493]	4409	/STEP 1 copy of Required vertices/
[5120]	4410
[18493]	4411	int i;
	4412	for (i=0;i<Gh.nbv;i++) if (Gh[i].Required()) NbVerticesOnGeomVertex++;
	4413	printf("\n");
	4414	if(NbVerticesOnGeomVertex >= maxnbv){
[22085]	4415	delete [] bcurve;
[18493]	4416	_error_("too many vertices on geometry: " << NbVerticesOnGeomVertex << " >= " << maxnbv);
	4417	}
[5120]	4418
[18493]	4419	VerticesOnGeomVertex = new VertexOnGeom[ NbVerticesOnGeomVertex];
	4420	VertexOnBThVertex = new VertexOnVertex[NbVerticesOnGeomVertex];
[5120]	4421
[18493]	4422	//At this point there is NO vertex but vertices should have been allocated by Init
	4423	_assert_(vertices);
	4424	for (i=0;i<Gh.nbv;i++){
	4425	if (Gh[i].Required()) {//Gh vertices Required
	4426	vertices[nbv] =Gh[i];
	4427	vertices[nbv].i=I2(0,0);
	4428	Gh[i].MeshVertexHook = vertices + nbv;// save Geom -> Th
	4429	VerticesOnGeomVertex[nbv]= VertexOnGeom(vertices[nbv],Gh[i]);
[5120]	4430	nbv++;
	4431	}
[18493]	4432	else Gh[i].MeshVertexHook=0;
	4433	}
	4434	for (i=0;i<BTh.NbVerticesOnGeomVertex;i++){
	4435	VertexOnGeom &vog=BTh.VerticesOnGeomVertex[i];
	4436	if (vog.IsRequiredVertex()){
	4437	GeomVertex* gv=vog;
	4438	BamgVertex *bv = vog;
	4439	_assert_(gv->MeshVertexHook); // use of Geom -> Th
	4440	VertexOnBThVertex[NbVertexOnBThVertex++]=VertexOnVertex(gv->MeshVertexHook,bv);
	4441	gv->MeshVertexHook->m = bv->m; // for taking the metric of the background mesh
	4442	}
[5120]	4443	}
[18493]	4444	_assert_(NbVertexOnBThVertex==NbVerticesOnGeomVertex); /This might be due to MaxCornerAngle too small/
[5120]	4445
[18493]	4446	/STEP 2: reseed boundary edges/
[5120]	4447
[18493]	4448	// find the begining of the curve in BTh
	4449	Gh.UnMarkEdges();
	4450	int bfind=0;
	4451	for (int i=0;i<Gh.nbcurves;i++) bcurve[i]=-1;
[5120]	4452
[18493]	4453	/Loop over the backgrounf mesh BTh edges/
	4454	for (int iedge=0;iedge<BTh.nbe;iedge++){
	4455	Edge &ei = BTh.edges[iedge];
[5120]	4456
[18493]	4457	/Loop over the 2 vertices of the current edge/
	4458	for(int je=0;je<2;je++){
[5120]	4459
[18493]	4460	/* If one of the vertex is required we are in a new curve*/
	4461	if (ei[je].GeomEdgeHook->IsRequiredVertex()){
[5120]	4462
[18493]	4463	/Get curve number/
	4464	int nc=ei.GeomEdgeHook->CurveNumber;
[13622]	4465
[18493]	4466	//_printf_("Dealing with curve number " << nc << "\n");
	4467	//_printf_("edge on geometry is same as GhCurve? " << (ei.GeomEdgeHook==Gh.curves[nc].FirstEdge \|\| ei.GeomEdgeHook==Gh.curves[nc].LastEdge)?"yes":"no\n");
	4468	//if(ei.GeomEdgeHook==Gh.curves[nc].FirstEdge \|\| ei.GeomEdgeHook==Gh.curves[nc].LastEdge){
	4469	// _printf_("Do we have the right extremity? curve first vertex -> " << ((GeomVertex )ei[je].GeomEdgeHook==&(*Gh.curves[nc].FirstEdge)[Gh.curves[nc].FirstVertexIndex])?"yes":"no\n");
	4470	// _printf_("Do we have the right extremity? curve last vertex -> " << ((GeomVertex )ei[je].GeomEdgeHook==&(*Gh.curves[nc].LastEdge)[Gh.curves[nc].LastVertexIndex])?"yes":"no\n");
	4471	//}
	4472	//BUG FIX from original bamg
	4473	/Check that we are on the same edge and right vertex (0 or 1) /
	4474	if(ei.GeomEdgeHook==Gh.curves[nc].FirstEdge && (GeomVertex )ei[je].GeomEdgeHook==&(*Gh.curves[nc].FirstEdge)[Gh.curves[nc].FirstVertexIndex]){
	4475	bcurve[nc]=iedge*2+je;
	4476	bfind++;
	4477	}
	4478	else if ((ei.GeomEdgeHook==Gh.curves[nc].LastEdge && (GeomVertex )ei[je].GeomEdgeHook==&(*Gh.curves[nc].LastEdge)[Gh.curves[nc].LastVertexIndex]) && bcurve[nc]==-1){
	4479	bcurve[nc]=iedge*2+je;
	4480	bfind++;
	4481	}
	4482	}
	4483	}
	4484	}
	4485	if (bfind!=Gh.nbcurves){
	4486	delete [] bcurve;
	4487	_error_("problem generating number of curves (" << Gh.nbcurves << " found in the geometry but " << bfind << " curve found in the mesh)");
	4488	}
[5120]	4489
[18493]	4490	// method in 2 + 1 step
	4491	// 0.0) compute the length and the number of vertex to do allocation
	4492	// 1.0) recompute the length
	4493	// 1.1) compute the vertex
[5120]	4494
[18493]	4495	long nbex=0,NbVerticesOnGeomEdgex=0;
	4496	for (int step=0; step <2;step++){
[5120]	4497
[18493]	4498	long NbOfNewPoints=0;
	4499	long NbOfNewEdge=0;
	4500	long iedge;
	4501	Gh.UnMarkEdges();
	4502	double L=0;
[5120]	4503
[18493]	4504	/Go through all geometrical curve/
	4505	for (int icurve=0;icurve<Gh.nbcurves;icurve++){
[5120]	4506
[18493]	4507	/Get edge and vertex (index) of background mesh on this curve/
	4508	iedge=bcurve[icurve]/2;
	4509	int jedge=bcurve[icurve]%2;
[5120]	4510
[18493]	4511	/Get edge of Bth with index iedge/
	4512	Edge &ei = BTh.edges[iedge];
[13622]	4513
[18493]	4514	/Initialize variables/
	4515	double Lstep=0; // step between two points (phase==1)
	4516	long NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1)
[5120]	4517
[18493]	4518	/Do phase 0 to step/
	4519	for(int phase=0;phase<=step;phase++){
[5120]	4520
[18493]	4521	/Current curve pointer/
	4522	Curve *curve= Gh.curves+icurve;
[5120]	4523
[18493]	4524	/Get index of current curve/
	4525	int icurveequi= Gh.GetId(curve);
[5120]	4526
[18493]	4527	/For phase 0, check that we are at the begining of the curve only/
	4528	if(phase==0 && icurveequi!=icurve) continue;
[5120]	4529
[18493]	4530	int k0=jedge,k1;
	4531	Edge* pe= BTh.edges+iedge;
	4532	int iedgeequi=bcurve[icurveequi]/2;
	4533	int jedgeequi=bcurve[icurveequi]%2;
[5120]	4534
[18493]	4535	int k0equi=jedgeequi,k1equi;
	4536	Edge * peequi= BTh.edges+iedgeequi;
	4537	GeomEdge *ongequi = peequi->GeomEdgeHook;
[5120]	4538
[18493]	4539	double sNew=Lstep;// abscisse of the new points (phase==1)
	4540	L=0;// length of the curve
	4541	long i=0;// index of new points on the curve
	4542	GeomVertex * GA0 = (peequi)[k0equi].GeomEdgeHook;
	4543	BamgVertex *A0;
	4544	A0 = GA0->MeshVertexHook; // the vertex in new mesh
	4545	BamgVertex *A1;
	4546	VertexOnGeom *GA1;
	4547	Edge* PreviousNewEdge = 0;
[5120]	4548
[18493]	4549	// New Curve phase
	4550	_assert_(A0-vertices>=0 && A0-vertices<nbv);
	4551	if(ongequi->Required()){
	4552	GeomVertex GA1 = (*peequi)[1-k0equi].GeomEdgeHook;
	4553	A1 = GA1->MeshVertexHook; //
	4554	}
	4555	else {
	4556	for(;;){
	4557	Edge &ee=*pe;
	4558	Edge &eeequi=*peequi;
	4559	k1 = 1-k0; // next vertex of the edge
	4560	k1equi= 1 - k0equi;
	4561	_assert_(pe && ee.GeomEdgeHook);
	4562	ee.GeomEdgeHook->SetMark();
	4563	BamgVertex & v0=ee[0], & v1=ee[1];
	4564	R2 AB=(R2)v1-(R2)v0;
	4565	double L0=L,LAB;
	4566	LAB=LengthInterpole(v0.m,v1.m,AB);
	4567	L+= LAB;
[5120]	4568
[18493]	4569	if (phase){
	4570	// computation of the new points for the given curve
	4571	while ((i!=NbCreatePointOnCurve) && sNew<=L) {
[5120]	4572
[18493]	4573	//some checks
	4574	_assert_(sNew>=L0);
	4575	_assert_(LAB);
	4576	_assert_(vertices && nbv<maxnbv);
	4577	_assert_(edges && nbe<nbex);
	4578	_assert_(VerticesOnGeomEdge && NbVerticesOnGeomEdge<NbVerticesOnGeomEdgex);
[5120]	4579
[18493]	4580	// new vertex on edge
	4581	A1=vertices+nbv++;
	4582	GA1=VerticesOnGeomEdge+NbVerticesOnGeomEdge;
	4583	Edge* e = edges + nbe++;
	4584	double se= (sNew-L0)/LAB;
	4585	if (se<0 \|\| se>=1.000000001){
	4586	_error_("Problem creating point on a boundary: se=" << se << " should be in [0 1]");
[5120]	4587	}
[18493]	4588	se = abscisseInterpole(v0.m,v1.m,AB,se,1);
	4589	if (se<0 \|\| se>1){
	4590	_error_("Problem creating point on a boundary: se=" << se << " should be in [0 1]");
	4591	}
	4592	se = k1 ? se : 1. - se;
	4593	se = k1==k1equi ? se : 1. - se;
	4594	VertexOnBThEdge[NbVerticesOnGeomEdge++] = VertexOnEdge(A1,&eeequi,se); // save
	4595	ongequi=Gh.ProjectOnCurve(eeequi,se,A1,GA1);
	4596	A1->ReferenceNumber = eeequi.ReferenceNumber;
	4597	e->GeomEdgeHook = ongequi;
	4598	e->v[0]=A0;
	4599	e->v[1]=A1;
	4600	e->ReferenceNumber = eeequi.ReferenceNumber;
	4601	e->adj[0]=PreviousNewEdge;
[5120]	4602
[18493]	4603	if (PreviousNewEdge) PreviousNewEdge->adj[1]=e;
	4604	PreviousNewEdge=e;
	4605	A0=A1;
	4606	sNew += Lstep;
	4607	if (++i== NbCreatePointOnCurve) break;
[5120]	4608	}
[5401]	4609	}
[18493]	4610
	4611	//some checks
	4612	_assert_(ee.GeomEdgeHook->CurveNumber==ei.GeomEdgeHook->CurveNumber);
	4613	if (ee[k1].GeomEdgeHook->IsRequiredVertex()) {
	4614	_assert_(eeequi[k1equi].GeomEdgeHook->IsRequiredVertex());
	4615	GeomVertex * GA1 = *eeequi[k1equi].GeomEdgeHook;
	4616	A1=GA1->MeshVertexHook;// the vertex in new mesh
	4617	_assert_(A1-vertices>=0 && A1-vertices<nbv);
	4618	break;
	4619	}
	4620	if (!ee.adj[k1]) {
	4621	_error_("adj edge " << BTh.GetId(ee) << ", nbe=" << nbe << ", Gh.vertices=" << Gh.vertices);
	4622	}
	4623	pe = ee.adj[k1]; // next edge
	4624	k0 = pe->Intersection(ee);
	4625	peequi= eeequi.adj[k1equi]; // next edge
	4626	k0equi=peequi->Intersection(eeequi);
	4627	}// for(;;) end of the curve
[5401]	4628	}
[18493]	4629
	4630	if (phase){ // construction of the last edge
	4631	Edge* e=edges + nbe++;
	4632	e->GeomEdgeHook = ongequi;
	4633	e->v[0]=A0;
	4634	e->v[1]=A1;
	4635	e->ReferenceNumber = peequi->ReferenceNumber;
	4636	e->adj[0]=PreviousNewEdge;
	4637	e->adj[1]=0;
	4638	if (PreviousNewEdge) PreviousNewEdge->adj[1]=e;
	4639	PreviousNewEdge = e;
	4640
	4641	_assert_(i==NbCreatePointOnCurve);
	4642	}
	4643
	4644	if (!phase) { //
	4645	long NbSegOnCurve = Max((long)(L+0.5),(long) 1);// nb of seg
	4646	Lstep = L/NbSegOnCurve;
	4647	NbCreatePointOnCurve = NbSegOnCurve-1;
	4648	NbOfNewEdge += NbSegOnCurve;
	4649	NbOfNewPoints += NbCreatePointOnCurve;
	4650	}
[18492]	4651	}
[18493]	4652	}// end of curve loop
[5120]	4653
[18493]	4654	//Allocate memory
	4655	if(step==0){
	4656	if(nbv+NbOfNewPoints > maxnbv) {
	4657	_error_("too many vertices on geometry: " << nbv+NbOfNewPoints << " >= " << maxnbv);
	4658	}
	4659	edges = new Edge[NbOfNewEdge];
	4660	nbex = NbOfNewEdge;
	4661	if(NbOfNewPoints) {
	4662	VerticesOnGeomEdge = new VertexOnGeom[NbOfNewPoints];
	4663	NbVertexOnBThEdge = NbOfNewPoints;
	4664	VertexOnBThEdge = new VertexOnEdge[NbOfNewPoints];
	4665	NbVerticesOnGeomEdgex = NbOfNewPoints;
	4666	}
	4667	NbOfNewPoints =0;
	4668	NbOfNewEdge = 0;
[5120]	4669	}
	4670	}
[18493]	4671	_assert_(nbe!=0);
	4672	delete [] bcurve;
[5120]	4673
	4674	//Insert points inside existing triangles
[15100]	4675	if (verbose>4) _printf_(" -- current number of vertices = " << nbv << "\n");
[15104]	4676	if (verbose>3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n");
	4677	if (verbose>3) _printf_(" Inserting boundary points\n");
[22085]	4678	Insert();
[5120]	4679
	4680	//Force the boundary
[15104]	4681	if (verbose>3) _printf_(" Forcing boundaries\n");
[5120]	4682	ForceBoundary();
	4683
	4684	//Extract SubDomains
[15104]	4685	if (verbose>3) _printf_(" Extracting subdomains\n");
[5120]	4686	FindSubDomain();
	4687
[15104]	4688	if (verbose>3) _printf_(" Inserting internal points\n");
[5120]	4689	NewPoints(BTh,bamgopts,KeepVertices) ;
[15100]	4690	if (verbose>4) _printf_(" -- current number of vertices = " << nbv << "\n");
[5120]	4691	}
[12365]	4692	/}}}/
[22085]	4693	int Mesh::RandomNumber(int max){/{{{/
	4694	/* Generate a random number from 0 to max-1 using linear congruential
	4695	* random number generator*/
[3913]	4696
[22085]	4697	this->randomseed = (this->randomseed * 1366l + 150889l) % 714025l;
	4698	return this->randomseed / (714025l / max + 1);
	4699	} /}}}/
	4700	int Mesh::ForceEdge(BamgVertex &a, BamgVertex & b,AdjacentTriangle & taret) { /{{{/
[18064]	4701	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ForceEdge)/
[3913]	4702
[18064]	4703	int NbSwap =0;
	4704	if (!a.t \|\| !b.t){ // the 2 vertex is in a mesh
	4705	_error_("!a.t \|\| !b.t");
	4706	}
	4707	int k=0;
	4708	taret=AdjacentTriangle(0,0); // erreur
[3913]	4709
[18064]	4710	AdjacentTriangle tta(a.t,EdgesVertexTriangle[a.IndexInTriangle][0]);
	4711	BamgVertex v1, v2 = tta.EdgeVertex(0),*vbegin =v2;
	4712	// we turn around a in the direct direction
[3913]	4713
[22085]	4714	long long det2 = v2 ? det(*v2,a,b): -1 , det1;
[18064]	4715	if(v2) // normal case
	4716	det2 = det(*v2,a,b);
	4717	else { // no chance infini vertex try the next
	4718	tta= Previous(Adj(tta));
	4719	v2 = tta.EdgeVertex(0);
	4720	vbegin =v2;
	4721	if (!v2){
	4722	_error_("!v2");
	4723	}
	4724	det2 = det(*v2,a,b);
[3913]	4725	}
	4726
[18064]	4727	while (v2 != &b) {
	4728	AdjacentTriangle tc = Previous(Adj(tta));
	4729	v1 = v2;
	4730	v2 = tc.EdgeVertex(0);
	4731	det1 = det2;
	4732	det2 = v2 ? det(*v2,a,b): det2;
[3913]	4733
[18064]	4734	if((det1 < 0) && (det2 >0)) {
	4735	// try to force the edge
	4736	BamgVertex * va = &a, *vb = &b;
	4737	tc = Previous(tc);
	4738	if (!v1 \|\| !v2){
	4739	_error_("!v1 \|\| !v2");
	4740	}
[22085]	4741	long long detss = 0,l=0;
	4742	while ((this->SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap)))
[18064]	4743	if(l++ > 10000000) {
	4744	_error_("Loop in forcing Egde, nb de swap=" << NbSwap << ", nb of try swap (" << l << ") too big");
	4745	}
	4746	BamgVertex aa = tc.EdgeVertex(0), bb = tc.EdgeVertex(1);
	4747	if (((aa == &a ) && (bb == &b)) \|\|((bb == &a ) && (aa == &b))){
	4748	tc.SetLock();
	4749	a.Optim(1,0);
	4750	b.Optim(1,0);
	4751	taret = tc;
	4752	return NbSwap;
	4753	}
	4754	else
	4755	{
	4756	taret = tc;
	4757	return -2; // error boundary is crossing
	4758	}
[3913]	4759	}
[18064]	4760	tta = tc;
	4761	k++;
	4762	if (k>=2000){
	4763	_error_("k>=2000");
[3913]	4764	}
[18064]	4765	if ( vbegin == v2 ) return -1;// error
[3913]	4766	}
[18064]	4767
	4768	tta.SetLock();
	4769	taret=tta;
	4770	a.Optim(1,0);
	4771	b.Optim(1,0);
	4772	return NbSwap;
[3913]	4773	}
[18064]	4774	/}}}/
[22085]	4775	int Mesh::SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb ,AdjacentTriangle & tt1,long long & dets1, long long & detsa,long long & detsb, int & NbSwap) {/{{{/
[3913]	4776	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SwapForForcingEdge)/
	4777	// l'arete ta coupe l'arete pva pvb
	4778	// de cas apres le swap sa coupe toujours
	4779	// on cherche l'arete suivante
	4780	// on suppose que detsa >0 et detsb <0
	4781	// attention la routine echange pva et pvb
	4782
	4783	if(tt1.Locked()) return 0; // frontiere croise
	4784
[5143]	4785	AdjacentTriangle tt2 = Adj(tt1);
[3913]	4786	Triangle t1=tt1,t2=tt2;// les 2 triangles adjacent
	4787	short a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles
	4788	if ( a1<0 \|\| a1>=3 ){
[13036]	4789	_error_("a1<0 \|\| a1>=3");
[3913]	4790	}
	4791
[5120]	4792	BamgVertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]];
	4793	BamgVertex & s1= (*t1)[OppositeVertex[a1]];
	4794	BamgVertex & s2= (*t2)[OppositeVertex[a2]];
[3913]	4795
[22085]	4796	long long dets2 = det(pva,pvb,s2);
	4797	long long det1=t1->det , det2=t2->det ;
	4798	long long detT = det1+det2;
[3913]	4799	if ((det1<=0 ) \|\| (det2<=0)){
[13036]	4800	_error_("(det1<=0 ) \|\| (det2<=0)");
[3913]	4801	}
	4802	if ( (detsa>=0) \|\| (detsb<=0) ){ // [a,b] cut infinite line va,bb
[13036]	4803	_error_("(detsa>=0) \|\| (detsb<=0)");
[3913]	4804	}
[22085]	4805	long long ndet1 = bamg::det(s1,sa,s2);
	4806	long long ndet2 = detT - ndet1;
[3913]	4807
	4808	int ToSwap =0; //pas de swap
	4809	if ((ndet1 >0) && (ndet2 >0))
	4810	{ // on peut swaper
	4811	if ((dets1 <=0 && dets2 <=0) \|\| (dets2 >=0 && detsb >=0))
	4812	ToSwap =1;
	4813	else // swap alleatoire
[22085]	4814	if (this->RandomNumber(1))
[3913]	4815	ToSwap =2;
	4816	}
	4817	if (ToSwap) NbSwap++,
	4818	bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2);
	4819
	4820	int ret=1;
	4821
	4822	if (dets2 < 0) {// haut
	4823	dets1 = ToSwap ? dets1 : detsa ;
	4824	detsa = dets2;
	4825	tt1 = Previous(tt2) ;}
	4826	else if (dets2 > 0){// bas
	4827	dets1 = ToSwap ? dets1 : detsb ;
	4828	detsb = dets2;
	4829	//xxxx tt1 = ToSwap ? tt1 : Next(tt2);
	4830	if(!ToSwap) tt1 = Next(tt2);
	4831	}
[5148]	4832	else { // changement de direction
[3913]	4833	ret = -1;
	4834	Exchange(pva,pvb);
	4835	Exchange(detsa,detsb);
	4836	Exchange(dets1,dets2);
	4837	Exchange(tt1,tt2);
	4838	dets1=-dets1;
	4839	dets2=-dets2;
	4840	detsa=-detsa;
	4841	detsb=-detsb;
	4842
[16521]	4843	if(ToSwap){
[18064]	4844	if (dets2 < 0) {// haut
	4845	dets1 = (ToSwap ? dets1 : detsa) ;
	4846	detsa = dets2;
	4847	tt1 = Previous(tt2) ;}
	4848	else if(dets2 > 0){// bas
	4849	dets1 = (ToSwap ? dets1 : detsb) ;
	4850	detsb = dets2;
	4851	if(!ToSwap) tt1 = Next(tt2);
	4852	}
	4853	else {// on a fin ???
	4854	tt1 = Next(tt2);
	4855	ret =0;}
[16521]	4856	}
[3913]	4857
	4858	}
	4859	return ret;
	4860	}
[12365]	4861	/}}}/
[22085]	4862
	4863	/Intermediary/
	4864	AdjacentTriangle CloseBoundaryEdge(I2 A,Triangle t, double &a,double &b) {/{{{*/
	4865	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/CloseBoundaryEdge)/
	4866
	4867	int k=(t)(0) ? (( (t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0;
	4868	int dir=0;
	4869	if (k<0){
	4870	_error_("k<0");
	4871	}
	4872	int kkk=0;
	4873	long long IJ_IA,IJ_AJ;
	4874	AdjacentTriangle edge(t,OppositeEdge[k]);
	4875	for (;;edge = dir >0 ? Next(Adj(Next(edge))) : Previous(Adj(Previous(edge)))) {
	4876	kkk++;
	4877	if (kkk>=1000){
	4878	_error_("kkk>=1000");
	4879	}
	4880	BamgVertex &vI = *edge.EdgeVertex(0);
	4881	BamgVertex &vJ = *edge.EdgeVertex(1);
	4882	I2 I=vI, J=vJ, IJ= J-I;
	4883	IJ_IA = (IJ ,(A-I));
	4884	if (IJ_IA<0) {
	4885	if (dir>0) {a=1;b=0;return edge;}// change of signe => I
	4886	else {dir=-1;
	4887	continue;}};// go in direction i
	4888	IJ_AJ = (IJ ,(J-A));
	4889	if (IJ_AJ<0) {
	4890	if(dir<0) {a=0;b=1;return edge;}
	4891	else {dir = 1;
	4892	continue;}}// go in direction j
	4893	double IJ2 = IJ_IA + IJ_AJ;
	4894	if (IJ2==0){
	4895	_error_("IJ2==0");
	4896	}
	4897	a= IJ_AJ/IJ2;
	4898	b= IJ_IA/IJ2;
	4899	return edge;
	4900	}
	4901	}
	4902	/}}}/
	4903	void swap(Triangle t1,short a1, Triangle t2,short a2, BamgVertex s1,BamgVertex s2,long long det1,long long det2){ /{{{/
	4904	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/swap)/
	4905	// --------------------------------------------------------------
	4906	// short a2=aa[a];// les 2 numero de l arete dans les 2 triangles
	4907	//
	4908	// sb sb
	4909	// / \| \ / \ !
	4910	// as1/ \| \ /a2 \ !
	4911	// / \| \ / t2 \ !
	4912	// s1 /t1 \| t2 \s2 --> s1 /___as2___\s2 !
	4913	// \ a1\|a2 / \ as1 /
	4914	// \ \| / \ t1 /
	4915	// \ \| / as2 \ a1/
	4916	// \ \| / \ /
	4917	// sa sa
	4918	// -------------------------------------------------------------
	4919	int as1 = NextEdge[a1];
	4920	int as2 = NextEdge[a2];
	4921	int ap1 = PreviousEdge[a1];
	4922	int ap2 = PreviousEdge[a2];
	4923	(*t1)(VerticesOfTriangularEdge[a1][1]) = s2 ; // avant sb
	4924	(*t2)(VerticesOfTriangularEdge[a2][1]) = s1 ; // avant sa
	4925	// mise a jour des 2 adjacences externes
	4926	AdjacentTriangle taas1 = t1->Adj(as1),
	4927	taas2 = t2->Adj(as2),
	4928	tas1(t1,as1), tas2(t2,as2),
	4929	ta1(t1,a1),ta2(t2,a2);
	4930	// externe haut gauche
	4931	taas1.SetAdj2(ta2, taas1.GetAllFlag_UnSwap());
	4932	// externe bas droite
	4933	taas2.SetAdj2(ta1, taas2.GetAllFlag_UnSwap());
	4934	// remove the Mark UnMarkSwap
	4935	t1->SetUnMarkUnSwap(ap1);
	4936	t2->SetUnMarkUnSwap(ap2);
	4937	// interne
	4938	tas1.SetAdj2(tas2);
	4939
	4940	t1->det = det1;
	4941	t2->det = det2;
	4942
	4943	t1->SetSingleVertexToTriangleConnectivity();
	4944	t2->SetSingleVertexToTriangleConnectivity();
	4945	} // end swap
	4946	/}}}/
[3913]	4947	}

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

Context Navigation

source: issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Mesh.cpp@ 22085

Download in other formats: