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

Line
1	#include <cstdio>
2	#include <cstring>
3	#include <cmath>
4	#include <ctime>
5
6	#include "../shared/shared.h"
7	#include "./bamgobjects.h"
8	#include "./det.h"
9
10	namespace bamg {
11
12	/Constructors/Destructors/
13	Mesh::Mesh(BamgGeom* bamggeom,BamgMesh* bamgmesh, BamgOpts* bamgopts):Gh((new Geometry())),BTh(this){ /{{{/
14
15	/Initialize fields/
16	Init(0);
17
18	/Read Geometry if provided/
19	if(bamggeom->Edges) {
20	Gh.ReadGeometry(bamggeom,bamgopts);
21	Gh.PostRead();
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/
30	_printf_("WARNING: mesh present but no geometry found. Reconstructing...\n");
31	BuildGeometryFromMesh(bamgopts);
32	Gh.PostRead();
33	}
34
35	/Set integer coordinates/
36	SetIntCoor();
37
38	/Fill holes and generate mesh properties/
39	ReconstructExistingMesh();
40	}
41	/}}}/
42	Mesh::Mesh(int* index,double* x,double* y,int nods,int nels):Gh((new Geometry())),BTh(this){/{{{/
43
44	Init(0);
45	ReadMesh(index,x,y,nods,nels);
46	SetIntCoor();
47	ReconstructExistingMesh();
48	}
49	/}}}/
50	Mesh::Mesh(double* x,double* y,int nods):Gh((new Geometry())),BTh(this){/{{{/
51	Triangulate(x,y,nods);
52	}
53	/}}}/
54	Mesh::Mesh(const Mesh & Tho,const int flag ,const int bb,BamgOpts* bamgopts) : Gh((new Geometry())), BTh(this) {/{{{/
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++;
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));
79	if ( reft[itadj] >=0 && !flag[itadj])
80	{ nbNewBedge++;
81	refv[Tho.GetId(t[VerticesOfTriangularEdge[0][0]])]=bb[i];
82	refv[Tho.GetId(t[VerticesOfTriangularEdge[0][1]])]=bb[i];
83	}
84	itadj=Tho.GetId(t.TriangleAdj(1));
85	if ( reft[itadj] >=0 && !flag[itadj])
86	{ nbNewBedge++;
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));
90	if ( reft[itadj] >=0 && !flag[itadj])
91	{ nbNewBedge++;
92	refv[Tho.GetId(t[VerticesOfTriangularEdge[2][0]])]=bb[i];
93	refv[Tho.GetId(t[VerticesOfTriangularEdge[2][1]])]=bb[i];}
94	}
95	k=0;
96	for (i=0;i<Tho.nbv;i++){
97	if (kk[i]>=0) kk[i]=k++;
98	}
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");
102	long imaxnbv =k;
103	Init(imaxnbv);
104	for (i=0;i<Tho.nbv;i++)
105	if (kk[i]>=0)
106	{
107	vertices[nbv] = Tho.vertices[i];
108	if (!vertices[nbv].GetReferenceNumber())
109	vertices[nbv].ReferenceNumber = refv[i];
110	nbv++;
111	}
112	if (imaxnbv != nbv){
113	delete [] kk;
114	delete [] refv;
115	_error_("imaxnbv != nbv");
116	}
117	for (i=0;i<Tho.nbt;i++)
118	if(reft[i] >=0 && flag[i]){
119	const Triangle & t = Tho.triangles[i];
120	int i0 = Tho.GetId(t[0]);
121	int i1 = Tho.GetId(t[1]);
122	int i2 = Tho.GetId(t[2]);
123	if(i0<0 \|\| i1<0 \|\| i2<0){
124	delete [] refv;
125	_error_("i0<0 \|\| i1<0 \|\| i2< 0");
126	}
127	if(i0>=Tho.nbv \|\| i1>=Tho.nbv \|\| i2>=Tho.nbv){
128	delete [] refv;
129	_error_("i0>=Tho.nbv \|\| i1>=Tho.nbv \|\| i2>=Tho.nbv");
130	}
131	triangles[nbt] = Triangle(this,kk[i0],kk[i1],kk[i2]);
132	triangles[nbt].color = Tho.subdomains[reft[i]].ReferenceNumber;
133	nbt++;
134	}
135	if (nbt==0 && nbv==0) {
136	_error_("All triangles have been removed");
137	}
138	delete [] kk;
139	delete [] reft;
140	delete [] refv;
141	BuildGeometryFromMesh(bamgopts);
142	Gh.PostRead();
143	SetIntCoor();
144	ReconstructExistingMesh();
145
146	/Final checks/
147	_assert_(kt==nbt);
148	_assert_(nbsubdomains);
149	_assert_(subdomains[0].head && subdomains[0].head->link);
150	}
151	/}}}/
152	Mesh::Mesh(Mesh & Th,Geometry * pGh,Mesh * pBth,long maxnbv_in)/{{{/
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++;
156	maxnbv_in = Max(maxnbv_in,Th.nbv);
157	long i;
158	// do all the allocation to be sure all the pointer existe
159
160	Init(maxnbv_in);// to make the allocation
161	// copy of triangles
162	nbv = Th.nbv;
163	nbt = Th.nbt;
164	nbe = Th.nbe;
165	nbsubdomains = Th.nbsubdomains;
166	nbtout = Th.nbtout;
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];
192	if(nbsubdomains)
193	subdomains = new SubDomain[nbsubdomains];
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);
203	for(i=0;i<nbsubdomains;i++)
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	}
212	/}}}/
213	Mesh::Mesh(long imaxnbv,Mesh & BT,BamgOpts* bamgopts,int keepBackVertices) :Gh(BT.Gh),BTh(BT) {/{{{/
214	this->Init(imaxnbv);
215	TriangulateFromGeom1(bamgopts,keepBackVertices);
216	}
217	/}}}/
218	Mesh::Mesh(long imaxnbv,Geometry & G,BamgOpts* bamgopts):Gh(G),BTh(this){/{{{*/
219	Init(imaxnbv);
220	TriangulateFromGeom0(bamgopts);
221	}
222	/}}}/
223	Mesh::~Mesh() {/{{{/
224	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Triangles)/
225
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;
236
237	if (Gh.NbRef>0) Gh.NbRef--;
238	else if (Gh.NbRef==0) delete &Gh;
239	if(&BTh != this){
240	if (BTh.NbRef>0) BTh.NbRef--;
241	else if (BTh.NbRef==0) delete &BTh;
242	}
243	Init(0); // set all to zero
244	}
245	/}}}/
246
247	/IO/
248	void Mesh::ReadMesh(int* index,double* x,double* y,int nods,int nels){/{{{/
249
250	long i1,i2,i3;
251	long i;
252	Metric M1(1);
253	int verbose=0;
254	bool* nodeflags=NULL;
255
256	nbv=nods;
257	maxnbv=nbv;
258	nbt=nels;
259
260	//Vertices
261	if (verbose) _printf_("Reading vertices (" << nbv << ")\n");
262	vertices=xNew<BamgVertex>(nbv);
263	orderedvertices=xNew<BamgVertex*>(nbv);
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	}
271	maxnbt=2*maxnbv-2; // for filling The Holes and quadrilaterals
272
273	//Triangles
274	if (verbose) _printf_("Reading triangles (" << nbt << ")\n");
275	triangles =new Triangle[maxnbt]; //we cannot allocate only nbt triangles since
276	nodeflags=xNew<bool>(nbv);
277	for(i=0;i<nbv;i++) nodeflags[i]=false;
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;
286	nodeflags[i1]=nodeflags[i2]=nodeflags[i3]=true;
287	}
288
289	/Recreate geometry: /
290	if (verbose) _printf_("Building Geometry\n");
291	BuildGeometryFromMesh();
292	if (verbose) _printf_("Completing geometry\n");
293	Gh.PostRead();
294
295	/Check that there is no orphan/
296	bool isorphan=false;
297	for(i=0;i<nbv;i++){
298	if(!nodeflags[i]){
299	_printf_("Vertex " << i+1 << " does not belong to any element\n");
300	isorphan=true;
301	}
302	}
303	if(isorphan) _error_("Orphan found in mesh, see ids above");
304
305	/Clean up/
306	xDelete<bool>(nodeflags);
307	}
308	/}}}/
309	void Mesh::ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts){/{{{/
310
311	int verbose;
312	double Hmin = HUGE_VAL; // the infinie value
313	long i1,i2,i3;
314	long i,j;
315	Metric M1(1);
316
317	verbose=bamgopts->verbose;
318
319	nbv=bamgmesh->VerticesSize[0];
320	maxnbv=nbv;
321	nbt=bamgmesh->TrianglesSize[0];
322
323	//Vertices
324	if(bamgmesh->Vertices){
325	if(verbose>5) _printf_(" processing Vertices\n");
326
327	vertices=xNew<BamgVertex>(nbv);
328	orderedvertices=xNew<BamgVertex*>(nbv);
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	}
337	maxnbt=2*maxnbv-2; // for filling The Holes and quadrilaterals
338	}
339	else{
340	if(verbose>5) _error_("no Vertices found in the initial mesh");
341	}
342
343	//Triangles
344	if(bamgmesh->Triangles){
345	if(verbose>5) _printf_(" processing Triangles\n");
346	triangles =new Triangle[maxnbt]; //we cannot allocate only nbt triangles since
347	//other triangles will be added for each edge
348	for (i=0;i<nbt;i++){
349	Triangle &t=triangles[i];
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{
358	if(verbose>5) _error_("no Triangles found in the initial mesh");
359	}
360
361	//VerticesOnGeomEdge
362	if(bamgmesh->VerticesOnGeomEdge){
363	if(verbose>5) _printf_(" processing VerticesOnGeomEdge\n");
364	NbVerticesOnGeomEdge=bamgmesh->VerticesOnGeomEdgeSize[0];
365	VerticesOnGeomEdge= new VertexOnGeom[NbVerticesOnGeomEdge] ;
366	for (i=0;i<NbVerticesOnGeomEdge;i++){
367	long i1,i2;
368	double s;
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];
372	VerticesOnGeomEdge[i]=VertexOnGeom(vertices[i1],Gh.edges[i2],s);
373	}
374	}
375
376	//VerticesOnGeomVertex
377	if(bamgmesh->VerticesOnGeomVertexSize[0]){
378	if(verbose>5) _printf_(" processing VerticesOnGeomVertex\n");
379	NbVerticesOnGeomVertex=bamgmesh->VerticesOnGeomVertexSize[0];
380	VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex] ;
381	for (i=0;i<NbVerticesOnGeomVertex;i++){
382	long i1,i2;
383	i1=(long)bamgmesh->VerticesOnGeomVertex[i*2+0]-1; //for C indexing
384	i2=(long)bamgmesh->VerticesOnGeomVertex[i*2+1]-1; //for C indexing
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
394	if(verbose>5) _printf_(" processing Edges\n");
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
404	edges[i].ReferenceNumber=(long)bamgmesh->Edges[i*3+2];
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++) {
435	BamgVertex *v=edges[i].v[j];
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;
443	_assert_(v==edges[i0 ].v[j0]);
444	edges[i ].adj[j ] =edges +i0;
445	edges[i0].adj[j0] =edges +i ;
446	v->color = -3;
447	}
448	}
449	}
450	}
451
452	//EdgeOnGeomEdge
453	if(bamgmesh->EdgesOnGeomEdge){
454	if(verbose>5) _printf_(" processing EdgesOnGeomEdge\n");
455	int i1,i2,i,j;
456	i2=bamgmesh->EdgesOnGeomEdgeSize[0];
457	for (i1=0;i1<i2;i1++) {
458	i=(int)bamgmesh->EdgesOnGeomEdge[i1*2+0]-1; //C indexing
459	j=(int)bamgmesh->EdgesOnGeomEdge[i1*2+1]-1; //C indexing
460	//Check value
461	if(!(i>=0 && j>=0 && i<nbe && j<Gh.nbe)) {
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 << "]");
463	}
464	edges[i].GeomEdgeHook=Gh.edges+j;
465	}
466	}
467
468	//SubDomain
469	if(bamgmesh->SubDomains){
470	long i3,head,direction;
471	if(verbose>5) _printf_(" processing SubDomains\n");
472	nbsubdomains=bamgmesh->SubDomainsSize[0];
473	subdomains = new SubDomain [ nbsubdomains ];
474	for (i=0;i<nbsubdomains;i++) {
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];
478	if (i3!=3) _error_("Bad Subdomain definition: first number should be 3");
479	if (head<0 \|\| head>=nbt) _error_("Bad Subdomain definition: head should in [1 " << nbt << "] (triangle number)");
480	subdomains[i].head = triangles+head;
481	subdomains[i].direction = direction;
482	subdomains[i].ReferenceNumber = i3;
483	}
484	}
485
486	}
487	/}}}/
488	void Mesh::WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts){/{{{/
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
510	head_1=xNew<int>(nbv);
511	next_1=xNew<int>(3*nbt);
512	connectivitysize_1=xNew<int>(nbv);
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++){
523	int v=GetId(triangles[i][j]); //jth vertex of the ith triangle
524	if (k>3*nbt-1 \|\| k<0) _error_("k = " << k << ", nbt = " << nbt);
525	next_1[k]=head_1[v];
526	if (v>nbv-1 \|\| v<0) _error_("v = " << v << ", nbv = " << nbv);
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/
541	if(verbose>5) _printf_(" writing Vertices\n");
542	bamgmesh->VerticesSize[0]=nbv;
543	bamgmesh->VerticesSize[1]=3;
544	if(nbv){
545	bamgmesh->Vertices=xNew<double>(3*nbv);
546	bamgmesh->PreviousNumbering=xNew<double>(nbv);
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;
550	bamgmesh->Vertices[i*3+2]=vertices[i].GetReferenceNumber();
551	bamgmesh->PreviousNumbering[i]=vertices[i].PreviousNumber;
552	}
553	}
554
555	/Edges/
556	if(verbose>5) _printf_(" writing Edges\n");
557	bamgmesh->EdgesSize[0]=nbe;
558	bamgmesh->EdgesSize[1]=3;
559	int NumIssmSegments=0;
560	if (nbe){
561	bamgmesh->Edges=xNew<double>(3*nbe);
562	for (i=0;i<nbe;i++){
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
565	bamgmesh->Edges[i*3+2]=edges[i].ReferenceNumber;
566	if(edges[i].GeomEdgeHook){
567	NumIssmSegments++;
568	}
569	}
570	}
571
572	/Element edges/
573	if(verbose>5) _printf_(" writing element edges\n");
574	SetOfEdges4* edge4=new SetOfEdges4(nbt*3,nbv);
575	double* elemedge=NULL;
576	elemedge=xNew<double>(3*nbt);
577	for (i=0;i<3*nbt;i++) elemedge[i]=-2.;//will become -1
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++) {
583	i1=GetId(triangles[i][VerticesOfTriangularEdge[j][0]]);
584	i2=GetId(triangles[i][VerticesOfTriangularEdge[j][1]]);
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	}
599	bamgmesh->IssmEdgesSize[0]=edge4->nb();
600	bamgmesh->IssmEdgesSize[1]=4;
601	bamgmesh->IssmEdges=xNew<double>(4*edge4->nb());
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
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
611	}
612	else{
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
615	}
616	found=true;
617	break;
618	}
619	}
620	_assert_(found);
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
623	}
624	//clean up
625	delete edge4;
626	xDelete<double>(elemedge);
627
628	/IssmSegments/
629	if(verbose>5) _printf_(" writing IssmSegments\n");
630	bamgmesh->IssmSegmentsSize[0]=NumIssmSegments;
631	bamgmesh->IssmSegmentsSize[1]=4;
632	bamgmesh->IssmSegments=xNew<double>(4*NumIssmSegments);
633	num=0;
634	for (i=0;i<nbe;i++){
635	if(edges[i].GeomEdgeHook){
636	//build segment
637	int i1=GetId(edges[i][0]);
638	int i2=GetId(edges[i][1]);
639	bool stop=false;
640	for(j=head_1[i1];j!=-1;j=next_1[j]){
641	for(k=0;k<3;k++){
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
646	bamgmesh->IssmSegments[num*4+2]=(int)j/3+1; //back to M indexing
647	bamgmesh->IssmSegments[num*4+3]=edges[i].ReferenceNumber;
648	num+=1;
649	stop=true;
650	break;
651	}
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
655	bamgmesh->IssmSegments[num*4+2]=(int)j/3+1; //back to M indexing
656	bamgmesh->IssmSegments[num*4+3]=edges[i].ReferenceNumber;
657	num+=1;
658	stop=true;
659	break;
660	}
661	}
662	}
663	if(stop) break;
664	}
665	if (!stop){
666	_error_("Element holding segment [" << i1+1 << " " << i2+1 << "] not found...");
667	}
668	}
669	}
670
671	/Triangles/
672	if(verbose>5) _printf_(" writing Triangles\n");
673	k=nbInT;
674	num=0;
675	bamgmesh->TrianglesSize[0]=k;
676	bamgmesh->TrianglesSize[1]=4;
677	if (k){
678	bamgmesh->Triangles=xNew<double>(4*k);
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) )){
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
686	bamgmesh->Triangles[num*4+3]=subdomains[reft[i]].ReferenceNumber;
687	num=num+1;
688	}
689	}
690	}
691
692	/SubDomains/
693	if(verbose>5) _printf_(" writing SubDomains\n");
694	bamgmesh->SubDomainsSize[0]=nbsubdomains;
695	bamgmesh->SubDomainsSize[1]=4;
696	if (nbsubdomains){
697	bamgmesh->SubDomains=xNew<double>(4*nbsubdomains);
698	for (i=0;i<nbsubdomains;i++){
699	bamgmesh->SubDomains[i*4+0]=3;
700	bamgmesh->SubDomains[i*4+1]=reft[GetId(subdomains[i].head)];
701	bamgmesh->SubDomains[i*4+2]=1;
702	bamgmesh->SubDomains[i*4+3]=subdomains[i].ReferenceNumber;
703	}
704	}
705
706	/SubDomainsFromGeom/
707	if(verbose>5) _printf_(" writing SubDomainsFromGeom\n");
708	bamgmesh->SubDomainsFromGeomSize[0]=Gh.nbsubdomains;
709	bamgmesh->SubDomainsFromGeomSize[1]=4;
710	if (Gh.nbsubdomains){
711	bamgmesh->SubDomainsFromGeom=xNew<double>(4*Gh.nbsubdomains);
712	for (i=0;i<Gh.nbsubdomains;i++){
713	bamgmesh->SubDomainsFromGeom[i*4+0]=2;
714	bamgmesh->SubDomainsFromGeom[i*4+1]=GetId(subdomains[i].edge)+1; //back to Matlab indexing
715	bamgmesh->SubDomainsFromGeom[i*4+2]=subdomains[i].direction;
716	bamgmesh->SubDomainsFromGeom[i*4+3]=Gh.subdomains[i].ReferenceNumber;
717	}
718	}
719
720	/VerticesOnGeomVertex/
721	if(verbose>5) _printf_(" writing VerticesOnGeomVertex\n");
722	bamgmesh->VerticesOnGeomVertexSize[0]=NbVerticesOnGeomVertex;
723	bamgmesh->VerticesOnGeomVertexSize[1]=2;
724	if (NbVerticesOnGeomVertex){
725	bamgmesh->VerticesOnGeomVertex=xNew<double>(2*NbVerticesOnGeomVertex);
726	for (i=0;i<NbVerticesOnGeomVertex;i++){
727	VertexOnGeom &v=VerticesOnGeomVertex[i];
728	_assert_(v.OnGeomVertex());
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
731	}
732	}
733
734	/VertexOnGeomEdge/
735	if(verbose>5) _printf_(" writing VerticesOnGeomEdge\n");
736	bamgmesh->VerticesOnGeomEdgeSize[0]=NbVerticesOnGeomEdge;
737	bamgmesh->VerticesOnGeomEdgeSize[1]=3;
738	if (NbVerticesOnGeomEdge){
739	bamgmesh->VerticesOnGeomEdge=xNew<double>(3*NbVerticesOnGeomEdge);
740	for (i=0;i<NbVerticesOnGeomEdge;i++){
741	const VertexOnGeom &v=VerticesOnGeomEdge[i];
742	if (!v.OnGeomEdge()){
743	_error_("A vertices supposed to be OnGeomEdge is actually not");
744	}
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
748	}
749	}
750
751	/EdgesOnGeomEdge/
752	if(verbose>5) _printf_(" writing EdgesOnGeomEdge\n");
753	k=0;
754	for (i=0;i<nbe;i++){
755	if (edges[i].GeomEdgeHook) k=k+1;
756	}
757	bamgmesh->EdgesOnGeomEdgeSize[0]=k;
758	bamgmesh->EdgesOnGeomEdgeSize[1]=2;
759	if (k){
760	bamgmesh->EdgesOnGeomEdge=xNew<double>(2*(int)k);
761	int count=0;
762	for (i=0;i<nbe;i++){
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
766	count=count+1;
767	}
768	}
769	}
770
771	/Element Connectivity/
772	if(verbose>5) _printf_(" writing Element connectivity\n");
773	bamgmesh->ElementConnectivitySize[0]=nbt-nbtout;
774	bamgmesh->ElementConnectivitySize[1]=3;
775	bamgmesh->ElementConnectivity=xNew<double>(3*(nbt-nbtout));
776	for (i=0;i<3*(nbt-nbtout);i++) bamgmesh->ElementConnectivity[i]=NAN;
777	num=0;
778	for (i=0;i<nbt;i++){
779	if (reft[i]>=0){
780	for (j=0;j<3;j++){
781	k=GetId(triangles[i].TriangleAdj(j));
782	if (reft[k]>=0){
783	_assert_(3num+j<3(nbt-nbtout));
784	bamgmesh->ElementConnectivity[3*num+j]=k+1; // back to Matlab indexing
785	}
786	}
787	num+=1;
788	}
789	}
790
791	/ElementNodal Connectivity/
792	if(verbose>5) _printf_(" writing Nodal element connectivity\n");
793	bamgmesh->NodalElementConnectivitySize[0]=nbv;
794	bamgmesh->NodalElementConnectivitySize[1]=connectivitymax_1;
795	bamgmesh->NodalElementConnectivity=xNew<double>(connectivitymax_1*nbv);
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]){
800	_assert_(connectivitymax_1i+k < connectivitymax_1nbv);
801	bamgmesh->NodalElementConnectivity[connectivitymax_1*i+k]=floor((double)j/3)+1;
802	k++;
803	}
804	}
805
806	/Nodal Connectivity/
807	if(verbose>5) _printf_(" writing Nodal connectivity\n");
808	//chaining algorithm (again...)
809	int* head_2=NULL;
810	int* next_2=NULL;
811	int* connectivitysize_2=NULL;
812	int connectivitymax_2=0;
813	i1=bamgmesh->IssmEdgesSize[0];
814	i2=bamgmesh->IssmEdgesSize[1];
815	head_2=xNew<int>(nbv);
816	next_2=xNew<int>(2*i1);
817	connectivitysize_2=xNew<int>(nbv);
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++){
825	int v=(int)bamgmesh->IssmEdges[i*i2+j]-1; //back to C indexing
826	if (k>2i1-1 \|\| k<0) _error_("Index exceed matrix dimensions (k=" << k << " not in [0 " << 2i1-1 << "]");
827	next_2[k]=head_2[v];
828	if (v>nbv-1 \|\| v<0) _error_("Index exceed matrix dimensions (v=" << v << " not in [0 " << nbv-1 << "])");
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
838	connectivitymax_2++;//add last column for size
839	bamgmesh->NodalConnectivitySize[0]=nbv;
840	bamgmesh->NodalConnectivitySize[1]=connectivitymax_2;
841	bamgmesh->NodalConnectivity=xNew<double>(connectivitymax_2*nbv);
842	for (i=0;i<connectivitymax_2*nbv;i++) bamgmesh->NodalConnectivity[i]=0;
843	for (i=0;i<nbv;i++){
844	k=0;
845	for(j=head_2[i];j!=-1;j=next_2[j]){
846	_assert_(connectivitymax_2i+k < connectivitymax_2nbv);
847	num=(int)bamgmesh->IssmEdges[int(j/2)*i2+0];
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
850	bamgmesh->NodalConnectivity[connectivitymax_2i+k]=bamgmesh->IssmEdges[int(j/2)i2+1];
851	}
852	else{
853	bamgmesh->NodalConnectivity[connectivitymax_2*i+k]=num;
854	}
855	k++;
856	}
857	bamgmesh->NodalConnectivity[connectivitymax_2*(i+1)-1]=k;
858	}
859
860	/Cracked vertices/
861	if(verbose>5) _printf_(" writing Cracked vertices\n");
862	bamgmesh->CrackedVerticesSize[0]=NbCrackedVertices;
863	bamgmesh->CrackedVerticesSize[1]=2;
864	if (NbCrackedVertices){
865	bamgmesh->CrackedVertices=xNew<double>(2*NbCrackedVertices);
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/
873	if(verbose>5) _printf_(" writing Cracked vertices\n");
874	bamgmesh->CrackedEdgesSize[0]=NbCrackedEdges;
875	bamgmesh->CrackedEdgesSize[1]=4;
876	if (NbCrackedEdges){
877	bamgmesh->CrackedEdges=xNew<double>(2*NbCrackedEdges);
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
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);
891	delete [] reft;
892	delete [] numt;
893	}
894	/}}}/
895	void Mesh::ReadMetric(const BamgOpts* bamgopts) {/{{{/
896
897	/Intermediary/
898	int i,j;
899
900	if(bamgopts->verbose>3) _printf_(" processing metric\n");
901	double hmin = Max(bamgopts->hmin,MinimalHmin());
902	double hmax = Min(bamgopts->hmax,MaximalHmax());
903	double coef = bamgopts->coeff;
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);
922	EigenMetric Vp(M/coef);
923
924	Vp.Maxh(hmax);
925	Vp.Minh(hmin);
926	vertices[i].m = Vp;
927	}
928	}
929	}
930	}
931	/}}}/
932	void Mesh::WriteMetric(BamgOpts* bamgopts) {/{{{/
933	int i;
934	xDelete<double>(bamgopts->metric);
935	bamgopts->metric=xNew<double>(3*nbv);
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	}
942	/}}}/
943	void Mesh::WriteIndex(int** pindex,int* pnels){/{{{/
944
945	/Intermediary/
946	int i,k;
947
948	/output/
949	int* index=NULL;
950	int num=0;
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){
960	index=xNew<int>(3*k);
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) )){
964	/Remove triangles that have a bad aspect ratio/
965	//if(t.Anisotropy()<2 & t.Length()<1.e+5){
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;
970	//}
971	}
972	}
973	}
974
975	/Assign output pointers/
976	*pindex=index;
977	*pnels=num;
978	}
979	/}}}/
980
981	/Methods/
982	void Mesh::AddGeometryMetric(BamgOpts* bamgopts){/{{{/
983	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/IntersectGeomMetric)/
984
985	/Get options/
986	double anisomax = bamgopts->anisomax;
987	double errg = bamgopts->errg;
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
995	_assert_(hmax>0);
996
997	//errC cannot be higher than 1
998	if(errC>1) errC=1;
999
1000	//Set all vertices to "on"
1001	SetVertexFieldOn();
1002
1003	//loop over all the vertices on edges
1004	for(int i=0;i<nbe;i++){
1005	for (int j=0;j<2;j++){
1006
1007	BamgVertex V;
1008	VertexOnGeom GV;
1009	Gh.ProjectOnCurve(edges[i],ss[j],V,GV);
1010
1011	GeomEdge* eg = GV;
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){
1023	_error_("ht<=0 \|\| hn<=0");
1024	}
1025	EigenMetric Vp(1/(htht),1/(hnhn),tg);
1026	Metric MVp(Vp);
1027	edges[i][j].m.IntersectWith(MVp);
1028	}
1029	}
1030	// the problem is for the vertex on vertex
1031	}
1032	/}}}/
1033	void Mesh::AddMetric(BamgOpts* bamgopts){/{{{/
1034	// Hessiantype = 0 => H is computed using double L2 projection
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{
1047	_error_("Hessiantype " << Hessiantype << " not supported yet (1->use Green formula, 0-> double L2 projection)");
1048	}
1049	}
1050	/}}}/
1051	void Mesh::AddVertex( BamgVertex &s,Triangle* t, Icoor2* det3){/{{{/
1052	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Add)/
1053	// -------------------------------
1054	// s2
1055	// !
1056	// /\|\ !
1057	// / \| \ !
1058	// / \| \ !
1059	// tt1 / \| \ tt0 !
1060	// / \|s \ !
1061	// / . \ !
1062	// / . ` \ !
1063	// / . ` \ !
1064	// ---------------- !
1065	// s0 tt2 s1
1066	//-------------------------------
1067
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/
1076	BamgVertex* s0=t->GetVertex(0);
1077	BamgVertex* s1=t->GetVertex(1);
1078	BamgVertex* s2=t->GetVertex(2);
1079
1080	//determinant of t
1081	Icoor2 detOld=t->det;
1082
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*/
1086	int infvertexindex = s0 ? ((s1? (s2?-1:2):1)):0;
1087
1088	//some checks
1089	if(((infvertexindex <0 ) && (detOld <0)) \|\| ((infvertexindex >=0) && (detOld >0)) ){
1090	_error_("inconsistent configuration (Contact ISSM developers)");
1091	}
1092
1093	// if det3 does not exist, build it
1094	if (!det3){
1095	//allocate
1096	det3 = det3local;
1097	//if no infinite vertex
1098	if (infvertexindex<0 ) {
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());}
1102	else {
1103	// one of &s1 &s2 &s0 is NULL
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()) ;
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
1115	if (nbzerodet>0){
1116	/s is on an edge/
1117	if (nbzerodet==1) {
1118	iedge = OppositeEdge[izerodet];
1119	AdjacentTriangle ta = t->Adj(iedge);
1120
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 ) {
1125	// add in outside triangle
1126	AddVertex(s,( Triangle *)ta);
1127	return;
1128	}
1129	}
1130	}
1131	else{
1132	_error_("Cannot add a vertex more than once. Check duplicates");
1133	}
1134	}
1135
1136	// remove de MarkUnSwap edge
1137	t->SetUnMarkUnSwap(0);
1138	t->SetUnMarkUnSwap(1);
1139	t->SetUnMarkUnSwap(2);
1140
1141	tt[0]= t;
1142	tt[1]= &triangles[nbt++];
1143	tt[2]= &triangles[nbt++];
1144
1145	if (nbt>maxnbt) _error_("Not enough triangles");
1146
1147	tt[1]=tt[2]=*t;
1148	tt[0]->link=tt[1];
1149	tt[1]->link=tt[2];
1150
1151	(*tt[0])(OppositeVertex[0])=&s;
1152	(*tt[1])(OppositeVertex[1])=&s;
1153	(*tt[2])(OppositeVertex[2])=&s;
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
1172	tt[0]->SetSingleVertexToTriangleConnectivity();
1173	tt[1]->SetSingleVertexToTriangleConnectivity();
1174	tt[2]->SetSingleVertexToTriangleConnectivity();
1175
1176	// swap if the point s is on a edge
1177	if(izerodet>=0) {
1178	int rswap=tt[izerodet]->swap(iedge);
1179
1180	if (!rswap) {
1181	_error_("swap the point s is on a edge");
1182	}
1183	}
1184
1185	}/}}}/
1186	void Mesh::BoundAnisotropy(double anisomax,double hminaniso) {/{{{/
1187	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/BoundAnisotropy)/
1188
1189	long int verbose=0;
1190	double lminaniso = 1/ (Max(hminaniso*hminaniso,1e-100));
1191
1192	//display info
1193	if (verbose > 1) _printf_(" BoundAnisotropy by " << anisomax << "\n");
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++){
1201	EigenMetric Vp(vertices[i]);
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
1208	if (verbose>2){
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
1219	if (verbose>2){
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");
1222	}
1223	}
1224	/}}}/
1225	void Mesh::BuildGeometryFromMesh(BamgOpts* bamgopts){/{{{/
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;
1232	int verbose=0;
1233	double cutoffradian=10*Pi/180;
1234
1235	/Recover options/
1236	if (bamgopts){
1237	verbose=bamgopts->verbose;
1238	cutoffradian=bamgopts->MaxCornerAngle*Pi/180;
1239	}
1240
1241	//display info
1242	if (verbose>1) _printf_(" construction of the geometry from the 2d mesh\n");
1243
1244	//check that the mesh is not empty
1245	if (nbt<=0 \|\| nbv <=0 ) {
1246	_error_("nbt or nbv is negative (Mesh empty?)");
1247	}
1248
1249	//Gh is the geometry of the mesh (this), initialize MaxCornerAngle
1250	if (cutoffradian>=0) Gh.MaxCornerAngle = cutoffradian;
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++){
1263	edge4->SortAndAdd(GetId(edges[i][0]),GetId(edges[i][1]));
1264	}
1265	//check that there is no double edge
1266	if (nbe != edge4->nb()){
1267	delete [] st;
1268	_error_("Some Double edge in the mesh, the number is " << nbe << ", nbe4=" << edge4->nb());
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)
1278	long k =edge4->SortAndAdd(GetId(triangles[i][VerticesOfTriangularEdge[j][0]]), GetId(triangles[i][VerticesOfTriangularEdge[j][1]]));
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))){
1288	_error_("problem in Geometry reconstruction: an edge on boundary is duplicated (double element?)");
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 {
1302	_printf_("The edge (" << GetId(triangles[i][VerticesOfTriangularEdge[j][0]]) << "," << GetId(triangles[i][VerticesOfTriangularEdge[j][1]]) << ") belongs to more than 2 triangles (" << k << ")\n");
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");
1306	_error_("An edge belongs to more than 2 triangles");
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
1316	if(verbose>5) {
1317	_printf_(" info on Mesh:\n");
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");
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
1353	if(verbose>4) _printf_(" Construction of the edges " << nbe << "\n");
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]];
1374	edges[add].GeomEdgeHook=NULL;
1375	//if already existed
1376	if (i<nbeold){
1377	edges[add].ReferenceNumber=edgessave[i].ReferenceNumber;
1378	edges[add].GeomEdgeHook=edgessave[i].GeomEdgeHook; // HACK to get required edges
1379	_printf_("oh no...\n");
1380	}
1381	else
1382	edges[add].ReferenceNumber=Min(edges[add].v[0]->GetReferenceNumber(),edges[add].v[1]->GetReferenceNumber());
1383	}
1384	}
1385
1386	//check that we have been through all edges
1387	if (k!=nbe){
1388	_error_("problem in edge construction process: k!=nbe (should not happen)");
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
1414	BamgVertex* v=edges[i].v[j];
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]){
1433	_error_("v!=edges[i0 ].v[j0]: this should not happen as the vertex belongs to this edge");
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
1448	//check that nbsubdomains is empty
1449	if (nbsubdomains){
1450	_error_("nbsubdomains should be 0");
1451	}
1452	nbsubdomains=0;
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
1468	colorT[it]=nbsubdomains;
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
1482	if ( ! t->Locked(j) && tt && (colorT[jt = GetId(tt)] == -1) && ( tt->color==kolor)) {
1483	colorT[jt]=nbsubdomains;
1484	st[++level]=jt;
1485	st[++level]=0;
1486	k++;
1487	}
1488	}
1489	else level-=2;
1490	}
1491	nbsubdomains++;
1492	}
1493	}
1494	if (verbose> 3) _printf_(" The Number of sub domain = " << nbsubdomains << "\n");
1495
1496	//build subdomains
1497	long isd;
1498	subdomains = new SubDomain[nbsubdomains];
1499
1500	//initialize subdomains[isd].head as 0
1501	for (isd=0;isd<nbsubdomains;isd++) subdomains[isd].head =0;
1502
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;
1509	subdomains[isd].ReferenceNumber = triangles[it].color;
1510	subdomains[isd].direction = j; // hack
1511	subdomains[isd].edge = 0;
1512	k++;
1513	}
1514	}
1515	}
1516	//check that we have been through all subdomains
1517	if (k!= nbsubdomains){
1518	delete [] colorT;
1519	_error_("k!= nbsubdomains");
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++){
1531	for ( j=0;j<2;j++) colorV[GetId(edges[i][j])]=0;
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;
1542	Gh.vertices = new GeomVertex[k];
1543	Gh.edges = new GeomEdge[nbe];
1544	Gh.nbsubdomains = nbsubdomains;
1545	Gh.subdomains = new GeomSubDomain[nbsubdomains];
1546	if (verbose>3) _printf_(" number of vertices = " << Gh.nbv << "\n number of edges = " << Gh.nbe << "\n");
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){
1555	BamgVertex & v = Gh.vertices[j];
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
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));
1579	if (Gh.coefIcoor<=0){
1580	delete [] colorV;
1581	_error_("Gh.coefIcoor<=0 in infered Geometry (this should not happen)");
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
1596	long i0 = GetId(edges[i][0]);
1597	long i1 = GetId(edges[i][1]);
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
1604	Gh.edges[i].type = 0;
1605
1606	Gh.edges[i].tg[0]=R2();
1607	Gh.edges[i].tg[1]=R2();
1608
1609	bool required= edges[i].GeomEdgeHook;
1610	if(required) kreq++;
1611	edges[i].GeomEdgeHook = Gh.edges + i;
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);
1628	Gh.edges[i].ReferenceNumber = edges[i].ReferenceNumber;
1629
1630	k = edge4->SortAndAdd(i0,i1);
1631	if (k != i){
1632	delete [] len;
1633	delete [] colorV;
1634	_error_("problem in Edge4 construction: k != i");
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
1649	for (i=0;i<nbsubdomains;i++){
1650	it = GetId(subdomains[i].head);
1651	j = subdomains[i].direction;
1652	long i0 = GetId(triangles[it][VerticesOfTriangularEdge[j][0]]);
1653	long i1 = GetId(triangles[it][VerticesOfTriangularEdge[j][1]]);
1654	k = edge4->SortAndFind(i0,i1);
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	}
1662
1663	delete edge4;
1664	delete [] colorV;
1665
1666	//unset adj
1667	for (i=0;i<nbt;i++){
1668	for (j=0;j<3;j++){
1669	triangles[i].SetAdj2(j,0,triangles[i].GetAllflag(j));
1670	}
1671	}
1672
1673	}
1674	/}}}/
1675	void Mesh::BuildMetric0(BamgOpts* bamgopts){/{{{/
1676
1677	/Options/
1678	double* s=NULL;
1679	long nbsol;
1680	int verbose;
1681
1682	int i,j,k,iA,iB,iC;
1683	int iv;
1684
1685	/Recover options/
1686	verbose=bamgopts->verbose;
1687
1688	/Get and process fields/
1689	s=bamgopts->field;
1690	nbsol=bamgopts->fieldSize[1];
1691
1692	/Check size/
1693	if (bamgopts->fieldSize[0] != nbv) _error_("'field' should have " << nbv << " rows");
1694
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
1714	if(verbose>1) {
1715	_printf_(" Construction of Metric: number of field: " << nbsol << " (nbt=" << nbt << ", nbv=" << nbv << ")\n");
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;
1720	head_s=xNew<int>(nbv);
1721	int* next_p=NULL;
1722	next_p=xNew<int>(3*nbt);
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++){
1763	k=GetId(triangles[i][j]);
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
1787	if(verbose>2) _printf_(" Solution " << nusol << ", Min = " << smin << ", Max = " << smax << ", Delta = " << sdelta << "\n");
1788
1789	//skip constant field
1790	if (sdelta < 1.0e-10*Max(absmax,1e-20)){
1791	_printf_(" Solution " << nusol << " is constant, skipping...\n");
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
1802	iA = GetId(triangles[i][0]);
1803	iB = GetId(triangles[i][1]);
1804	iC = GetId(triangles[i][2]);
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
1831	iA = GetId(triangles[i][0]);
1832	iB = GetId(triangles[i][1]);
1833	iC = GetId(triangles[i][2]);
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
1861	xDelete<int>(head_s);
1862	xDelete<int>(next_p);
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	}
1878	/}}}/
1879	void Mesh::BuildMetric1(BamgOpts* bamgopts){/{{{/
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;
1886	int verbose;
1887
1888	/Recover options/
1889	verbose=bamgopts->verbose;
1890	NbJacobi=bamgopts->nbjacobi;
1891
1892	/Get and process fields/
1893	s=bamgopts->field;
1894	nbsol=bamgopts->fieldSize[1];
1895
1896	/Check size/
1897	if (bamgopts->fieldSize[0] != nbv) _error_("'field' should have " << nbv << " rows");
1898
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
1915	if(verbose>1) {
1916	_printf_(" Construction of Metric: number of field: " << nbsol << " (nbt=" << nbt << ", nbv=" << nbv << ")\n");
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
1941	iA = GetId(t[0]);
1942	iB = GetId(t[1]);
1943	iC = GetId(t[2]);
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
1958	OnBoundary[GetId(t[VerticesOfTriangularEdge[j][0]])]=1;
1959	OnBoundary[GetId(t[VerticesOfTriangularEdge[j][1]])]=1;
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
1999	if(verbose>2) _printf_(" Solution " << nusol << ", Min = " << smin << ", Max = " << smax << ", Delta = " << sdelta << ", number of fields = " << nbsol << "\n");
2000
2001	//skip constant field
2002	if (sdelta < 1.0e-10*Max(absmax,1e-20) ){
2003	if (verbose>2) _printf_(" Solution " << nusol << " is constant, skipping...\n");
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
2031	iA = GetId(triangles[i][0]);
2032	iB = GetId(triangles[i][1]);
2033	iC = GetId(triangles[i][2]);
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]);
2098	EigenMetric Vp(M);
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++) {
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");
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
2121	iA = GetId(triangles[i][0]);
2122	iB = GetId(triangles[i][1]);
2123	iC = GetId(triangles[i][2]);
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
2134	iA = GetId(triangles[i][0]);
2135	iB = GetId(triangles[i][1]);
2136	iC = GetId(triangles[i][2]);
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	}
2170	/}}}/
2171	void Mesh::CrackMesh(BamgOpts* bamgopts) {/{{{/
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++){
2184	if(edges[i].GeomEdgeHook->Cracked()) k++;
2185	}
2186
2187	//Return if no edge is cracked
2188	if(k==0) return;
2189	if (verbose>4) _printf_(" number of Cracked Edges = " << k << "\n");
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++){
2203	if(edges[i].GeomEdgeHook->Cracked()){
2204
2205	//Fill edges fields of CrackedEdges
2206	CrackedEdges[k ].E =edges[i].GeomEdgeHook;
2207	CrackedEdges[k++].e1=&edges[i];
2208
2209	//Get number of the two vertices on the edge
2210	i1=GetId(edges[i][0]);
2211	i2=GetId(edges[i][1]);
2212	_assert_(i1>=0 && i1<nbv && i2>=0 && i2<nbv);
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){
2222	delete [] splitvertex;
2223	_error_("Crossing rifts not supported yet");
2224	}
2225	}
2226	}
2227	_assert_(k==NbCrackedEdges);
2228
2229	//Add new vertices
2230	if (verbose>4) _printf_(" number of Cracked Vertices = " << NbCrackedVertices << "\n");
2231	if (NbCrackedVertices){
2232	CrackedVertices=xNew<long>(2*NbCrackedVertices);
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	}
2241	_assert_(num==NbCrackedVertices);
2242	}
2243	delete [] splitvertex;
2244
2245	//Now, find the triangles that hold a cracked edge
2246	CreateSingleVertexToTriangleConnectivity();
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
2253	i1=GetId((*CrackedEdges[i].e1)[0]);
2254	i2=GetId((*CrackedEdges[i].e1)[1]);
2255
2256	//find a triangle holding the vertex i1 (first vertex of the ith cracked edge)
2257	Triangle* tbegin=vertices[i1].t;
2258	k=vertices[i1].IndexInTriangle;//local number of i in triangle tbegin
2259	_assert_(GetId((*tbegin)[k])==GetId(vertices[i1]));
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
2263	AdjacentTriangle ta(tbegin,EdgesVertexTriangle[k][0]);
2264	count=0;
2265	do {
2266	for(j=0;j<3;j++){
2267	//Find the position of i1 in the triangle index
2268	if (GetId((*ta.t)[j])==i1){
2269	j1=j;
2270	break;
2271	}
2272	}
2273	for(j=0;j<3;j++){
2274	//Check wether i2 is also in the triangle index
2275	if (GetId((*ta.t)[j])==i2){
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	}
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");
2300	ta = Next(ta).Adj();
2301	if (count++>50) _error_("Maximum number of iteration exceeded");
2302	}while ((tbegin != ta));
2303	}
2304
2305	//Check EdgeFlag
2306	for(i=0;i<NbCrackedEdges;i++){
2307	if (Edgeflags[i]!=2){
2308	_error_("A problem occured: at least one crack edge (number " << i+1 << ") does not belong to 2 elements");
2309	}
2310	}
2311	delete [] Edgeflags;
2312
2313	//Reset BamgVertex to On
2314	SetVertexFieldOn();
2315
2316	}
2317	/}}}/
2318	void Mesh::Echo(void) {/{{{/
2319
2320	int i;
2321
2322	_printf_("Mesh Echo:\n");
2323	_printf_(" nbv = " << nbv << "\n");
2324	_printf_(" nbt = " << nbt << "\n");
2325	_printf_(" nbe = " << nbe << "\n");
2326	_printf_(" index:\n");
2327	for (i=0;i<nbt;i++){
2328	_printf_(" " << setw(4) << i+1 << ": ["
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");
2332	}
2333	_printf_(" coordinates:\n");
2334	for (i=0;i<nbv;i++){
2335	_printf_(" " << setw(4) << i+1 << ": [" << vertices[i].r.x << " " << vertices[i].r.y << "] and [" << vertices[i].i.x << " " << vertices[i].i.y << "]\n");
2336	}
2337
2338	}
2339	/}}}/
2340	void Mesh::ForceBoundary() {/{{{/
2341	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ForceBoundary)/
2342
2343	long int verbose=2;
2344	int k=0;
2345	int nbfe=0,nbswp=0,Nbswap=0;
2346
2347	//display
2348	if (verbose > 2) _printf_(" ForceBoundary nb of edge: " << nbe << "\n");
2349
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	}
2357
2358	//Force Edges
2359	AdjacentTriangle ta(0,0);
2360	for (int i = 0; i < nbe; i++){
2361
2362	//Force edge i
2363	nbswp = ForceEdge(edges[i][0],edges[i][1],ta);
2364	if (nbswp<0) k++;
2365	else Nbswap += nbswp;
2366
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]));
2370	}
2371	}
2372
2373	if (k!=0) {
2374	_error_("There are " << k << " lost edges, the boundary might be crossing");
2375	}
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	}
2381	/}}}/
2382	void Mesh::FindSubDomain(int OutSide) {/{{{/
2383	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FindSubDomain)/
2384
2385	long int verbose=0;
2386
2387	if (verbose >2){
2388	if (OutSide) _printf_(" Find all external sub-domain\n");
2389	else _printf_(" Find all internal sub-domain\n");
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;
2434	nbtout = 0;
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];
2443	nbtout--; // on fait un coup de trop.
2444	while (t){
2445	nbtout++;
2446	t1=t;
2447	t=t->link;
2448	t1->link=0;
2449	}
2450	}
2451	}
2452	it++;} // end while (it<nbt)
2453	if (nbt == nbtout \|\| !NbSubDomTot) {
2454	delete [] HeapArete;
2455	_error_("The boundary is not close: all triangles are outside");
2456	}
2457
2458	delete [] HeapArete;
2459	delete [] HeapTriangle;
2460
2461	if (OutSide\|\| !Gh.subdomains \|\| !Gh.nbsubdomains )
2462	{ // No geom sub domain
2463	long i;
2464	if (subdomains) delete [] subdomains;
2465	subdomains = new SubDomain[ NbSubDomTot];
2466	nbsubdomains= NbSubDomTot;
2467	for ( i=0;i<nbsubdomains;i++) {
2468	subdomains[i].head=NULL;
2469	subdomains[i].ReferenceNumber=i+1;
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 {
2484	mark[GetId(t)]=k;
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)
2490	if (k!=nbsubdomains){
2491	delete [] mark;
2492	_error_("k!=nbsubdomains");
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
2499	long nbk = nbsubdomains;
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);
2505	long kl = ta ? mark[GetId(ta)] : -2;
2506	long kr = mark[it];
2507	if(kr !=kl) {
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;
2516	}
2517	}
2518	long j=0;
2519	for ( i=0;i<nbsubdomains;i++)
2520	if((-subdomains[i].ReferenceNumber) %2) { // good
2521	if(i != j)
2522	Exchange(subdomains[i],subdomains[j]);
2523	j++;}
2524	else{
2525	t= subdomains[i].head;
2526	while (t){
2527	nbtout++;
2528	t1=t;
2529	t=t->link;
2530	t1->link=0;
2531	}//while (t)
2532	}
2533	if(verbose>4) _printf_(" Number of removes subdomains (OutSideMesh) = " << nbsubdomains-j << "\n");
2534	nbsubdomains=j;
2535	}
2536
2537	delete [] mark;
2538
2539	}
2540	else{ // find the head for all subdomains
2541	if (Gh.nbsubdomains != nbsubdomains && subdomains)
2542	delete [] subdomains, subdomains=0;
2543	if (! subdomains )
2544	subdomains = new SubDomain[ Gh.nbsubdomains];
2545	nbsubdomains =Gh.nbsubdomains;
2546	CreateSingleVertexToTriangleConnectivity();
2547	long * mark = new long[nbt];
2548	Edge **GeomEdgetoEdge = MakeGeomEdgeToEdge();
2549
2550	for (it=0;it<nbt;it++)
2551	mark[it]=triangles[it].link ? -1 : -2;
2552	long inew =0;
2553	for (int i=0;i<nbsubdomains;i++) {
2554	GeomEdge &eg = *Gh.subdomains[i].edge;
2555	subdomains[i].ReferenceNumber = Gh.subdomains[i].ReferenceNumber;
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)];
2559	_assert_(&e);
2560	BamgVertex * v0 = e(0),*v1 = e(1);
2561	Triangle *t = v0->t;
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;
2566	subdomains[i].edge = &e;
2567	_assert_(t && direction);
2568
2569	AdjacentTriangle ta(t,EdgesVertexTriangle[v0->IndexInTriangle][0]);// previous edges
2570
2571	while (1) {
2572	_assert_(v0==ta.EdgeVertex(1));
2573	if (ta.EdgeVertex(0) == v1) { // ok we find the edge
2574	if (direction>0)
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) {
2579	_error_("bad definition of SubSomain " << i);
2580	}
2581	long it = GetId(t);
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++;
2593	if (mark[GetId(tt)]>=0){
2594	_error_("mark[GetId(tt)]>=0");
2595	}
2596	mark[GetId(tt)]=i;
2597	tt=tt->link;
2598	} while (tt!=t);
2599	break;
2600	}
2601	ta = Previous(Adj(ta));
2602	if(t == (Triangle *) ta) {
2603	_error_("bad definition of SubSomain " << i);
2604	}
2605	}
2606	}
2607
2608	if (inew < nbsubdomains) {
2609	if (verbose>5) _printf_("WARNING: " << nbsubdomains-inew << " SubDomains are being removed\n");
2610	nbsubdomains=inew;}
2611
2612	for (it=0;it<nbt;it++)
2613	if ( mark[it] ==-1 )
2614	nbtout++,triangles[it].link =0;
2615	delete [] GeomEdgetoEdge;
2616	delete [] mark;
2617
2618	}
2619	nbtout=0;
2620	for (it=0;it<nbt;it++)
2621	if(!triangles[it].link) nbtout++;
2622	}
2623	/}}}/
2624	long Mesh::GetId(const Triangle & t) const { /{{{/
2625	return &t - triangles;
2626	}
2627	/}}}/
2628	long Mesh::GetId(const Triangle * t) const { /{{{/
2629	return t - triangles;
2630	}
2631	/}}}/
2632	long Mesh::GetId(const BamgVertex & t) const { /{{{/
2633	return &t - vertices;
2634	}
2635	/}}}/
2636	long Mesh::GetId(const BamgVertex * t) const { /{{{/
2637	return t - vertices;
2638	}
2639	/}}}/
2640	long Mesh::GetId(const Edge & t) const { /{{{/
2641	return &t - edges;
2642	}
2643	/}}}/
2644	long Mesh::GetId(const Edge * t) const { /{{{/
2645	return t - edges;
2646	}
2647	/}}}/
2648	void Mesh::Init(long maxnbv_in) {/{{{/
2649	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/PreInit)/
2650
2651	/* initialize random seed: */
2652	srand(19999999);
2653
2654	/Initialize fields/
2655	NbRef=0;
2656	quadtree=NULL;
2657	nbv=0;
2658	nbt=0;
2659	nbe=0;
2660	edges=NULL;
2661	nbsubdomains=0;
2662	subdomains=NULL;
2663	maxnbv=maxnbv_in;
2664	maxnbt=2 *maxnbv_in-2;
2665	NbVertexOnBThVertex=0;
2666	VertexOnBThVertex=NULL;
2667	NbVertexOnBThEdge=0;
2668	VertexOnBThEdge=NULL;
2669	NbCrackedVertices=0;
2670	CrackedVertices =NULL;
2671	NbCrackedEdges =0;
2672	CrackedEdges =NULL;
2673	NbVerticesOnGeomVertex=0;
2674	VerticesOnGeomVertex=NULL;
2675	NbVerticesOnGeomEdge=0;
2676	VerticesOnGeomEdge=NULL;
2677
2678	/Allocate if maxnbv_in>0/
2679	if(maxnbv_in){
2680	vertices=new BamgVertex[maxnbv];
2681	_assert_(vertices);
2682	orderedvertices=new BamgVertex* [maxnbv];
2683	_assert_(orderedvertices);
2684	triangles=new Triangle[maxnbt];
2685	_assert_(triangles);
2686	}
2687	else{
2688	vertices=NULL;
2689	orderedvertices=NULL;
2690	triangles=NULL;
2691	maxnbt=0;
2692	}
2693	}
2694	/}}}/
2695	void Mesh::Insert(bool random) {/{{{/
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/
2704	long int verbose=2;
2705
2706	//Display info
2707	if (verbose>2) _printf_(" Insert initial " << nbv << " vertices\n");
2708
2709	//Compute integer coordinates for the existing vertices
2710	SetIntCoor();
2711
2712	/*Now we want to build a list (orderedvertices) of the vertices in a random
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	*
2723	* for all j in [0 nbv[, there is a unique i in [0 nbv[ such that k_i=j
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*/
2739
2740	//Get Prime number
2741	const long PrimeNumber= BigPrimeNumber(nbv);
2742	int temp = rand(); if(!random) temp = 756804110;
2743	int k0=temp%nbv;
2744
2745	//Build orderedvertices
2746	for (i=0; i<nbv; i++){
2747	orderedvertices[i]=&vertices[k0=(k0+PrimeNumber)%nbv];
2748	}
2749
2750	/Modify orderedvertices such that the first 3 vertices form a triangle/
2751
2752	//get first vertex i such that [0,1,i] are not aligned
2753	for (i=2; det(orderedvertices[0]->i,orderedvertices[1]->i,orderedvertices[i]->i)==0;){
2754	//if i is higher than nbv, it means that all the determinants are 0,
2755	//all vertices are aligned!
2756	if (++i>=nbv) _error_("all the vertices are aligned");
2757	}
2758	// exchange i et 2 in "orderedvertices" so that
2759	// the first 3 vertices are not aligned (real triangle)
2760	Exchange(orderedvertices[2], orderedvertices[i]);
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
2765	BamgVertex v0=orderedvertices[0], v1=orderedvertices[1];
2766
2767	nbt = 2;
2768	triangles[0](0) = NULL;//infinite vertex
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
2786	triangles[0].SetSingleVertexToTriangleConnectivity();
2787	triangles[1].SetSingleVertexToTriangleConnectivity();
2788
2789	triangles[0].link=&triangles[1];
2790	triangles[1].link=&triangles[0];
2791
2792	//build quadtree
2793	if (!quadtree) quadtree = new BamgQuadtree(this,0);
2794	quadtree->Add(*v0);
2795	quadtree->Add(*v1);
2796
2797	/Now, add the vertices One by One/
2798	long NbSwap=0;
2799	if (verbose>3) _printf_(" Begining of insertion process...\n");
2800
2801	for (int icount=2; icount<nbv; icount++) {
2802
2803	//Get new vertex
2804	BamgVertex *newvertex=orderedvertices[icount];
2805
2806	//Find the triangle in which newvertex is located
2807	Icoor2 det3[3];
2808	Triangle* tcvi = TriangleFindFromCoord(newvertex->i,det3); //(newvertex->i = integer coordinates)
2809
2810	//Add newvertex to the quadtree
2811	quadtree->Add(*newvertex);
2812
2813	//Add newvertex to the existing mesh
2814	AddVertex(*newvertex,tcvi,det3);
2815
2816	//Make the mesh Delaunay around newvertex by swaping the edges
2817	NbSwap += newvertex->Optim(1,0);
2818	}
2819
2820	//Display info
2821	if (verbose>3) {
2822	_printf_(" NbSwap of insertion: " << NbSwap << "\n");
2823	_printf_(" NbSwap/nbv: " << NbSwap/nbv << "\n");
2824	}
2825	}
2826	/}}}/
2827	long Mesh::InsertNewPoints(long nbvold,long & NbTSwap,bool random) {/{{{/
2828	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/InsertNewPoints)/
2829
2830	long int verbose=0;
2831	double seuil= 1.414/2 ;// for two close point
2832	long i;
2833	long NbSwap=0;
2834	Icoor2 det3[3];
2835
2836	//number of new points
2837	const long nbvnew=nbv-nbvold;
2838
2839	//display info if required
2840	if (verbose>5) _printf_(" Try to Insert " << nbvnew << " new points\n");
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
2848	int temp = rand(); if(!random) temp = 1098566905;
2849	long k3 = temp%nbvnew;
2850	//loop over the new points
2851	for (int is3=0; is3<nbvnew; is3++){
2852	long j=nbvold +(k3 = (k3+PrimeNumber)%nbvnew);
2853	long i=nbvold+is3;
2854	orderedvertices[i]= vertices + j;
2855	orderedvertices[i]->ReferenceNumber=i;
2856	}
2857
2858	// for all the new point
2859	long iv=nbvold;
2860	for(i=nbvold;i<nbv;i++){
2861	BamgVertex &vi=*orderedvertices[i];
2862	vi.i=R2ToI2(vi.r);
2863	vi.r=I2ToR2(vi.i);
2864	double hx,hy;
2865	vi.m.Box(hx,hy);
2866	int hi=(int) (hxcoefIcoor),hj=(int) (hycoefIcoor);
2867	if(!quadtree->TooClose(&vi,seuil,hi,hj)){
2868	// a good new point
2869	BamgVertex &vj = vertices[iv];
2870	long j=vj.ReferenceNumber;
2871	if (&vj!=orderedvertices[j]){
2872	_error_("&vj!= orderedvertices[j]");
2873	}
2874	if(i!=j){
2875	Exchange(vi,vj);
2876	Exchange(orderedvertices[j],orderedvertices[i]);
2877	}
2878	vj.ReferenceNumber=0;
2879	Triangle *tcvj=TriangleFindFromCoord(vj.i,det3);
2880	if (tcvj && !tcvj->link){
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");
2884	tcvj->Echo();
2885	_error_("problem inserting point in InsertNewPoints (tcvj=" << tcvj << " and tcvj->link=" << tcvj->link << ")");
2886	}
2887	quadtree->Add(vj);
2888	AddVertex(vj,tcvj,det3);
2889	NbSwap += vj.Optim(1);
2890	iv++;
2891	}
2892	else{
2893	vi.PreviousNumber = 0;
2894	}
2895	}
2896	if (verbose>3) {
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");
2900	}
2901	nbv = iv;
2902
2903	for (i=nbvold;i<nbv;i++) NbSwap += vertices[i].Optim(1);
2904	if (verbose>3) _printf_(" NbSwap=" << NbSwap << "\n");
2905
2906	NbTSwap += NbSwap ;
2907	return nbv-nbvold;
2908	}
2909	/}}}/
2910	Edge** Mesh::MakeGeomEdgeToEdge() {/{{{/
2911	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeGeomEdgeToEdge)/
2912
2913	if (!Gh.nbe){
2914	_error_("!Gh.nbe");
2915	}
2916	Edge *e= new Edge [Gh.nbe];
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;
2924	GeomEdge *GeomEdgeHook = ei->GeomEdgeHook;
2925	e[Gh.GetId(GeomEdgeHook)] = ei;
2926	}
2927	for ( i=0;i<nbe ; i++)
2928	for (int ii=0;ii<2;ii++) {
2929	Edge * ei = edges+i;
2930	GeomEdge *GeomEdgeHook = ei->GeomEdgeHook;
2931	int j= ii;
2932	while (!(*GeomEdgeHook)[j].Required()) {
2933	Adj(GeomEdgeHook,j); // next geom edge
2934	j=1-j;
2935	if (e[Gh.GetId(GeomEdgeHook)]) break; // optimisation
2936	e[Gh.GetId(GeomEdgeHook)] = ei;
2937	}
2938	}
2939
2940	int kk=0;
2941	for ( i=0;i<Gh.nbe ; i++){
2942	if (!e[i]){
2943	kk++;
2944	if(kk<10) _printf_("BUG: the geometrical edge " << i << " is on no edge curve\n");
2945	}
2946	}
2947	if(kk) _error_("See above");
2948
2949	return e;
2950	}
2951	/}}}/
2952	void Mesh::MakeBamgQuadtree() { /{{{/
2953	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeBamgQuadtree)/
2954	if(!quadtree) quadtree = new BamgQuadtree(this);
2955	}
2956	/}}}/
2957	double Mesh::MaximalHmax() {/{{{/
2958	return Max(pmax.x-pmin.x,pmax.y-pmin.y);
2959	}
2960	/}}}/
2961	void Mesh::MaxSubDivision(double maxsubdiv) {/{{{/
2962	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/MaxSubDivision)/
2963
2964	long int verbose=0;
2965
2966	const double maxsubdiv2 = maxsubdiv*maxsubdiv;
2967	if(verbose>1) _printf_(" Limit the subdivision of a edges in the new mesh by " << maxsubdiv << "\n");
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++){
2975	Triangle *ptt=t.TriangleAdj(j);
2976	Triangle &tt = *ptt;
2977	if ( (!ptt \|\| it < GetId(tt)) && ( tt.link \|\| t.link)){
2978	BamgVertex &v0 = t[VerticesOfTriangularEdge[j][0]];
2979	BamgVertex &v1 = t[VerticesOfTriangularEdge[j][1]];
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){
2985	R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M
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){
2998	R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M
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	}
3010	if(verbose>3){
3011	_printf_(" number of metric changes = " << nbchange << ", maximum number of subdivision of a edges before change = " << pow(lmax,0.5) << "\n");
3012	}
3013	}
3014	/}}}/
3015	Metric Mesh::MetricAt(const R2 & A) const { /{{{/
3016	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MetricAt)/
3017
3018	I2 a = R2ToI2(A);
3019	Icoor2 deta[3];
3020	Triangle * t =TriangleFindFromCoord(a,deta);
3021	if (t->det <0) { // outside
3022	double ba,bb;
3023	AdjacentTriangle edge= CloseBoundaryEdge(a,t,ba,bb) ;
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	}
3034	/}}}/
3035	double Mesh::MinimalHmin() {/{{{/
3036	return 2.0/coefIcoor;
3037	}
3038	/}}}/
3039	BamgVertex* Mesh::NearestVertex(int i,int j) {/{{{/
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){/{{{/
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;
3053	long* first_np_or_next_t=new long[maxnbt];
3054	Triangle* t=NULL;
3055
3056	/Recover options/
3057	int verbose=bamgopts->verbose;
3058
3059	/First, insert old points if requested/
3060	if(KeepVertices && (&Bh != this) && (nbv+Bh.nbv< maxnbv)){
3061	if (verbose>5) _printf_(" Inserting initial mesh points\n");
3062	bool pointsoutside = false;
3063	for(i=0;i<Bh.nbv;i++){
3064	BamgVertex &bv=Bh[i];
3065	/Do not insert if the point is outside/
3066	long long det3[3];
3067	Triangle* tcvj=TriangleFindFromCoord(bv.i,det3);
3068	if(tcvj->det<0 \|\| !tcvj->link){
3069	pointsoutside = true;
3070	continue;
3071	}
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){
3083	vertices[nbv].r = bv.r;
3084	vertices[nbv].PreviousNumber = i+1;
3085	vertices[nbv++].m = bv.m;
3086	}
3087	}
3088	if(pointsoutside) _printf_("WARNING: One or more points of the initial mesh fall outside of the geometric boundary\n");
3089	Bh.CreateSingleVertexToTriangleConnectivity();
3090	InsertNewPoints(nbvold,NbTSwap,bamgopts->random);
3091	}
3092	else Bh.CreateSingleVertexToTriangleConnectivity();
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;
3100	if (verbose>5) _printf_(" Big loop\n");
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 ){
3114	_error_("Index problem in NewPoints (i=" << i << " not in [0 " << nbt-1 << "])");
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++){
3121	AdjacentTriangle tj(t,j);
3122	BamgVertex &vA = *tj.EdgeVertex(0);
3123	BamgVertex &vB = *tj.EdgeVertex(1);
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
3130	AdjacentTriangle tadjj = t->Adj(j);
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;
3138	k=GetId(ta);
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);
3144	lIntTria.NewPoints(vertices,nbv,maxnbv);
3145	} // end loop for each edge
3146	}// for triangle
3147
3148	if (!InsertNewPoints(nbvold,NbTSwap,bamgopts->random)) break;
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++){
3154	BamgVertex* s = vertices + i;
3155	AdjacentTriangle ta(s->t, EdgesVertexTriangle[s->IndexInTriangle][1]);
3156	Triangle* tbegin= (Triangle*) ta;
3157	long kt;
3158	do {
3159	kt = GetId((Triangle*) ta);
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){
3165	_error_("ta.EdgeVertex(0)!=s");
3166	}
3167	ta = Next(Adj(ta));
3168	} while ( (tbegin != (Triangle*) ta));
3169	}
3170
3171	}while(nbv!=nbvold);
3172	delete [] first_np_or_next_t;
3173
3174	long NbSwapf =0;
3175	for(i=0;i<nbv;i++) NbSwapf += vertices[i].Optim(0);
3176	}/}}}/
3177	GeomEdge* Mesh::ProjectOnCurve( Edge & BhAB, BamgVertex & vA, BamgVertex & vB,/{{{/
3178	double theta,BamgVertex & R,VertexOnEdge & BR,VertexOnGeom & GR) {
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;
3184	BamgVertex * pvA=&vA, * pvB=&vB;
3185	if (vA.IndexInTriangle == IsVertexOnVertex){
3186	pA=vA.BackgroundVertexHook;
3187	}
3188	else if (vA.IndexInTriangle == IsVertexOnEdge){
3189	pA=vA.BackgroundEdgeHook->be;
3190	tA=vA.BackgroundEdgeHook->abcisse;
3191	}
3192	else {
3193	_error_("ProjectOnCurve On BamgVertex " << BTh.GetId(vA) << " forget call to SetVertexFieldOnBTh");
3194	}
3195
3196	if (vB.IndexInTriangle == IsVertexOnVertex){
3197	pB=vB.BackgroundVertexHook;
3198	}
3199	else if(vB.IndexInTriangle == IsVertexOnEdge){
3200	pB=vB.BackgroundEdgeHook->be;
3201	tB=vB.BackgroundEdgeHook->abcisse;
3202	}
3203	else {
3204	_error_("ProjectOnCurve On BamgVertex " << BTh.GetId(vB) << " forget call to SetVertexFieldOnBTh");
3205	}
3206	Edge * e = &BhAB;
3207	if (!pA \|\| !pB \|\| !e){
3208	_error_("!pA \|\| !pB \|\| !e");
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){
3214	_error_("e<BTh.edges \|\| e>=BTh.edges+BTh.nbe");
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
3221	if( vA.IndexInTriangle == IsVertexOnEdge)
3222	{ // find the start edge
3223	e = vA.BackgroundEdgeHook->be;
3224
3225	}
3226	else if (vB.IndexInTriangle == IsVertexOnEdge)
3227	{
3228	theta = 1-theta;
3229	Exchange(tA,tB);
3230	Exchange(pA,pB);
3231	Exchange(pvA,pvB);
3232	Exchange(A,B);
3233	e = vB.BackgroundEdgeHook->be;
3234
3235	}
3236	else{ // do the search by walking
3237	_error_("case not supported yet");
3238	}
3239
3240	// find the direction of walking with direction of edge and pA,PB;
3241	R2 AB=B-A;
3242
3243	double cosE01AB = (( (R2) (e)[1] - (R2) (e)[0] ) , AB);
3244	int kkk=0;
3245	int direction = (cosE01AB>0) ? 1 : 0;
3246
3247	// double l=0; // length of the edge AB
3248	double abscisse = -1;
3249
3250	for (int step=0;step<2;step++){
3251	// 2 times algo:
3252	// 1 for computing the length l
3253	// 2 for find the vertex
3254	int iii;
3255	BamgVertex v0=pvA,v1;
3256	Edge neee,eee;
3257	double lg =0; // length of the curve
3258	double te0;
3259	// we suppose take the curve's abcisse
3260	for ( eee=e,iii=direction,te0=tA;
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;
3265	_assert_(kkk<100);
3266	_assert_(eee);
3267	double lg0 = lg;
3268	double dp = LengthInterpole(v0->m,v1->m,(R2) v1 - (R2) v0);
3269	lg += dp;
3270	if (step && abscisse <= lg) { // ok we find the geom edge
3271	double sss = (abscisse-lg0)/dp;
3272	double thetab = te0(1-sss)+ sssiii;
3273	_assert_(thetab>=0 && thetab<=1);
3274	BR = VertexOnEdge(&R,eee,thetab);
3275	return Gh.ProjectOnCurve(*eee,thetab,R,GR);
3276	}
3277	}
3278	// we find the end
3279	if (v1 != pvB){
3280	if (( void*) v1 == pB)
3281	tB = iii;
3282
3283	double lg0 = lg;
3284	_assert_(eee);
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;
3293	_assert_(thetab>=0 && thetab<=1);
3294	BR = VertexOnEdge(&R,eee,thetab);
3295	return Gh.ProjectOnCurve(*eee,thetab,R,GR);
3296	}
3297	}
3298	abscisse = lg*theta;
3299
3300	}
3301	_error_("Big bug...");
3302	return 0; // just for the compiler
3303	}
3304	/}}}/
3305	void Mesh::ReconstructExistingMesh(){/{{{/
3306	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FillHoleInMesh)/
3307
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	*/
3315
3316	/Intermediary/
3317	int verbose=0;
3318
3319	// generation of the integer coordinate
3320
3321	// find extrema coordinates of vertices pmin,pmax
3322	long i;
3323	if(verbose>2) _printf_(" Reconstruct mesh of " << nbv << " vertices\n");
3324
3325	//initialize orderedvertices
3326	_assert_(orderedvertices);
3327	for (i=0;i<nbv;i++) orderedvertices[i]=0;
3328
3329	//Initialize nbsubdomains
3330	nbsubdomains =0;
3331
3332	/* generation of triangles adjacency*/
3333
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	}
3343
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++){
3349
3350	//Add current triangle edge to edge4
3351	long k =edge4->SortAndAdd(GetId(triangles[i][VerticesOfTriangularEdge[j][0]]),GetId(triangles[i][VerticesOfTriangularEdge[j][1]]));
3352
3353	long invisible=triangles[i].Hidden(j);
3354
3355	//If the edge has not been added to st, add it
3356	if(st[k]==-1) st[k]=3*i+j;
3357
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)));
3362
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);
3366
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	}
3370
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	}
3375	}
3376	}
3377
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	}
3387
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;
3397	}
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	}
3407	}
3408	}
3409	}
3410	if(k) {
3411	_error_(k << " boundary edges (from the geometry) are not defined as mesh edges");
3412	}
3413
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	}
3421
3422	Triangle* savetriangles=triangles;
3423	long savenbt=nbt;
3424	long savemaxnbt=maxnbt;
3425	SubDomain* savesubdomains=subdomains;
3426	subdomains=0;
3427
3428	long Nbtriafillhole=2*nbvb;
3429	Triangle* triafillhole=new Triangle[Nbtriafillhole];
3430	triangles = triafillhole;
3431
3432	nbt=2;
3433	maxnbt= Nbtriafillhole;
3434
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]);
3442
3443	/Reconstruct mesh beginning with 2 triangles/
3444	BamgVertex * v0=orderedvertices[0], *v1=orderedvertices[1];
3445
3446	triangles[0](0) = NULL; // Infinite vertex
3447	triangles[0](1) = v0;
3448	triangles[0](2) = v1;
3449
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);
3459
3460	triangles[0].det = -1; // boundary triangles
3461	triangles[1].det = -1; // boundary triangles
3462
3463	triangles[0].SetSingleVertexToTriangleConnectivity();
3464	triangles[1].SetSingleVertexToTriangleConnectivity();
3465
3466	triangles[0].link=&triangles[1];
3467	triangles[1].link=&triangles[0];
3468
3469	if (!quadtree) delete quadtree; //ReInitialise;
3470	quadtree = new BamgQuadtree(this,0);
3471	quadtree->Add(*v0);
3472	quadtree->Add(*v1);
3473
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	}
3484
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	}
3495	}
3496	}
3497	if(nbloss) {
3498	_error_("we lost " << nbloss << " existing edges other " << knbe);
3499	}
3500
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);
3518
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	}
3524	}
3525	}
3526	}
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	}
3536	}
3537	_assert_(savenbt+NbTfillHoll<=savemaxnbt);
3538
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	}
3546	}
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	}
3562	}
3563	}
3564
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();
3576
3577	delete edge4;
3578	delete [] st;
3579	for (i=0;i<nbv;i++) quadtree->Add(vertices[i]);
3580
3581	SetVertexFieldOn();
3582
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");
3601
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	}
3605	}
3606	}
3607	}
3608	}
3609	}
3610	/}}}/
3611	void Mesh::TrianglesRenumberBySubDomain(bool justcompress){/{{{/
3612	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingTheTriangleBySubDomain)/
3613
3614	long *renu= new long[nbt];
3615	Triangle t0,t,*te=triangles+nbt;
3616	long k=0,it,i,j;
3617
3618	for ( it=0;it<nbt;it++)
3619	renu[it]=-1; // outside triangle
3620	for ( i=0;i<nbsubdomains;i++)
3621	{
3622	t=t0=subdomains[i].head;
3623	if (!t0){ // not empty sub domain
3624	_error_("!t0");
3625	}
3626	do {
3627	long kt = GetId(t);
3628	if (kt<0 \|\| kt >= nbt ){
3629	_error_("kt<0 \|\| kt >= nbt");
3630	}
3631	if (renu[kt]!=-1){
3632	_error_("renu[kt]!=-1");
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){
3649	_error_("k != nbt");
3650	}
3651	// do the change on all the pointeur
3652	for ( it=0;it<nbt;it++)
3653	triangles[it].Renumbering(triangles,te,renu);
3654
3655	for ( i=0;i<nbsubdomains;i++)
3656	subdomains[i].head=triangles+renu[GetId(subdomains[i].head)];
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	}
3677	/}}}/
3678	void Mesh::SetIntCoor(const char * strfrom) {/{{{/
3679	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SetIntCoor)/
3680
3681	/Set integer coordinate for existing vertices/
3682
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;
3696
3697	//Compute coefIcoor
3698	int MaxICoord = 1073741823; //2^30 - 1 = =111...111 (29 times one)
3699	coefIcoor= (MaxICoord)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
3700	if (coefIcoor<=0){
3701	_error_("coefIcoor should be positive, a problem in the geometry is likely");
3702	}
3703
3704	// generation of integer coord
3705	for (i=0;i<nbv;i++) {
3706	vertices[i].i = R2ToI2(vertices[i].r);
3707	}
3708
3709	// computation of the det
3710	int number_of_errors=0;
3711	for (i=0;i<nbt;i++) {
3712	BamgVertex* v0 = triangles[i](0);
3713	BamgVertex* v1 = triangles[i](1);
3714	BamgVertex* v2 = triangles[i](2);
3715
3716	//If this is not a boundary triangle
3717	if (v0 && v1 && v2){
3718
3719	/Compute determinant/
3720	triangles[i].det= det(v0->GetIntegerCoordinates(),v1->GetIntegerCoordinates(),v2->GetIntegerCoordinates());
3721
3722	/Check that determinant is positive/
3723	if (triangles[i].det <=0){
3724
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	}
3730	}
3731	}
3732
3733	//else, set as -1
3734	else triangles[i].det=-1;
3735	}
3736
3737	if (number_of_errors) _error_("Fatal error: some triangles have negative areas, see above");
3738	}
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)/
3742
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)
3753	for ( i=0;i<nbv;i++)
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));
3770	for ( i=0;i<nbv;i++)
3771	if (tstart[i] != &vide) // not a boundary vertex
3772	NbSwap += vertices[i].Optim(1);
3773	if (verbose>3) _printf_(" move max = " << pow(delta,0.5) << ", iteration = " << k << ", nb of swap = " << NbSwap << "\n");
3774	}
3775
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)/
3782
3783	long int verbose=0;
3784
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);
3793
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");
3807	}
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)");
3815	}
3816	BamgVertex *pvj = (ta.EdgeVertex(0));
3817	BamgVertex & vj = *pvj;
3818	if(pvj){
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) {
3826	double hi = ll/vi.m.Length(Aij.x,Aij.y);
3827	double hj = ll/vj.m.Length(Aij.x,Aij.y);
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
3839	{
3840	double li = vi.m.Length(Aij.x,Aij.y);
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;
3844	}
3845	}
3846	if ( &vj == pvj0 ) break;
3847	}
3848	}
3849	Head0 = Head1;
3850	Head1 = -1;
3851	Exchange(first_np_or_next_t0,first_np_or_next_t1);
3852	}
3853	if(verbose>2) _printf_(" number of iterations = " << kch << "\n");
3854	delete [] first_np_or_next_t0;
3855	delete [] first_np_or_next_t1;
3856	}
3857	/}}}/
3858	long Mesh::SplitInternalEdgeWithBorderVertices(){/{{{/
3859	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SplitInternalEdgeWithBorderVertices)/
3860
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)){
3871
3872	Triangle *ptt = t.TriangleAdj(j);
3873	Triangle &tt = *ptt;
3874
3875	if( ptt && tt.link && it < GetId(tt))
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++;
3887	}
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);
3900
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	}
3907
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++;
3915	}
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	}
3921	}
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;
3926	}
3927	/}}}/
3928	I2 Mesh::R2ToI2(const R2 & P) const {/{{{/
3929	return I2( (int) (coefIcoor(P.x-pmin.x)),(int) (coefIcoor(P.y-pmin.y)) );
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)/
3938
3939	Triangle * t=0;
3940	int j,jp,jn,jj;
3941	int counter;
3942
3943	/Get starting triangle. Take tsart if provided/
3944	if (tstart) t=tstart;
3945
3946	/Else find the closest Triangle using the quadtree/
3947	else {
3948
3949	/Check that the quadtree does exist/
3950	if (!quadtree) _error_("no starting triangle provided and no quadtree available");
3951
3952	/Call NearestVertex/
3953	BamgVertex *a = quadtree->NearestVertex(B.x,B.y) ;
3954
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);
3959
3960	/Get starting triangle/
3961	t = a->t;
3962	_assert_(t>=triangles && t<triangles+nbt);
3963	}
3964
3965	Icoor2 detop ;
3966
3967	/initialize number of test triangle/
3968	counter=0;
3969
3970	/The initial triangle might be outside/
3971	while (t->det < 0){
3972
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	}
3985
3986	jj=0;
3987	detop = det((t)(VerticesOfTriangularEdge[jj][0]),(t)(VerticesOfTriangularEdge[jj][1]),B);
3988
3989	while(t->det>0){
3990
3991	/Increase counter/
3992	if (++counter>=10000) _error_("Maximum number of iteration reached (threshold = " << counter << ").");
3993
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];
4000
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
4009
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	}
4021
4022	if (t->det<0) // outside triangle
4023	det3[0]=det3[1]=det3[2]=-1,det3[OppositeVertex[jj]]=detop;
4024	return t;
4025	}
4026	/}}}/
4027	void Mesh::TriangleIntNumbering(long* renumbering){/{{{/
4028
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;
4035	}
4036	/}}}/
4037	long Mesh::TriangleReferenceList(long* reft) const {/{{{/
4038	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ConsRefTriangle)/
4039
4040	Triangle t0,t;
4041	long k=0, num;
4042
4043	//initialize all triangles as -1 (outside)
4044	for (int it=0;it<nbt;it++) reft[it]=-1;
4045
4046	//loop over all subdomains
4047	for (int i=0;i<nbsubdomains;i++){
4048
4049	//first triangle of the subdomain i
4050	t=t0=subdomains[i].head;
4051
4052	//check that the subdomain is not empty
4053	if (!t0){ _error_("At least one subdomain is empty");}
4054
4055	//loop
4056	do{
4057	k++;
4058
4059	//get current triangle number
4060	num = GetId(t);
4061
4062	//check that num is in [0 nbt[
4063	_assert_(num>=0 && num<nbt);
4064
4065	//reft of this triangle is the subdomain number
4066	reft[num]=i;
4067
4068	} while (t0 != (t=t->link));
4069	//stop when all triangles of subdomains have been tagged
4070
4071	}
4072	return k;
4073	}
4074	/}}}/
4075	void Mesh::Triangulate(double* x,double* y,int nods){/{{{/
4076
4077	int verbose=0;
4078	int i;
4079	Metric M1(1);
4080
4081	/Initialize mesh/
4082	Init(nods);//this resets nbv to 0
4083	nbv=nods;
4084
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);
4098	}
4099	/}}}/
4100	void Mesh::TriangulateFromGeom0(BamgOpts* bamgopts){/{{{/
4101	/Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/GeomToTriangles0)/
4102	/Generate mesh from geometry/
4103
4104	/Intermediaries/
4105	int i,k;
4106	int nbcurves = 0;
4107	int NbNewPoints,NbEdgeCurve;
4108	double lcurve,lstep,s;
4109	const int MaxSubEdge = 10;
4110
4111	R2 AB;
4112	GeomVertex a, b;
4113	BamgVertex va,vb;
4114	GeomEdge *e;
4115
4116	// add a ref to GH to make sure that it is not destroyed by mistake
4117	Gh.NbRef++;
4118
4119	/Get options/
4120	int verbose=bamgopts->verbose;
4121
4122	//build background mesh flag (1 if background, else 0)
4123	bool background=(&BTh != this);
4124
4125	/Build VerticesOnGeomVertex/
4126
4127	//Compute the number of geometrical vertices that we are going to use to mesh
4128	for (i=0;i<Gh.nbv;i++){
4129	if (Gh[i].Required()) NbVerticesOnGeomVertex++;
4130	}
4131	//allocate
4132	VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];
4133	if(NbVerticesOnGeomVertex >= maxnbv) _error_("too many vertices on geometry: " << NbVerticesOnGeomVertex << " >= " << maxnbv);
4134	_assert_(nbv==0);
4135	//Build VerticesOnGeomVertex
4136	for (i=0;i<Gh.nbv;i++){
4137	/* Add vertex only if required*/
4138	if (Gh[i].Required()) {//Gh vertices Required
4139
4140	//Add the vertex
4141	_assert_(nbv<maxnbv);
4142	vertices[nbv]=Gh[i];
4143
4144	//Add pointer from geometry (Gh) to vertex from mesh (Th)
4145	Gh[i].MeshVertexHook=vertices+nbv;
4146
4147	//Build VerticesOnGeomVertex for current point
4148	VerticesOnGeomVertex[nbv]=VertexOnGeom(*Gh[i].MeshVertexHook,Gh[i]);
4149
4150	//nbv increment
4151	nbv++;
4152	}
4153	}
4154
4155	/Build VerticesOnGeomEdge/
4156
4157	//check that edges is still empty (Init)
4158	_assert_(!edges);
4159
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
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 */
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();
4173	nbcurves=0;
4174
4175	//go through the edges of the geometry
4176	for (i=0;i<Gh.nbe;i++){
4177
4178	//ei = current Geometrical edge
4179	GeomEdge &ei=Gh.edges[i];
4180
4181	//loop over the two vertices of the edge ei
4182	for(int j=0;j<2;j++) {
4183
4184	/Take only required vertices (corner->beginning of a new curve)/
4185	if (!ei.Mark() && ei[j].Required()){
4186
4187	long nbvend=0;
4188	Edge* PreviousNewEdge=NULL;
4189	lstep = -1;
4190
4191	/If Edge is required (do that only once for the 2 vertices)/
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);
4200
4201	//check that edges has been allocated
4202	_assert_(edges);
4203	edges[nbe].v[0]=a->MeshVertexHook;
4204	edges[nbe].v[1]=b->MeshVertexHook;;
4205	edges[nbe].ReferenceNumber = e->ReferenceNumber;
4206	edges[nbe].GeomEdgeHook = e;
4207	edges[nbe].adj[0] = 0;
4208	edges[nbe].adj[1] = 0;
4209	nbe++;
4210	}
4211	}
4212	}
4213
4214	/If Edge is not required: we are on a curve/
4215	else {
4216	for (int kstep=0;kstep<=step;kstep++){
4217	//kstep=0, compute number of edges (discretize curve)
4218	//kstep=1 create the points and edge
4219	PreviousNewEdge=0;
4220	NbNewPoints=0;
4221	NbEdgeCurve=0;
4222	if (nbvend>=maxnbv) _error_("maximum number of vertices too low! Check the domain outline or increase maxnbv");
4223	lcurve =0;
4224	s = lstep; //-1 initially, then length of each sub edge
4225
4226	/reminder: i = edge number, j=[0;1] vertex index in edge/
4227	k=j; // k = vertex index in edge (0 or 1)
4228	e=&ei; // e = reference of current edge
4229	a=ei(k); // a = pointer toward the kth vertex of the current edge
4230	va = a->MeshVertexHook; // va = pointer toward mesh vertex associated
4231	e->SetMark(); // Mark edge
4232
4233	/Loop until we reach the end of the curve/
4234	for(;;){
4235	k = 1-k; // other vertx index of the curve
4236	b = (*e)(k); // b = pointer toward the other vertex of the current edge
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
4239	Metric MB = background ? BTh.MetricAt(b->r) :b->m ; //Get metric associated to B
4240	double ledge = (MA.Length(AB.x,AB.y) + MB.Length(AB.x,AB.y))/2; //Get edge length in metric
4241
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). */
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
4253	if (ledge < 1.5){
4254	//if ledge < 1.5 (between one and 2), take the edge as is
4255	lSubEdge[0] = ledge;
4256	}
4257	//else, divide the edge
4258	else {
4259	//compute number of subedges (division of the edge), Maximum is 10
4260	NbSubEdge = Min( MaxSubEdge, (int) (ledge +0.5));
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)*/
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);
4274	MBs= background ? BTh.MetricAt(B) : Metric(1-x,MA,x,MB);
4275	AB = A-B;
4276	lSubEdge[kk]=(ledge+=(MAs.Length(AB.x,AB.y)+MBs.Length(AB.x,AB.y))/2);
4277	}
4278	}
4279
4280	double lcurveb = lcurve+ledge;
4281
4282	/Now, create corresponding points/
4283	while(s>=lcurve && s<=lcurveb && nbv<nbvend){
4284
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
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	}
4305	_assert_(kk1!=kk0);
4306
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;
4316
4317	// new vertex on edge
4318	vb = &vertices[nbv++];
4319	vb->m = Metric(aa,a->m,bb,b->m);
4320	vb->ReferenceNumber = e->ReferenceNumber;
4321	double abcisse = k ? bb : aa;
4322	vb->r = e->F(abcisse);
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;
4329	edges[nbe].ReferenceNumber =e->ReferenceNumber;
4330	edges[nbe].GeomEdgeHook = e;
4331	edges[nbe].adj[0] = PreviousNewEdge;
4332	if(PreviousNewEdge) PreviousNewEdge->adj[1]=&edges[nbe];
4333	PreviousNewEdge=edges+nbe;
4334	nbe++;
4335	va = vb;
4336	}
4337
4338	/We just added one edge to the curve: Go to the next one/
4339	lcurve = lcurveb;
4340	e->SetMark();
4341	a=b;
4342
4343	/If b is required, we are on a new curve->break/
4344	if (b->Required()) break;
4345	int kprev=k;
4346	k = e->AdjVertexIndex[kprev];// next vertices
4347	e = e->Adj[kprev];
4348	_assert_(e);
4349	}// for(;;)
4350	vb = b->MeshVertexHook;
4351
4352	/Number of edges in the last disretized curve/
4353	NbEdgeCurve = Max((long) (lcurve +0.5), (long) 1);
4354	/Number of internal vertices in the last disretized curve/
4355	NbNewPoints = NbEdgeCurve-1;
4356	if(!kstep){
4357	NbVerticesOnGeomEdge0 += NbNewPoints;
4358	nbcurves++;
4359	}
4360	nbvend=nbv+NbNewPoints;
4361	lstep = lcurve / NbEdgeCurve; //approximately one
4362	}// end of curve --
4363	if (edges) { // last edges of the curves
4364	edges[nbe].v[0]=va;
4365	edges[nbe].v[1]=vb;
4366	edges[nbe].ReferenceNumber = e->ReferenceNumber;
4367	edges[nbe].GeomEdgeHook = e;
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) {
4379	_assert_(!edges);
4380	_assert_(!VerticesOnGeomEdge);
4381
4382	edges = new Edge[nbex=nbe];
4383	if(NbVerticesOnGeomEdge0) VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge0];
4384
4385	// do the vertex on a geometrical vertex
4386	_assert_(VerticesOnGeomEdge \|\| NbVerticesOnGeomEdge0==0);
4387	NbVerticesOnGeomEdge0 = NbVerticesOnGeomEdge;
4388	}
4389	else{
4390	_assert_(NbVerticesOnGeomEdge==NbVerticesOnGeomEdge0);
4391	}
4392	}
4393
4394	//Insert points inside existing triangles
4395	if (verbose>4) _printf_(" -- current number of vertices = " << nbv << "\n");
4396	if (verbose>3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n");
4397	if (verbose>3) _printf_(" Inserting boundary points\n");
4398	Insert(bamgopts->random);
4399
4400	//Force the boundary
4401	if (verbose>3) _printf_(" Forcing boundaries\n");
4402	ForceBoundary();
4403
4404	//Extract SubDomains
4405	if (verbose>3) _printf_(" Extracting subdomains\n");
4406	FindSubDomain();
4407
4408	if (verbose>3) _printf_(" Inserting internal points\n");
4409	NewPoints(*this,bamgopts,0) ;
4410	if (verbose>4) _printf_(" -- current number of vertices = " << nbv << "\n");
4411	}
4412	/}}}/
4413	void Mesh::TriangulateFromGeom1(BamgOpts* bamgopts,int KeepVertices){ /{{{/
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
4419	Gh.NbRef++;// add a ref to Gh
4420
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	*************************************************************************/
4437
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
4441
4442	//Initialize new mesh
4443	BTh.SetVertexFieldOn();
4444	int* bcurve = new int[Gh.nbcurves]; //
4445
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 */
4452
4453	NbVerticesOnGeomVertex=0;
4454	NbVerticesOnGeomEdge=0;
4455
4456	/STEP 1 copy of Required vertices/
4457
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	}
4464
4465	VerticesOnGeomVertex = new VertexOnGeom[ NbVerticesOnGeomVertex];
4466	VertexOnBThVertex = new VertexOnVertex[NbVerticesOnGeomVertex];
4467
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]);
4476	nbv++;
4477	}
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	}
4489	}
4490	_assert_(NbVertexOnBThVertex==NbVerticesOnGeomVertex); /This might be due to MaxCornerAngle too small/
4491
4492	/STEP 2: reseed boundary edges/
4493
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;
4498
4499	/Loop over the backgrounf mesh BTh edges/
4500	for (int iedge=0;iedge<BTh.nbe;iedge++){
4501	Edge &ei = BTh.edges[iedge];
4502
4503	/Loop over the 2 vertices of the current edge/
4504	for(int je=0;je<2;je++){
4505
4506	/* If one of the vertex is required we are in a new curve*/
4507	if (ei[je].GeomEdgeHook->IsRequiredVertex()){
4508
4509	/Get curve number/
4510	int nc=ei.GeomEdgeHook->CurveNumber;
4511
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	}
4535
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
4540
4541	long nbex=0,NbVerticesOnGeomEdgex=0;
4542	for (int step=0; step <2;step++){
4543
4544	long NbOfNewPoints=0;
4545	long NbOfNewEdge=0;
4546	long iedge;
4547	Gh.UnMarkEdges();
4548	double L=0;
4549
4550	/Go through all geometrical curve/
4551	for (int icurve=0;icurve<Gh.nbcurves;icurve++){
4552
4553	/Get edge and vertex (index) of background mesh on this curve/
4554	iedge=bcurve[icurve]/2;
4555	int jedge=bcurve[icurve]%2;
4556
4557	/Get edge of Bth with index iedge/
4558	Edge &ei = BTh.edges[iedge];
4559
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)
4563
4564	/Do phase 0 to step/
4565	for(int phase=0;phase<=step;phase++){
4566
4567	/Current curve pointer/
4568	Curve *curve= Gh.curves+icurve;
4569
4570	/Get index of current curve/
4571	int icurveequi= Gh.GetId(curve);
4572
4573	/For phase 0, check that we are at the begining of the curve only/
4574	if(phase==0 && icurveequi!=icurve) continue;
4575
4576	int k0=jedge,k1;
4577	Edge* pe= BTh.edges+iedge;
4578	int iedgeequi=bcurve[icurveequi]/2;
4579	int jedgeequi=bcurve[icurveequi]%2;
4580
4581	int k0equi=jedgeequi,k1equi;
4582	Edge * peequi= BTh.edges+iedgeequi;
4583	GeomEdge *ongequi = peequi->GeomEdgeHook;
4584
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;
4594
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;
4614
4615	if (phase){
4616	// computation of the new points for the given curve
4617	while ((i!=NbCreatePointOnCurve) && sNew<=L) {
4618
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);
4625
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]");
4633	}
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;
4648
4649	if (PreviousNewEdge) PreviousNewEdge->adj[1]=e;
4650	PreviousNewEdge=e;
4651	A0=A1;
4652	sNew += Lstep;
4653	if (++i== NbCreatePointOnCurve) break;
4654	}
4655	}
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
4674	}
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	}
4697	}
4698	}// end of curve loop
4699
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;
4715	}
4716	}
4717	_assert_(nbe!=0);
4718	delete [] bcurve;
4719
4720	//Insert points inside existing triangles
4721	if (verbose>4) _printf_(" -- current number of vertices = " << nbv << "\n");
4722	if (verbose>3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n");
4723	if (verbose>3) _printf_(" Inserting boundary points\n");
4724	Insert(bamgopts->random);
4725
4726	//Force the boundary
4727	if (verbose>3) _printf_(" Forcing boundaries\n");
4728	ForceBoundary();
4729
4730	//Extract SubDomains
4731	if (verbose>3) _printf_(" Extracting subdomains\n");
4732	FindSubDomain();
4733
4734	if (verbose>3) _printf_(" Inserting internal points\n");
4735	NewPoints(BTh,bamgopts,KeepVertices) ;
4736	if (verbose>4) _printf_(" -- current number of vertices = " << nbv << "\n");
4737	}
4738	/}}}/
4739
4740	/Intermediary/
4741	AdjacentTriangle CloseBoundaryEdge(I2 A,Triangle t, double &a,double &b) {/{{{*/
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){
4747	_error_("k<0");
4748	}
4749	int kkk=0;
4750	Icoor2 IJ_IA,IJ_AJ;
4751	AdjacentTriangle edge(t,OppositeEdge[k]);
4752	for (;;edge = dir >0 ? Next(Adj(Next(edge))) : Previous(Adj(Previous(edge)))) {
4753	kkk++;
4754	if (kkk>=1000){
4755	_error_("kkk>=1000");
4756	}
4757	BamgVertex &vI = *edge.EdgeVertex(0);
4758	BamgVertex &vJ = *edge.EdgeVertex(1);
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){
4772	_error_("IJ2==0");
4773	}
4774	a= IJ_AJ/IJ2;
4775	b= IJ_IA/IJ2;
4776	return edge;
4777	}
4778	}
4779	/}}}/
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)/
4782
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
4789
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
4793
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);
4805	}
4806
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;
4813
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	}
4839	}
4840	tta = tc;
4841	k++;
4842	if (k>=2000){
4843	_error_("k>=2000");
4844	}
4845	if ( vbegin == v2 ) return -1;// error
4846	}
4847
4848	tta.SetLock();
4849	taret=tta;
4850	a.Optim(1,0);
4851	b.Optim(1,0);
4852	return NbSwap;
4853	}
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);
4891
4892	t1->det = det1;
4893	t2->det = det2;
4894
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) {/{{{/
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
4909	AdjacentTriangle tt2 = Adj(tt1);
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 ){
4913	_error_("a1<0 \|\| a1>=3");
4914	}
4915
4916	BamgVertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]];
4917	BamgVertex & s1= (*t1)[OppositeVertex[a1]];
4918	BamgVertex & s2= (*t2)[OppositeVertex[a2]];
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)){
4924	_error_("(det1<=0 ) \|\| (det2<=0)");
4925	}
4926	if ( (detsa>=0) \|\| (detsb<=0) ){ // [a,b] cut infinite line va,bb
4927	_error_("(detsa>=0) \|\| (detsb<=0)");
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	}
4956	else { // changement de direction
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
4967	if(ToSwap){
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;}
4980	}
4981
4982	}
4983	return ret;
4984	}
4985	/}}}/
4986	}

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

Download in other formats: