source: issm/trunk/src/c/classes/bamg/Geometry.cpp@ 13395

Last change on this file since 13395 was 13395, checked in by Mathieu Morlighem, 12 years ago

merged trunk-jpl and trunk for revision 13393

File size: 30.2 KB
RevLine 
[3913]1#include <cstdio>
2#include <cstring>
3#include <cmath>
4#include <ctime>
5
6#include "../objects.h"
7#include "../../include/include.h"
8#include "../../shared/Exceptions/exceptions.h"
9
10namespace bamg {
11
12 static const Direction NoDirOfSearch=Direction();
13
14 /*Constructors/Destructors*/
[12365]15 /*FUNCTION Geometry::Geometry(){{{*/
[5120]16 Geometry::Geometry(){
17 Init();
18 }
[5130]19 /*}}}*/
[12365]20 /*FUNCTION Geometry::Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts){{{*/
[5120]21 Geometry::Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts){
22 Init();
23 ReadGeometry(bamggeom,bamgopts);
[5397]24 PostRead();
[5120]25 }
[5130]26 /*}}}*/
[12365]27 /*FUNCTION Geometry::Geometry(const Geometry & Gh) (COPY operator){{{*/
[3913]28 Geometry::Geometry(const Geometry & Gh) {
29 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Geometry)*/
30
31 long i;
32 *this = Gh;
33 NbRef =0;
34 quadtree=0;
[5573]35 vertices = nbv ? new GeomVertex[nbv] : NULL;
36 edges = nbe ? new GeomEdge[nbe]:NULL;
[5148]37 curves= nbcurves ? new Curve[nbcurves]:NULL;
[5573]38 subdomains = nbsubdomains ? new GeomSubDomain[nbsubdomains]:NULL;
[3913]39 for (i=0;i<nbe;i++)
40 edges[i].Set(Gh.edges[i],Gh,*this);
[5148]41 for (i=0;i<nbcurves;i++)
[3913]42 curves[i].Set(Gh.curves[i],Gh,*this);
[5148]43 for (i=0;i<nbsubdomains;i++)
[3913]44 subdomains[i].Set(Gh.subdomains[i],Gh,*this);
45 }
[12365]46 /*}}}*/
47 /*FUNCTION Geometry::~Geometry(){{{*/
[3913]48 Geometry::~Geometry() {
49 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/~Geometry)*/
[12509]50 if(NbRef>0){ _printString_("Trying to delete geometry and NbRef>0, probably due to an error"); return;}
[5208]51 if(vertices) delete [] vertices; vertices=0;
52 if(edges) delete [] edges; edges=0;
53 if(quadtree) delete quadtree; quadtree=0;
54 if(curves) delete [] curves; curves=0;nbcurves=0;
[3913]55 if(subdomains) delete [] subdomains;subdomains=0;
[5120]56 Init();
[3913]57 }
[12365]58 /*}}}*/
[3913]59
60 /*IO*/
[12365]61 /*FUNCTION Geometry::ReadGeometry{{{*/
[3913]62 void Geometry::ReadGeometry(BamgGeom* bamggeom,BamgOpts* bamgopts){
63
64 int verbose;
[5120]65 nbv=0;
66 nbe=0;
[5148]67 nbcurves=0;
[3913]68
69 double Hmin = HUGE_VAL;// the infinie value
[5143]70 int i,j,k,n,i0,i1,i2,i3;
[3913]71
[5120]72 /*initialize some variables*/
[5143]73 verbose= bamgopts->verbose;
74 nbv = bamggeom->VerticesSize[0];
75 nbe = bamggeom->EdgesSize[0];
[3913]76
77 //some checks
[13036]78 if (bamggeom->Vertices==NULL) _error_("the domain provided does not contain any vertex");
79 if (bamggeom->Edges==NULL) _error_("the domain provided does not contain any edge");
[3913]80
81 //Vertices
82 if (bamggeom->Vertices){
[12509]83 if(verbose>5) _printLine_(" processing Vertices");
[13036]84 if (bamggeom->VerticesSize[1]!=3) _error_("Vertices should have 3 columns");
[5573]85 vertices = new GeomVertex[nbv];
[3913]86 for (i=0;i<nbv;i++) {
87 vertices[i].r.x=(double)bamggeom->Vertices[i*3+0];
88 vertices[i].r.y=(double)bamggeom->Vertices[i*3+1];
89 vertices[i].ReferenceNumber=(long)bamggeom->Vertices[i*3+2];
90 vertices[i].DirOfSearch=NoDirOfSearch;
91 vertices[i].color =0;
[5124]92 vertices[i].type=0;
[3913]93 }
[5120]94 /*find domain extrema (pmin,pmax) that will define the square box used for by the quadtree*/
[3913]95 pmin = vertices[0].r;
96 pmax = vertices[0].r;
97 for (i=0;i<nbv;i++) {
98 pmin.x = Min(pmin.x,vertices[i].r.x);
99 pmin.y = Min(pmin.y,vertices[i].r.y);
100 pmax.x = Max(pmax.x,vertices[i].r.x);
101 pmax.y = Max(pmax.y,vertices[i].r.y);
102 }
[5120]103 /*Offset pmin and pmax to avoid round-off errors*/
104 R2 offset = (pmax-pmin)*0.05;
105 pmin -= offset;
106 pmax += offset;
107 /*coefIcoor is the coefficient used for integer coordinates:
108 * (x-pmin.x)
109 * Icoor x = (2^30 -1) ------------
110 * D
111 * where D is the longest side of the domain (direction x or y)
112 * so that (x-pmin.x)/D is in ]0 1[
[5177]113 *
114 * coefIcoor = (2^30 -1)/D
[5120]115 */
116 coefIcoor=(MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
[13036]117 if(coefIcoor<=0) _error_("coefIcoor should be positive");
[3913]118 }
119 else{
[13036]120 _error_("No BamgVertex provided");
[3913]121 }
122
123 //Edges
124 if (bamggeom->Edges){
[5143]125 R2 zerovector(0,0);
126 double* verticeslength=NULL;
[3913]127
[12509]128 if(verbose>5) _printLine_(" processing Edges");
[13036]129 if (bamggeom->EdgesSize[1]!=3) _error_("Edges should have 3 columns");
[5573]130 edges = new GeomEdge[nbe];
[3913]131
[5143]132 //initialize verticeslength (sum of the lengths of the edges holding vertex)
133 verticeslength = new double[nbv];
134 for(i=0;i<nbv;i++) verticeslength[i]=0;
[3913]135
[5143]136 /*Loop over the edges*/
[3913]137 for (i=0;i<nbe;i++){
138
139 i1=(int)bamggeom->Edges[i*3+0]-1; //-1 for C indexing
140 i2=(int)bamggeom->Edges[i*3+1]-1; //-1 for C indexing
[5143]141 edges[i].v[0]= vertices + i1; //pointer toward vertex i1 (=&vertices[i1])
142 edges[i].v[1]= vertices + i2; //pointer toward vertex i2
[5148]143 edges[i].ReferenceNumber=(long)bamggeom->Edges[i*3+2];
[3913]144
145 //get length of edge
[5143]146 R2 x12=vertices[i2].r-vertices[i1].r;
[3913]147 double l12=sqrt((x12,x12));
148 Hmin=Min(Hmin,l12);
149
150 //initialize other fields
[5143]151 edges[i].tg[0]=zerovector;
152 edges[i].tg[1]=zerovector;
[5401]153 edges[i].AdjVertexIndex[0] = edges[i].AdjVertexIndex[1] = -1;
[5143]154 edges[i].Adj[0] = edges[i].Adj[1] = NULL;
[5146]155 edges[i].type = 0;
[3913]156
157 //Cracked?
[5148]158 if (edges[i].ReferenceNumber!=1) edges[i].SetCracked();
[3913]159
160 //prepare metric
161 vertices[i1].color++;
162 vertices[i2].color++;
[5143]163 verticeslength[i1] += l12;
164 verticeslength[i2] += l12;
[3913]165 }
166
[5397]167 // definition the default of the given mesh size
[5143]168 for (i=0;i<nbv;i++) {
169 if (vertices[i].color > 0)
170 vertices[i].m=Metric(verticeslength[i] /(double) vertices[i].color);
171 else
172 vertices[i].m=Metric(Hmin);
173 }
174 delete [] verticeslength;
[3913]175
176 }
177 else{
[13036]178 _error_("No edges provided");
[3913]179 }
180
181 //hVertices
[5187]182 if(bamgopts->hVertices && bamgopts->hVerticesSize[0]==nbv){
[12509]183 if(verbose>5) _printLine_(" processing hVertices");
[3913]184 for (i=0;i< nbv;i++){
[13255]185 if (!xIsNan<IssmPDouble>(bamgopts->hVertices[i])){
[5187]186 vertices[i].m=Metric((double)bamgopts->hVertices[i]);
[3913]187 }
188 }
189 }
190
191 //MetricVertices
[5187]192 if(bamgopts->metric && bamgopts->metric[0]==nbv){
[12509]193 if(verbose>5) _printLine_(" processing MetricVertices");
[3913]194 for (i=0;i< nbv;i++) {
[5187]195 vertices[i].m = Metric((double)bamgopts->metric[i*3+0],(double)bamgopts->metric[i*3+1],(double)bamgopts->metric[i*3+2]);
[3913]196 }
197 }
198
[5091]199 //MaxCornerAngle
200 if (bamgopts->MaxCornerAngle){
[12509]201 if(verbose>5) _printLine_(" processing MaxCornerAngle");
[5091]202 MaxCornerAngle=bamgopts->MaxCornerAngle*Pi/180;
[3913]203 }
204
205 //TangentAtEdges
206 if (bamggeom->TangentAtEdges){
[12509]207 if(verbose>5) _printString_(" processing TangentAtEdges");
[13036]208 if (bamggeom->TangentAtEdgesSize[1]!=4) _error_("TangentAtEdges should have 4 columns");
[3913]209 int n,i,j,k;
210 R2 tg;
211
212 n=bamggeom->TangentAtEdgesSize[0];
213 for (k=0;k<n;k++) {
214 i=(int)bamggeom->TangentAtEdges[k*4+0]-1; //for C indexing
215 j=(int)bamggeom->TangentAtEdges[k*4+1]-1; //for C indexing
216 tg.x=bamggeom->TangentAtEdges[k*4+2];
217 tg.y=bamggeom->TangentAtEdges[k*4+3];
[13036]218 if (i<0 || i>=nbe) _error_("TangentAtEdges first index exceeds matrix dimension");
219 if (j!=0 && j!=1) _error_("TangentAtEdges second index should be 1 or 2 only");
[3913]220 edges[i].tg[j] = tg;
221 }
222 }
223
224 //Corners
225 if(bamggeom->Corners){
[12509]226 if(verbose>5) _printString_(" processing Corners");
[13036]227 if (bamggeom->CornersSize[1]!=1) _error_("Corners should have 1 column");
[3913]228 n=bamggeom->CornersSize[0];
229 for (i=0;i<n;i++) {
230 j=(int)bamggeom->Corners[i]-1; //for C indexing
[13036]231 if (j>nbv-1 || j<0) _error_("Bad corner definition: should in [0 " << nbv << "]");
[5143]232 /*Required => at the same time SetRequired and SetCorner*/
[3913]233 vertices[j].SetCorner();
[5143]234 vertices[j].SetRequired();
235 }
[3913]236 }
237
238 //RequiredVertices
239 if(bamggeom->RequiredVertices){
[12509]240 if(verbose>5) _printLine_(" processing RequiredVertices");
[13036]241 if (bamggeom->RequiredVerticesSize[1]!=1) _error_("RequiredVertices should have 1 column");
[3913]242 n=bamggeom->RequiredVerticesSize[0];
243 for (i=0;i<n;i++) {
244 j=(int)bamggeom->RequiredVertices[i]-1; //for C indexing
[13036]245 if (j>nbv-1 || j<0) _error_("Bad RequiredVerticess definition: should in [0 " << nbv << "]");
[5143]246 vertices[j].SetRequired();
247 }
[3913]248 }
249
250 //RequiredEdges
251 if(bamggeom->RequiredEdges){
[12509]252 if(verbose>5) _printLine_(" processing RequiredEdges");
[13036]253 if (bamggeom->RequiredEdgesSize[1]!=1) _error_("RequiredEdges should have 1 column");
[3913]254 n=bamggeom->RequiredEdgesSize[0];
255 for (i=0;i<n;i++) {
256 j=(int)bamggeom->RequiredEdges[i]-1; //for C indexing
[13036]257 if (j>nbe-1 || j<0) _error_("Bad RequiredEdges definition: should in [0 " << nbe << "]");
[3913]258 edges[j].SetRequired();
259 }
260 }
261
262 //SubDomain
263 if(bamggeom->SubDomains){
[12509]264 if(verbose>5) _printLine_(" processing SubDomains");
[13036]265 if (bamggeom->SubDomainsSize[1]!=4) _error_("SubDomains should have 4 columns");
[5148]266 nbsubdomains=bamggeom->SubDomainsSize[0];
[5573]267 subdomains = new GeomSubDomain[nbsubdomains];
[5148]268 for (i=0;i<nbsubdomains;i++){
[3913]269 i0=(int)bamggeom->SubDomains[i*4+0];
270 i1=(int)bamggeom->SubDomains[i*4+1];
271 i2=(int)bamggeom->SubDomains[i*4+2];
272 i3=(int)bamggeom->SubDomains[i*4+3];
[13036]273 if (i0!=2) _error_("Bad Subdomain definition: first number should be 2 (for Edges)");
274 if (i1>nbe || i1<=0) _error_("Bad Subdomain definition: second number should in [1 " << nbe << "] (edge number)");
[3913]275 subdomains[i].edge=edges + (i1-1);
[5148]276 subdomains[i].direction = (int) i2;
277 subdomains[i].ReferenceNumber = i3;
[3913]278 }
279 }
280 }
[12365]281 /*}}}*/
282 /*FUNCTION Geometry::WriteGeometry{{{*/
[3913]283 void Geometry::WriteGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts){
284
285 int verbose;
286 int nbreq=0;
287 int nbreqv=0;
288 int nbtan=0;
289 int nbcracked=0;
290 int i,count;
291
292 /*Get options*/
293 verbose=bamgopts->verbose;
294
295 /*Vertices*/
[12509]296 if(verbose>5) _printLine_(" writing Vertices");
[3913]297 bamggeom->VerticesSize[0]=nbv;
298 bamggeom->VerticesSize[1]=3;
299 if (nbv){
[12456]300 bamggeom->Vertices=xNew<double>(3*nbv);
[3913]301 for (i=0;i<nbv;i++){
302 bamggeom->Vertices[i*3+0]=vertices[i].r.x;
303 bamggeom->Vertices[i*3+1]=vertices[i].r.y;
[5148]304 bamggeom->Vertices[i*3+2]=vertices[i].GetReferenceNumber();
[3913]305
306 //update counters
307 if (vertices[i].Required()) nbreqv++;
308 }
309 }
310
311 /*Edges*/
[12509]312 if(verbose>5) _printLine_(" writing Edges");
[3913]313 bamggeom->EdgesSize[0]=nbe;
314 bamggeom->EdgesSize[1]=3;
315 if (nbe){
[12456]316 bamggeom->Edges=xNew<double>(3*nbe);
[3913]317 for (i=0;i<nbe;i++){
[5149]318 bamggeom->Edges[i*3+0]=GetId(edges[i][0])+1; //back to Matlab indexing
319 bamggeom->Edges[i*3+1]=GetId(edges[i][1])+1; //back to Matlab indexing
[5148]320 bamggeom->Edges[i*3+2]=(double)edges[i].ReferenceNumber;
[3913]321
322 //update counters
323 if (edges[i].Required()) nbreq++;
324 if (edges[i].TgA() && edges[i][0].Corner()) nbtan++;
325 if (edges[i].TgB() && edges[i][1].Corner()) nbtan++;
326 }
327 }
328
329 /*RequiredEdges*/
[12509]330 if(verbose>5) _printLine_(" writing " << nbreq << " RequiredEdges");
[3913]331 bamggeom->RequiredEdgesSize[0]=nbreq;
332 bamggeom->RequiredEdgesSize[1]=1;
333 if (nbreq){
[12456]334 bamggeom->RequiredEdges=xNew<double>(1*nbreq);
[3913]335 count=0;
336 for (i=0;i<nbe;i++){
337 if (edges[i].Required()){
338 bamggeom->RequiredEdges[count]=i+1; //back to Matlab indexing
339 count=count+1;
340 }
341 }
342 }
343
344 //No corners
345
346 /*RequiredVertices*/
[12509]347 if(verbose>5) _printLine_(" writing " << nbreqv << " RequiredVertices");
[3913]348 bamggeom->RequiredVerticesSize[0]=nbreqv;
349 bamggeom->RequiredVerticesSize[1]=1;
350 if (nbreqv){
[12456]351 bamggeom->RequiredVertices=xNew<double>(1*nbreqv);
[3913]352 count=0;
353 for (i=0;i<nbv;i++){
354 if (vertices[i].Required()){
355 bamggeom->RequiredVertices[count]=i+1; //back to Matlab indexing
356 count=count+1;
357 }
358 }
359 }
360
361 /*SubDomains*/
[12509]362 if(verbose>5) _printLine_(" writing SubDomains");
[5148]363 bamggeom->SubDomainsSize[0]=nbsubdomains;
[3913]364 bamggeom->SubDomainsSize[1]=4;
[5148]365 if (nbsubdomains){
[12456]366 bamggeom->SubDomains=xNew<double>(4*nbsubdomains);
[5148]367 for (i=0;i<nbsubdomains;i++){
[3913]368 bamggeom->SubDomains[4*i+0]=2;
[5149]369 bamggeom->SubDomains[4*i+1]=GetId(subdomains[i].edge)+1; //back to Matlab indexing
[5148]370 bamggeom->SubDomains[4*i+2]=subdomains[i].direction;
371 bamggeom->SubDomains[4*i+3]=subdomains[i].ReferenceNumber;
[3913]372 }
373 }
374
375 /*TangentAtEdges*/
[12509]376 if(verbose>5) _printLine_(" writing TangentAtEdges");
[3913]377 bamggeom->TangentAtEdgesSize[0]=nbtan;
378 bamggeom->TangentAtEdgesSize[1]=4;
379 if (nbtan){
[12456]380 bamggeom->TangentAtEdges=xNew<double>(4*nbtan);
[3913]381 count=0;
382 for (i=0;i<nbe;i++){
383 if (edges[i].TgA() && edges[i][0].Corner()){
384 bamggeom->TangentAtEdges[4*i+0]=i+1; //back to Matlab indexing
385 bamggeom->TangentAtEdges[4*i+1]=1;
386 bamggeom->TangentAtEdges[4*i+2]=edges[i].tg[0].x;
387 bamggeom->TangentAtEdges[4*i+3]=edges[i].tg[0].y;
388 }
389 if (edges[i].TgB() && edges[i][1].Corner()){
390 bamggeom->TangentAtEdges[4*i+0]=i+1; //back to Matlab indexing
391 bamggeom->TangentAtEdges[4*i+1]=2;
392 bamggeom->TangentAtEdges[4*i+2]=edges[i].tg[1].x;
[5016]393 bamggeom->TangentAtEdges[4*i+3]=edges[i].tg[1].y;
[3913]394 }
395 count=count+1;
396 }
397 }
398 }
[12365]399 /*}}}*/
[3913]400
401 /*Methods*/
[12365]402 /*FUNCTION Geometry::Echo {{{*/
[5400]403 void Geometry::Echo(void){
404
[12509]405 _printLine_("Geometry:");
406 _printLine_(" nbv (number of vertices) : " << nbv);
407 _printLine_(" nbe (number of edges) : " << nbe);
408 _printLine_(" nbsubdomains: " << nbsubdomains);
409 _printLine_(" nbcurves: " << nbcurves);
410 _printLine_(" vertices: " << vertices);
411 _printLine_(" edges: " << edges);
412 _printLine_(" quadtree: " << quadtree);
413 _printLine_(" subdomains: " << subdomains);
414 _printLine_(" curves: " << curves);
415 _printLine_(" pmin (x,y): (" << pmin.x << " " << pmin.y << ")");
416 _printLine_(" pmax (x,y): (" << pmax.x << " " << pmax.y << ")");
417 _printLine_(" coefIcoor: " << coefIcoor);
418 _printLine_(" MaxCornerAngle: " << MaxCornerAngle);
[5400]419
420 return;
421 }
422 /*}}}*/
[12365]423 /*FUNCTION Geometry::Init{{{*/
[5400]424 void Geometry::Init(void){
425 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/EmptyGeometry)*/
426
427 NbRef=0;
428 nbv=0;
429 nbe=0;
430 quadtree=NULL;
431 curves=NULL;
432 edges=NULL;
433 vertices=NULL;
434 nbsubdomains=0;
435 nbcurves=0;
436 subdomains=NULL;
437 MaxCornerAngle = 10*Pi/180; //default is 10 degres
438 }
[12365]439 /*}}}*/
440 /*FUNCTION Geometry::MinimalHmin{{{*/
[5400]441 double Geometry::MinimalHmin() {
442 /* coeffIcoor = (2^30-1)/D
443 * We cannot go beyond hmin = D/2^30 because of the quadtree
444 * hmin is therefore approximately 2/coeffIcoor */
445 return 2.0/coefIcoor;
446 }/*}}}*/
[12365]447 /*FUNCTION Geometry::MaximalHmax{{{*/
[5400]448 double Geometry::MaximalHmax() {
449 return Max(pmax.x-pmin.x,pmax.y-pmin.y);
450 }/*}}}*/
[12365]451 /*FUNCTION Geometry::GetId(const GeomVertex &t){{{*/
[5573]452 long Geometry::GetId(const GeomVertex & t) const {
[5400]453 return &t - vertices;
454 }/*}}}*/
[12365]455 /*FUNCTION Geometry::GetId(const GeomVertex * t){{{*/
[5573]456 long Geometry::GetId(const GeomVertex * t) const {
[5400]457 return t - vertices;
458 }/*}}}*/
[12365]459 /*FUNCTION Geometry::GetId(const GeomEdge & t){{{*/
[5573]460 long Geometry::GetId(const GeomEdge & t) const {
[5400]461 return &t - edges;
462 }/*}}}*/
[12365]463 /*FUNCTION Geometry::GetId(const GeomEdge * t){{{*/
[5573]464 long Geometry::GetId(const GeomEdge * t) const {
[5400]465 return t - edges;
466 }/*}}}*/
[12365]467 /*FUNCTION Geometry::GetId(const Curve * c){{{*/
[5400]468 long Geometry::GetId(const Curve * c) const {
469 return c - curves;
470 }/*}}}*/
[12365]471 /*FUNCTION Geometry::Containing{{{*/
[5573]472 GeomEdge* Geometry::Containing(const R2 P, GeomEdge * start) const {
[5400]473 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Contening)*/
474
[5573]475 GeomEdge* on =start,* pon=0;
[5400]476 // walk with the cos on geometry
477 int counter=0;
478 while(pon != on){
479 counter++;
[6412]480 _assert_(counter<100);
[5400]481 pon = on;
482 R2 A= (*on)[0];
483 R2 B= (*on)[1];
484 R2 AB = B-A;
485 R2 AP = P-A;
486 R2 BP = P-B;
487 if ( (AB,AP) < 0)
488 on = on->Adj[0];
489 else if ( (AB,BP) > 0)
490 on = on->Adj[1];
491 else
492 return on;
493 }
494 return on;
495 }
[12365]496 /*}}}*/
497 /*FUNCTION Geometry::PostRead{{{*/
[5397]498 void Geometry::PostRead(){
[3913]499 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/AfterRead)*/
500
[5177]501 long i,j,k;
502 int jj;
503 long *head_v = new long[nbv];
504 long *next_p = new long[2*nbe];
505 float *eangle = new float[nbe];
506 double eps = 1e-20;
[12218]507 BamgQuadtree quadtree; // build quadtree to find duplicates
[5177]508 BamgVertex *v0 = vertices;
[5573]509 GeomVertex *v0g = (GeomVertex*) (void*)v0;
[3913]510
511 k=0;
512
513 //build quadtree for this geometry
514 for (i=0;i<nbv;i++){
515
516 /*build integer coordinates (non unique)
[5400]517 these coordinates are used by the quadtree to group
518 the vertices by groups of 5:
519 All the coordinates are transformed to ]0,1[^2
520 then, the integer coordinates are computed using
521 the transformation ]0,1[^2 -> [0 2^30-1[^2 for a quadtree of depth 30*/
[5180]522 vertices[i].i=R2ToI2(vertices[i].r);
[3913]523
[5177]524 /*find nearest vertex already present in the quadtree (NULL if empty)*/
[5120]525 BamgVertex* v=quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y);
[3913]526
[5177]527 /*if there is a vertex found that is to close to vertices[i] -> error*/
528 if( v && Norme1(v->r - vertices[i].r) < eps ){
[12509]529 _printLine_("reference numbers: " << v->ReferenceNumber << " " << vertices[i].ReferenceNumber);
530 _printLine_("Id: " << i+1);
[13385]531 _printLine_("Coords: ["<<v->r.x<<" "<<v->r.y<<"] ["<<vertices[i].r.x<<" "<<vertices[i].r.y<<"]");
532
[9306]533 delete [] next_p;
534 delete [] head_v;
535 delete [] eangle;
[13036]536 _error_("two points of the geometry are very closed to each other (see reference numbers above)");
[3913]537 }
538
[5177]539 /*Add vertices[i] to the quadtree*/
[5124]540 quadtree.Add(vertices[i]);
[3913]541 }
542
543 /* Here we use a powerful chaining algorithm
544 *
545 * 1. What is a chaining algorithm?
546 *
547 * If F is a function that goes from i in [0 n] to j in [0 m]
548 * and we want to compute the reciprocal function F-1 of F
549 * (what are the antecedents of a given j in Im(F) )
550 * We use 2 lists:
551 * head_F[j] that holds the head of lists
552 * next_F[i] that holds the list of elements that have the same image
553 *
554 * Example:
555 * i1, i2, ..., ip in [0,n] are all antecedents of a given j in [0 m]
556 * head_F[j] = ip
557 * next_F[ip]= ip-1
558 * ....
559 * next_F[i2]= i1
560 * next_F[i1]= -1 //end of the list
561 *
562 * Algorithm:
563 * for(j=0;j<m;j++) head_F[j] = -1 //initialization
564 * for(i=0;i<n;i++){
565 * j=F[i];
566 * next_F[i]= head_F[j];
567 * head_F[j]=i;
568 * }
569 *
570 * Then, we can go through all the elements that have for image j:
571 * for(i=head_F[j]; i!=-1; i=next_F[i])
572 * initialization of i by i=head_F[j]
573 * stop the loop when i=-1 (end of the chain)
574 * iterate using i=next_F[i] (next element that have for image j)
575 *
576 * 2. How to use this algorithm here?
577 *
578 * Here F is a function that associates two vertices v0 and v1 for a given edge E
579 * We want to build the reciprocal function: what are the edges that contains
580 * a vertex v?
581 * To do so, we use the same chaining algorithm but there is a difficulty
[5339]582 * coming from the fact that for F we have a couple of vertices and not one
583 * vertex.
584 * To overcome this difficulty, we use a global indexing exactly like in
[3913]585 * C/C++ so that
586 * a member of a 2-column-table can be described by one index p=i*2+j
587 * i=(int)p/2 line number of p
[5339]588 * j=p%2 column number of p
[3913]589 *
590 * Algorithm:
591 * for(i=0;i<nbv;i++) head_v[i] = -1 //initialization
592 * for(i=0;i<nbe;i++){
593 * for(j=0;j<2;j++){
594 * p=2*i+j;
595 * v=edges(i,j);
596 * next_p[p]= head_v[v];
597 * head_v[v]=p;
598 * }
599 * }
600 */
601
602 //initialize head_v as -1
603 for (i=0;i<nbv;i++) head_v[i]=-1;
604 k=0;
605 for (i=0;i<nbe;i++) {
606 //compute vector of edge i that goes from vertex 0 to vertex 1
607 R2 v10=edges[i].v[1]->r - edges[i].v[0]->r;
608 double lv10=Norme2(v10);
609 //check that its length is not 0
[9309]610 if(lv10==0){
611 delete [] next_p;
612 delete [] head_v;
613 delete [] eangle;
[13036]614 _error_("Length of edge " << i << " is 0");
[9309]615 }
[3913]616 //compute angle in [-Pi Pi]
617 eangle[i] = atan2(v10.y,v10.x);
618 //build chains head_v and next_p
619 for (j=0;j<2;j++){
[5149]620 long v=GetId(edges[i].v[j]);
[3913]621 next_p[k]=head_v[v];
622 head_v[v]=k++; //post increment: head_v[v]=k; and then k=k+1;
623 }
624 }
625
626 //sort head_v by order of increasing edges angle
627 for (i=0;i<nbv;i++) {
[5339]628 int exch=1,ord=0;
[3913]629
630 //exchange vertices position in head_v and next_p till tey are sorted
631 while (exch){
[5339]632 long *p=head_v+i;
633 long *po=p;
634 long n=*p;
[3913]635 register float angleold=-1000 ; // angle = - infinity
636 ord=0; exch=0;
637
[5339]638 // loop over the edges that contain the vertex i (till -1)
[3913]639 while (n >=0){
640 ord++;
[5339]641 long i1=n/2; // i1 = floor (n/2) -> Edge number
642 long j1=n%2; // j1 = 1 if n is odd -> Vertex index for this edge (0 or 1)
643 long* pn=next_p+n;
[3913]644
[5339]645 //Next vertex index
[3913]646 n=*pn;
647
648 //compute angle between horizontal axis and v0->v1
649 float angle = j1 ? OppositeAngle(eangle[i1]): eangle[i1];
650
651 //exchange if the current edge angle is smaller than the previous one
652 if (angleold > angle){
653 exch=1;
654 *pn=*po; // next_p[n] = n + 1
655 *po=*p; //
656 *p=n; // next_p[n+1] = n
657 po=pn; // po now point toward pn (invert next and current)
658 }
659
660 //else, continue going to the next edge position
661 else{ // to have : po -> p -> pn
662 angleold=angle; // update maximum angle
663 po=p; // po now point toward p (current position)
664 p=pn; // p now point toward pn (next position)
665 }
666 }
667 }
668
[5339]669 // angular test on current vertex to guess whether it is a corner (ord = number of edges holding i)
670 if(ord==2) {
[3913]671 long n1 = head_v[i];
672 long n2 = next_p[n1];
673 long i1 = n1/2, i2 = n2/2; // edge number
674 long j1 = n1%2, j2 = n2%2; // vertex in the edge
675 float angle1= j1 ? OppositeAngle(eangle[i1]) : eangle[i1];
676 float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2];
677 float da12 = Abs(angle2-angle1);
[5091]678 if (( da12 >= MaxCornerAngle ) && (da12 <= 2*Pi -MaxCornerAngle)) {
[3913]679 vertices[i].SetCorner() ;
680 }
[5339]681 // if the edge type/referencenumber a changing then is SetRequired();
[5146]682 if (edges[i1].type != edges[i2].type || edges[i1].Required()){
[3913]683 vertices[i].SetRequired();
684 }
[5148]685 if (edges[i1].ReferenceNumber != edges[i2].ReferenceNumber) {
[3913]686 vertices[i].SetRequired();
687 }
688 }
689 if(ord != 2) {
690 vertices[i].SetCorner();
691 }
692
[5339]693 /*close the list around the vertex to have a circular loop*/
[3913]694 long no=-1, ne = head_v[i];
695 while (ne >=0) ne = next_p[no=ne];
696 if(no>=0) next_p[no] = head_v[i];
697 }
698
[5339]699 /*Check that the list makes sense (we have all the time the same vertex)
700 * and build adjacent edges*/
[3913]701 k =0;
702 for (i=0;i<nbe;i++){
703 for (j=0;j<2;j++){
[5339]704
[3913]705 long n1 = next_p[k++];
706 long i1 = n1/2 ,j1=n1%2;
[5339]707
[13036]708 if( edges[i1].v[j1] != edges[i].v[j]) _error_("Problem while processing edges: check the edge list");
[5339]709
[3913]710 edges[i1].Adj[j1] = edges + i;
[5401]711 edges[i1].AdjVertexIndex[j1] = j;
[3913]712 }
713 }
714
[5397]715 /* generation of all the tangents*/
[3913]716 for (i=0;i<nbe;i++) {
[5339]717 R2 AB =edges[i].v[1]->r -edges[i].v[0]->r;// AB = vertex0 -> vertex1
718 double lAB=Norme2(AB); // Get length of AB
719 double ltg2[2]={0.0}; // initialize tangent
[3913]720
721 //loop over the 2 vertices of the edge
[5339]722 for (j=0;j<2;j++) {
723 R2 tg =edges[i].tg[j];
[3913]724 double ltg=Norme2(tg);
725
726 //by default, tangent=[0 0]
[5339]727 if(ltg==0){
[3913]728 //if the current vertex of the edge is not a corner
[5339]729 if(!edges[i].v[j]->Corner()){
730 /*The tangent is set as the vector between the
731 * previous and next vertices connected to current vertex
732 * normed by the edge length*/
[5401]733 tg = edges[i].v[1-j]->r - edges[i].Adj[j]->v[1-edges[i].AdjVertexIndex[j]]->r;
[3913]734 ltg= Norme2(tg);
735 tg = tg *(lAB/ltg);
736 ltg= lAB;
737 }
[5339]738 //else: a Corner no tangent => nothing to do
[3913]739 }
740 else{
741 //tangent has already been computed
742 tg = tg*(lAB/ltg),ltg=lAB;
743 }
744
[5339]745 ltg2[j] = ltg;
[3913]746
[5339]747 if ((tg,AB)<0) tg = -tg;
[3913]748
[5339]749 edges[i].tg[j]=tg;
[3913]750 }
751 if (ltg2[0]!=0) edges[i].SetTgA();
752 if (ltg2[1]!=0) edges[i].SetTgB();
753 }
754
[5397]755 /* generation of all curves (from corner to corner)*/
756 /*We proceed in 2 steps: first allocate, second build*/
[3913]757 for (int step=0;step<2;step++){
[5339]758
[5397]759 //unmark all edges
[3913]760 for (i=0;i<nbe;i++) edges[i].SetUnMark();
[5397]761 long nb_marked_edges=0;
[5339]762
[5397]763 //initialize number of curves
[5148]764 nbcurves = 0;
[4997]765
[5397]766 for (int level=0;level<2 && nb_marked_edges!=nbe;level++){
[5339]767 for (i=0;i<nbe;i++){
[5397]768
[5573]769 GeomEdge & ei=edges[i];
[5339]770 for(j=0;j<2;j++){
[5397]771 /*If current edge ei is unmarked and (level=1 or vertex i is required (corner)):
772 * we do have the first edge of a new curve*/
[5339]773 if (!ei.Mark() && (level || ei[j].Required())) {
774 int k0=j,k1;
[5573]775 GeomEdge *e=&ei;
776 GeomVertex *a=(*e)(k0); // begin
[5339]777 if(curves){
[5400]778 curves[nbcurves].FirstEdge=e;
779 curves[nbcurves].FirstVertexIndex=k0;
[5339]780 }
781 int nee=0;
[5397]782 for(;;){
[5339]783 nee++;
784 k1 = 1-k0; // next vertex of the edge
785 e->SetMark();
[5397]786 nb_marked_edges++;
[5339]787 e->CurveNumber=nbcurves;
[5573]788 GeomVertex *b=(*e)(k1);
[5397]789
790 //break if we have reached the other end of the curve
791 if (a==b || b->Required()){
[5400]792 if(curves){
793 curves[nbcurves].LastEdge=e;
794 curves[nbcurves].LastVertexIndex=k1;
795 }
[5397]796 break;
797 }
798 //else: go to next edge (adjacent)
799 else{
[5401]800 k0 = e->AdjVertexIndex[k1];// vertex in next edge
[5397]801 e = e->Adj[k1]; // next edge
802 }
[5339]803 }
804 nbcurves++;
805 if(level) a->SetRequired();
806 }
807 }
808 }
809 }
[6412]810 _assert_(nb_marked_edges && nbe);
[5397]811 //allocate if first step
812 if(step==0) curves=new Curve[nbcurves];
[3913]813 }
814
815 /*clean up*/
[4997]816 delete [] next_p;
817 delete [] head_v;
818 delete [] eangle;
[3913]819
820 }
[12365]821 /*}}}*/
822 /*FUNCTION Geometry::ProjectOnCurve {{{*/
[5573]823 GeomEdge* Geometry::ProjectOnCurve(const Edge &e,double s,BamgVertex &V,VertexOnGeom &GV) const {
[3913]824 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/ProjectOnCurve)*/
825 /*Add a vertex on an existing geometrical edge according to the metrics of the two vertices constituting the edge*/
826
827 double save_s=s;
828 int NbTry=0;
829
830retry:
831
832 s=save_s;
[5573]833 GeomEdge* on=e.GeomEdgeHook;
[3913]834 if (!on){
[13036]835 _error_("ProjectOnCurve error message: edge provided should be on geometry");
[3913]836 }
[5573]837 if (!e[0].GeomEdgeHook || !e[1].GeomEdgeHook){
[13036]838 _error_("ProjectOnCurve error message: at least one of the vertex of the edge provided is not on geometry");
[3913]839 }
840
841 //Get the two vertices of the edge
[5120]842 const BamgVertex &v0=e[0];
843 const BamgVertex &v1=e[1];
[3913]844
845 //Get position of V0, V1 and vector v0->v1
846 R2 V0=v0,V1=v1,V01=V1-V0;
847
848 //Get geometrical vertices corresponding to v0 and v1
[5573]849 VertexOnGeom vg0=*v0.GeomEdgeHook, vg1=*v1.GeomEdgeHook;
[3913]850
851 //build two pointers towrad current geometrical edge
[5573]852 GeomEdge *eg0=on, *eg1=on;
[3913]853
854 //Get edge direction and swap v0 and v1 if necessary
855 R2 Ag=(R2)(*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag;
856 int OppositeSens = (V01,AB)<0;
[5148]857 int direction0=0,direction1=1;
[3913]858 if (OppositeSens) s=1-s,Exchange(vg0,vg1),Exchange(V0,V1);
859
860 //Compute metric of new vertex as a linear interpolation of the two others
861 V.m=Metric(1.0-s,v0,s,v1);
862
863 const int mxe=100;
[5573]864 GeomEdge* ge[mxe+1];
[5148]865 int directionge[mxe+1];
[3913]866 double lge[mxe+1];
867 int bge=mxe/2,tge=bge;
[5573]868 ge[bge] = e.GeomEdgeHook;
[5148]869 directionge[bge]=1;
[3913]870
[5573]871 while (eg0!=(GeomEdge*)vg0 && (*eg0)(direction0)!=(GeomVertex*)vg0){
[3913]872 if (bge<=0) {
873 if(NbTry) {
[12509]874 _printLine_("Fatal Error: on the class Mesh before call Geometry::ProjectOnCurve");
875 _printLine_("That bug might come from:");
876 _printLine_(" 1) a mesh edge containing more than " << mxe/2 << " geometrical edges");
877 _printLine_(" 2) code bug : be sure that we call Mesh::SetVertexFieldOn() before");
878 _printLine_("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)");
[13036]879 _error_("see above");
[5130]880 }
[3913]881 NbTry++;
882 goto retry;
883 }
[5573]884 GeomEdge* tmpge = eg0;
[5148]885 ge[--bge] =eg0 = eg0->Adj[direction0];
[6412]886 _assert_(bge>=0 && bge<=mxe);
[5401]887 direction0 = 1-( directionge[bge] = tmpge->AdjVertexIndex[direction0]);
[3913]888 }
[5573]889 while (eg1 != (GeomEdge*) vg1 && (*eg1)(direction1) != (GeomVertex*) vg1) {
[3913]890 if(tge>=mxe ) {
[12509]891 _printLine_("WARNING: on the class Mesh before call Geometry::ProjectOnCurve is having issues (isn't it Eric?)");
[3913]892 NbTry++;
893 if (NbTry<2) goto retry;
[12509]894 _printLine_("Fatal Error: on the class Mesh before call Geometry::ProjectOnCurve");
895 _printLine_("That bug might come from:");
896 _printLine_(" 1) a mesh edge contening more than " << mxe/2 << " geometrical edges");
897 _printLine_(" 2) code bug : be sure that we call Mesh::SetVertexFieldOn() before");
898 _printLine_("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)");
[13036]899 _error_("see above");
[3913]900 }
[5573]901 GeomEdge* tmpge = eg1;
[5148]902 ge[++tge] =eg1 = eg1->Adj[direction1];
[5401]903 directionge[tge]= direction1 = 1-tmpge->AdjVertexIndex[direction1];
[6412]904 _assert_(tge>=0 && tge<=mxe);
[3913]905 }
906
907
[5573]908 if ((*eg0)(direction0)==(GeomVertex*)vg0)
[5148]909 vg0=VertexOnGeom(*(BamgVertex*) vg0,*eg0,direction0); //vg0 = absisce
[3913]910
[5573]911 if ((*eg1)(direction1)==(GeomVertex*)vg1)
[5148]912 vg1=VertexOnGeom(*(BamgVertex*) vg1,*eg1,direction1);
[3913]913
914 double sg;
915 if (eg0 == eg1) {
916 register double s0=vg0,s1=vg1;
917 sg = s0*(1.0-s) + s*s1;
918 on=eg0;
919 }
920 else {
921 R2 AA=V0,BB;
922 double s0,s1;
923 int i=bge;
924 double ll=0;
925 for(i=bge;i<tge;i++){
[6412]926 _assert_(i>=0 && i<=mxe);
[5148]927 BB = (*ge[i])[directionge[i]];
[3913]928 lge[i]=ll += Norme2(AA-BB);
929 AA=BB ;}
930 lge[tge]=ll+=Norme2(AA-V1);
931 // search the geometrical edge
[6412]932 _assert_(s<=1.0);
[3913]933 double ls= s*ll;
934 on =0;
935 s0 = vg0;
[5148]936 s1= directionge[bge];
[3913]937 double l0=0,l1;
938 i=bge;
939 while ( (l1=lge[i]) < ls ) {
[6412]940 _assert_(i>=0 && i<=mxe);
[5148]941 i++,s0=1-(s1=directionge[i]),l0=l1;
[3913]942 }
943 on=ge[i];
944 if (i==tge)
945 s1=vg1;
946
947 s =(ls-l0)/(l1-l0);
948 sg =s0*(1.0-s)+s*s1;
949 }
[6412]950 _assert_(on);
[3913]951 V.r= on->F(sg);
952 GV=VertexOnGeom(V,*on,sg);
953 return on;
954 }
[12365]955 /*}}}*/
956 /*FUNCTION Geometry::R2ToI2{{{*/
[5180]957 I2 Geometry::R2ToI2(const R2 & P) const {
[5177]958 /*coefIcoor is the coefficient used for integer coordinates:
959 * (x-pmin.x)
960 * Icoor x = (2^30 -1) ------------
961 * D
962 * where D is the longest side of the domain (direction x or y)
963 * so that (x-pmin.x)/D is in ]0 1[
964 *
965 * coefIcoor = (2^30 -1)/D
966 */
967 return I2( (Icoor1) (coefIcoor*(P.x-pmin.x)) ,(Icoor1) (coefIcoor*(P.y-pmin.y)) );
[5130]968 }/*}}}*/
[12365]969 /*FUNCTION Geometry::UnMarkEdges{{{*/
[5130]970 void Geometry::UnMarkEdges() {
971 for (int i=0;i<nbe;i++) edges[i].SetUnMark();
972 }/*}}}*/
[3913]973}
Note: See TracBrowser for help on using the repository browser.