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