Context Navigation

Mesh.cpp@ 22085

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

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

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

Context Navigation

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

Download in other formats: