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

Last change on this file was 22085, checked in by Eric.Larour, 8 years ago

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