Context Navigation

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

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

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

Download in other formats: