Context Navigation

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

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

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

Download in other formats: