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
Line 
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*/
15 /*FUNCTION Geometry::Geometry(){{{*/
16 Geometry::Geometry(){
17 Init();
18 }
19 /*}}}*/
20 /*FUNCTION Geometry::Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts){{{*/
21 Geometry::Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts){
22 Init();
23 ReadGeometry(bamggeom,bamgopts);
24 PostRead();
25 }
26 /*}}}*/
27 /*FUNCTION Geometry::Geometry(const Geometry & Gh) (COPY operator){{{*/
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;
35 vertices = nbv ? new GeomVertex[nbv] : NULL;
36 edges = nbe ? new GeomEdge[nbe]:NULL;
37 curves= nbcurves ? new Curve[nbcurves]:NULL;
38 subdomains = nbsubdomains ? new GeomSubDomain[nbsubdomains]:NULL;
39 for (i=0;i<nbe;i++)
40 edges[i].Set(Gh.edges[i],Gh,*this);
41 for (i=0;i<nbcurves;i++)
42 curves[i].Set(Gh.curves[i],Gh,*this);
43 for (i=0;i<nbsubdomains;i++)
44 subdomains[i].Set(Gh.subdomains[i],Gh,*this);
45 }
46 /*}}}*/
47 /*FUNCTION Geometry::~Geometry(){{{*/
48 Geometry::~Geometry() {
49 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/~Geometry)*/
50 if(NbRef>0){ _printString_("Trying to delete geometry and NbRef>0, probably due to an error"); return;}
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;
55 if(subdomains) delete [] subdomains;subdomains=0;
56 Init();
57 }
58 /*}}}*/
59
60 /*IO*/
61 /*FUNCTION Geometry::ReadGeometry{{{*/
62 void Geometry::ReadGeometry(BamgGeom* bamggeom,BamgOpts* bamgopts){
63
64 int verbose;
65 nbv=0;
66 nbe=0;
67 nbcurves=0;
68
69 double Hmin = HUGE_VAL;// the infinie value
70 int i,j,k,n,i0,i1,i2,i3;
71
72 /*initialize some variables*/
73 verbose= bamgopts->verbose;
74 nbv = bamggeom->VerticesSize[0];
75 nbe = bamggeom->EdgesSize[0];
76
77 //some checks
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");
80
81 //Vertices
82 if (bamggeom->Vertices){
83 if(verbose>5) _printLine_(" processing Vertices");
84 if (bamggeom->VerticesSize[1]!=3) _error_("Vertices should have 3 columns");
85 vertices = new GeomVertex[nbv];
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;
92 vertices[i].type=0;
93 }
94 /*find domain extrema (pmin,pmax) that will define the square box used for by the quadtree*/
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 }
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[
113 *
114 * coefIcoor = (2^30 -1)/D
115 */
116 coefIcoor=(MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
117 if(coefIcoor<=0) _error_("coefIcoor should be positive");
118 }
119 else{
120 _error_("No BamgVertex provided");
121 }
122
123 //Edges
124 if (bamggeom->Edges){
125 R2 zerovector(0,0);
126 double* verticeslength=NULL;
127
128 if(verbose>5) _printLine_(" processing Edges");
129 if (bamggeom->EdgesSize[1]!=3) _error_("Edges should have 3 columns");
130 edges = new GeomEdge[nbe];
131
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;
135
136 /*Loop over the edges*/
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
141 edges[i].v[0]= vertices + i1; //pointer toward vertex i1 (=&vertices[i1])
142 edges[i].v[1]= vertices + i2; //pointer toward vertex i2
143 edges[i].ReferenceNumber=(long)bamggeom->Edges[i*3+2];
144
145 //get length of edge
146 R2 x12=vertices[i2].r-vertices[i1].r;
147 double l12=sqrt((x12,x12));
148 Hmin=Min(Hmin,l12);
149
150 //initialize other fields
151 edges[i].tg[0]=zerovector;
152 edges[i].tg[1]=zerovector;
153 edges[i].AdjVertexIndex[0] = edges[i].AdjVertexIndex[1] = -1;
154 edges[i].Adj[0] = edges[i].Adj[1] = NULL;
155 edges[i].type = 0;
156
157 //Cracked?
158 if (edges[i].ReferenceNumber!=1) edges[i].SetCracked();
159
160 //prepare metric
161 vertices[i1].color++;
162 vertices[i2].color++;
163 verticeslength[i1] += l12;
164 verticeslength[i2] += l12;
165 }
166
167 // definition the default of the given mesh size
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;
175
176 }
177 else{
178 _error_("No edges provided");
179 }
180
181 //hVertices
182 if(bamgopts->hVertices && bamgopts->hVerticesSize[0]==nbv){
183 if(verbose>5) _printLine_(" processing hVertices");
184 for (i=0;i< nbv;i++){
185 if (!xIsNan<IssmPDouble>(bamgopts->hVertices[i])){
186 vertices[i].m=Metric((double)bamgopts->hVertices[i]);
187 }
188 }
189 }
190
191 //MetricVertices
192 if(bamgopts->metric && bamgopts->metric[0]==nbv){
193 if(verbose>5) _printLine_(" processing MetricVertices");
194 for (i=0;i< nbv;i++) {
195 vertices[i].m = Metric((double)bamgopts->metric[i*3+0],(double)bamgopts->metric[i*3+1],(double)bamgopts->metric[i*3+2]);
196 }
197 }
198
199 //MaxCornerAngle
200 if (bamgopts->MaxCornerAngle){
201 if(verbose>5) _printLine_(" processing MaxCornerAngle");
202 MaxCornerAngle=bamgopts->MaxCornerAngle*Pi/180;
203 }
204
205 //TangentAtEdges
206 if (bamggeom->TangentAtEdges){
207 if(verbose>5) _printString_(" processing TangentAtEdges");
208 if (bamggeom->TangentAtEdgesSize[1]!=4) _error_("TangentAtEdges should have 4 columns");
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];
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");
220 edges[i].tg[j] = tg;
221 }
222 }
223
224 //Corners
225 if(bamggeom->Corners){
226 if(verbose>5) _printString_(" processing Corners");
227 if (bamggeom->CornersSize[1]!=1) _error_("Corners should have 1 column");
228 n=bamggeom->CornersSize[0];
229 for (i=0;i<n;i++) {
230 j=(int)bamggeom->Corners[i]-1; //for C indexing
231 if (j>nbv-1 || j<0) _error_("Bad corner definition: should in [0 " << nbv << "]");
232 /*Required => at the same time SetRequired and SetCorner*/
233 vertices[j].SetCorner();
234 vertices[j].SetRequired();
235 }
236 }
237
238 //RequiredVertices
239 if(bamggeom->RequiredVertices){
240 if(verbose>5) _printLine_(" processing RequiredVertices");
241 if (bamggeom->RequiredVerticesSize[1]!=1) _error_("RequiredVertices should have 1 column");
242 n=bamggeom->RequiredVerticesSize[0];
243 for (i=0;i<n;i++) {
244 j=(int)bamggeom->RequiredVertices[i]-1; //for C indexing
245 if (j>nbv-1 || j<0) _error_("Bad RequiredVerticess definition: should in [0 " << nbv << "]");
246 vertices[j].SetRequired();
247 }
248 }
249
250 //RequiredEdges
251 if(bamggeom->RequiredEdges){
252 if(verbose>5) _printLine_(" processing RequiredEdges");
253 if (bamggeom->RequiredEdgesSize[1]!=1) _error_("RequiredEdges should have 1 column");
254 n=bamggeom->RequiredEdgesSize[0];
255 for (i=0;i<n;i++) {
256 j=(int)bamggeom->RequiredEdges[i]-1; //for C indexing
257 if (j>nbe-1 || j<0) _error_("Bad RequiredEdges definition: should in [0 " << nbe << "]");
258 edges[j].SetRequired();
259 }
260 }
261
262 //SubDomain
263 if(bamggeom->SubDomains){
264 if(verbose>5) _printLine_(" processing SubDomains");
265 if (bamggeom->SubDomainsSize[1]!=4) _error_("SubDomains should have 4 columns");
266 nbsubdomains=bamggeom->SubDomainsSize[0];
267 subdomains = new GeomSubDomain[nbsubdomains];
268 for (i=0;i<nbsubdomains;i++){
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];
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)");
275 subdomains[i].edge=edges + (i1-1);
276 subdomains[i].direction = (int) i2;
277 subdomains[i].ReferenceNumber = i3;
278 }
279 }
280 }
281 /*}}}*/
282 /*FUNCTION Geometry::WriteGeometry{{{*/
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*/
296 if(verbose>5) _printLine_(" writing Vertices");
297 bamggeom->VerticesSize[0]=nbv;
298 bamggeom->VerticesSize[1]=3;
299 if (nbv){
300 bamggeom->Vertices=xNew<double>(3*nbv);
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;
304 bamggeom->Vertices[i*3+2]=vertices[i].GetReferenceNumber();
305
306 //update counters
307 if (vertices[i].Required()) nbreqv++;
308 }
309 }
310
311 /*Edges*/
312 if(verbose>5) _printLine_(" writing Edges");
313 bamggeom->EdgesSize[0]=nbe;
314 bamggeom->EdgesSize[1]=3;
315 if (nbe){
316 bamggeom->Edges=xNew<double>(3*nbe);
317 for (i=0;i<nbe;i++){
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
320 bamggeom->Edges[i*3+2]=(double)edges[i].ReferenceNumber;
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*/
330 if(verbose>5) _printLine_(" writing " << nbreq << " RequiredEdges");
331 bamggeom->RequiredEdgesSize[0]=nbreq;
332 bamggeom->RequiredEdgesSize[1]=1;
333 if (nbreq){
334 bamggeom->RequiredEdges=xNew<double>(1*nbreq);
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*/
347 if(verbose>5) _printLine_(" writing " << nbreqv << " RequiredVertices");
348 bamggeom->RequiredVerticesSize[0]=nbreqv;
349 bamggeom->RequiredVerticesSize[1]=1;
350 if (nbreqv){
351 bamggeom->RequiredVertices=xNew<double>(1*nbreqv);
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*/
362 if(verbose>5) _printLine_(" writing SubDomains");
363 bamggeom->SubDomainsSize[0]=nbsubdomains;
364 bamggeom->SubDomainsSize[1]=4;
365 if (nbsubdomains){
366 bamggeom->SubDomains=xNew<double>(4*nbsubdomains);
367 for (i=0;i<nbsubdomains;i++){
368 bamggeom->SubDomains[4*i+0]=2;
369 bamggeom->SubDomains[4*i+1]=GetId(subdomains[i].edge)+1; //back to Matlab indexing
370 bamggeom->SubDomains[4*i+2]=subdomains[i].direction;
371 bamggeom->SubDomains[4*i+3]=subdomains[i].ReferenceNumber;
372 }
373 }
374
375 /*TangentAtEdges*/
376 if(verbose>5) _printLine_(" writing TangentAtEdges");
377 bamggeom->TangentAtEdgesSize[0]=nbtan;
378 bamggeom->TangentAtEdgesSize[1]=4;
379 if (nbtan){
380 bamggeom->TangentAtEdges=xNew<double>(4*nbtan);
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;
393 bamggeom->TangentAtEdges[4*i+3]=edges[i].tg[1].y;
394 }
395 count=count+1;
396 }
397 }
398 }
399 /*}}}*/
400
401 /*Methods*/
402 /*FUNCTION Geometry::Echo {{{*/
403 void Geometry::Echo(void){
404
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);
419
420 return;
421 }
422 /*}}}*/
423 /*FUNCTION Geometry::Init{{{*/
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 }
439 /*}}}*/
440 /*FUNCTION Geometry::MinimalHmin{{{*/
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 }/*}}}*/
447 /*FUNCTION Geometry::MaximalHmax{{{*/
448 double Geometry::MaximalHmax() {
449 return Max(pmax.x-pmin.x,pmax.y-pmin.y);
450 }/*}}}*/
451 /*FUNCTION Geometry::GetId(const GeomVertex &t){{{*/
452 long Geometry::GetId(const GeomVertex & t) const {
453 return &t - vertices;
454 }/*}}}*/
455 /*FUNCTION Geometry::GetId(const GeomVertex * t){{{*/
456 long Geometry::GetId(const GeomVertex * t) const {
457 return t - vertices;
458 }/*}}}*/
459 /*FUNCTION Geometry::GetId(const GeomEdge & t){{{*/
460 long Geometry::GetId(const GeomEdge & t) const {
461 return &t - edges;
462 }/*}}}*/
463 /*FUNCTION Geometry::GetId(const GeomEdge * t){{{*/
464 long Geometry::GetId(const GeomEdge * t) const {
465 return t - edges;
466 }/*}}}*/
467 /*FUNCTION Geometry::GetId(const Curve * c){{{*/
468 long Geometry::GetId(const Curve * c) const {
469 return c - curves;
470 }/*}}}*/
471 /*FUNCTION Geometry::Containing{{{*/
472 GeomEdge* Geometry::Containing(const R2 P, GeomEdge * start) const {
473 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Contening)*/
474
475 GeomEdge* on =start,* pon=0;
476 // walk with the cos on geometry
477 int counter=0;
478 while(pon != on){
479 counter++;
480 _assert_(counter<100);
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 }
496 /*}}}*/
497 /*FUNCTION Geometry::PostRead{{{*/
498 void Geometry::PostRead(){
499 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/AfterRead)*/
500
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;
507 BamgQuadtree quadtree; // build quadtree to find duplicates
508 BamgVertex *v0 = vertices;
509 GeomVertex *v0g = (GeomVertex*) (void*)v0;
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)
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*/
522 vertices[i].i=R2ToI2(vertices[i].r);
523
524 /*find nearest vertex already present in the quadtree (NULL if empty)*/
525 BamgVertex* v=quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y);
526
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 ){
529 _printLine_("reference numbers: " << v->ReferenceNumber << " " << vertices[i].ReferenceNumber);
530 _printLine_("Id: " << i+1);
531 _printLine_("Coords: ["<<v->r.x<<" "<<v->r.y<<"] ["<<vertices[i].r.x<<" "<<vertices[i].r.y<<"]");
532
533 delete [] next_p;
534 delete [] head_v;
535 delete [] eangle;
536 _error_("two points of the geometry are very closed to each other (see reference numbers above)");
537 }
538
539 /*Add vertices[i] to the quadtree*/
540 quadtree.Add(vertices[i]);
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
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
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
588 * j=p%2 column number of p
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
610 if(lv10==0){
611 delete [] next_p;
612 delete [] head_v;
613 delete [] eangle;
614 _error_("Length of edge " << i << " is 0");
615 }
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++){
620 long v=GetId(edges[i].v[j]);
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++) {
628 int exch=1,ord=0;
629
630 //exchange vertices position in head_v and next_p till tey are sorted
631 while (exch){
632 long *p=head_v+i;
633 long *po=p;
634 long n=*p;
635 register float angleold=-1000 ; // angle = - infinity
636 ord=0; exch=0;
637
638 // loop over the edges that contain the vertex i (till -1)
639 while (n >=0){
640 ord++;
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;
644
645 //Next vertex index
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
669 // angular test on current vertex to guess whether it is a corner (ord = number of edges holding i)
670 if(ord==2) {
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);
678 if (( da12 >= MaxCornerAngle ) && (da12 <= 2*Pi -MaxCornerAngle)) {
679 vertices[i].SetCorner() ;
680 }
681 // if the edge type/referencenumber a changing then is SetRequired();
682 if (edges[i1].type != edges[i2].type || edges[i1].Required()){
683 vertices[i].SetRequired();
684 }
685 if (edges[i1].ReferenceNumber != edges[i2].ReferenceNumber) {
686 vertices[i].SetRequired();
687 }
688 }
689 if(ord != 2) {
690 vertices[i].SetCorner();
691 }
692
693 /*close the list around the vertex to have a circular loop*/
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
699 /*Check that the list makes sense (we have all the time the same vertex)
700 * and build adjacent edges*/
701 k =0;
702 for (i=0;i<nbe;i++){
703 for (j=0;j<2;j++){
704
705 long n1 = next_p[k++];
706 long i1 = n1/2 ,j1=n1%2;
707
708 if( edges[i1].v[j1] != edges[i].v[j]) _error_("Problem while processing edges: check the edge list");
709
710 edges[i1].Adj[j1] = edges + i;
711 edges[i1].AdjVertexIndex[j1] = j;
712 }
713 }
714
715 /* generation of all the tangents*/
716 for (i=0;i<nbe;i++) {
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
720
721 //loop over the 2 vertices of the edge
722 for (j=0;j<2;j++) {
723 R2 tg =edges[i].tg[j];
724 double ltg=Norme2(tg);
725
726 //by default, tangent=[0 0]
727 if(ltg==0){
728 //if the current vertex of the edge is not a corner
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*/
733 tg = edges[i].v[1-j]->r - edges[i].Adj[j]->v[1-edges[i].AdjVertexIndex[j]]->r;
734 ltg= Norme2(tg);
735 tg = tg *(lAB/ltg);
736 ltg= lAB;
737 }
738 //else: a Corner no tangent => nothing to do
739 }
740 else{
741 //tangent has already been computed
742 tg = tg*(lAB/ltg),ltg=lAB;
743 }
744
745 ltg2[j] = ltg;
746
747 if ((tg,AB)<0) tg = -tg;
748
749 edges[i].tg[j]=tg;
750 }
751 if (ltg2[0]!=0) edges[i].SetTgA();
752 if (ltg2[1]!=0) edges[i].SetTgB();
753 }
754
755 /* generation of all curves (from corner to corner)*/
756 /*We proceed in 2 steps: first allocate, second build*/
757 for (int step=0;step<2;step++){
758
759 //unmark all edges
760 for (i=0;i<nbe;i++) edges[i].SetUnMark();
761 long nb_marked_edges=0;
762
763 //initialize number of curves
764 nbcurves = 0;
765
766 for (int level=0;level<2 && nb_marked_edges!=nbe;level++){
767 for (i=0;i<nbe;i++){
768
769 GeomEdge & ei=edges[i];
770 for(j=0;j<2;j++){
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*/
773 if (!ei.Mark() && (level || ei[j].Required())) {
774 int k0=j,k1;
775 GeomEdge *e=&ei;
776 GeomVertex *a=(*e)(k0); // begin
777 if(curves){
778 curves[nbcurves].FirstEdge=e;
779 curves[nbcurves].FirstVertexIndex=k0;
780 }
781 int nee=0;
782 for(;;){
783 nee++;
784 k1 = 1-k0; // next vertex of the edge
785 e->SetMark();
786 nb_marked_edges++;
787 e->CurveNumber=nbcurves;
788 GeomVertex *b=(*e)(k1);
789
790 //break if we have reached the other end of the curve
791 if (a==b || b->Required()){
792 if(curves){
793 curves[nbcurves].LastEdge=e;
794 curves[nbcurves].LastVertexIndex=k1;
795 }
796 break;
797 }
798 //else: go to next edge (adjacent)
799 else{
800 k0 = e->AdjVertexIndex[k1];// vertex in next edge
801 e = e->Adj[k1]; // next edge
802 }
803 }
804 nbcurves++;
805 if(level) a->SetRequired();
806 }
807 }
808 }
809 }
810 _assert_(nb_marked_edges && nbe);
811 //allocate if first step
812 if(step==0) curves=new Curve[nbcurves];
813 }
814
815 /*clean up*/
816 delete [] next_p;
817 delete [] head_v;
818 delete [] eangle;
819
820 }
821 /*}}}*/
822 /*FUNCTION Geometry::ProjectOnCurve {{{*/
823 GeomEdge* Geometry::ProjectOnCurve(const Edge &e,double s,BamgVertex &V,VertexOnGeom &GV) const {
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;
833 GeomEdge* on=e.GeomEdgeHook;
834 if (!on){
835 _error_("ProjectOnCurve error message: edge provided should be on geometry");
836 }
837 if (!e[0].GeomEdgeHook || !e[1].GeomEdgeHook){
838 _error_("ProjectOnCurve error message: at least one of the vertex of the edge provided is not on geometry");
839 }
840
841 //Get the two vertices of the edge
842 const BamgVertex &v0=e[0];
843 const BamgVertex &v1=e[1];
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
849 VertexOnGeom vg0=*v0.GeomEdgeHook, vg1=*v1.GeomEdgeHook;
850
851 //build two pointers towrad current geometrical edge
852 GeomEdge *eg0=on, *eg1=on;
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;
857 int direction0=0,direction1=1;
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;
864 GeomEdge* ge[mxe+1];
865 int directionge[mxe+1];
866 double lge[mxe+1];
867 int bge=mxe/2,tge=bge;
868 ge[bge] = e.GeomEdgeHook;
869 directionge[bge]=1;
870
871 while (eg0!=(GeomEdge*)vg0 && (*eg0)(direction0)!=(GeomVertex*)vg0){
872 if (bge<=0) {
873 if(NbTry) {
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)");
879 _error_("see above");
880 }
881 NbTry++;
882 goto retry;
883 }
884 GeomEdge* tmpge = eg0;
885 ge[--bge] =eg0 = eg0->Adj[direction0];
886 _assert_(bge>=0 && bge<=mxe);
887 direction0 = 1-( directionge[bge] = tmpge->AdjVertexIndex[direction0]);
888 }
889 while (eg1 != (GeomEdge*) vg1 && (*eg1)(direction1) != (GeomVertex*) vg1) {
890 if(tge>=mxe ) {
891 _printLine_("WARNING: on the class Mesh before call Geometry::ProjectOnCurve is having issues (isn't it Eric?)");
892 NbTry++;
893 if (NbTry<2) goto retry;
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)");
899 _error_("see above");
900 }
901 GeomEdge* tmpge = eg1;
902 ge[++tge] =eg1 = eg1->Adj[direction1];
903 directionge[tge]= direction1 = 1-tmpge->AdjVertexIndex[direction1];
904 _assert_(tge>=0 && tge<=mxe);
905 }
906
907
908 if ((*eg0)(direction0)==(GeomVertex*)vg0)
909 vg0=VertexOnGeom(*(BamgVertex*) vg0,*eg0,direction0); //vg0 = absisce
910
911 if ((*eg1)(direction1)==(GeomVertex*)vg1)
912 vg1=VertexOnGeom(*(BamgVertex*) vg1,*eg1,direction1);
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++){
926 _assert_(i>=0 && i<=mxe);
927 BB = (*ge[i])[directionge[i]];
928 lge[i]=ll += Norme2(AA-BB);
929 AA=BB ;}
930 lge[tge]=ll+=Norme2(AA-V1);
931 // search the geometrical edge
932 _assert_(s<=1.0);
933 double ls= s*ll;
934 on =0;
935 s0 = vg0;
936 s1= directionge[bge];
937 double l0=0,l1;
938 i=bge;
939 while ( (l1=lge[i]) < ls ) {
940 _assert_(i>=0 && i<=mxe);
941 i++,s0=1-(s1=directionge[i]),l0=l1;
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 }
950 _assert_(on);
951 V.r= on->F(sg);
952 GV=VertexOnGeom(V,*on,sg);
953 return on;
954 }
955 /*}}}*/
956 /*FUNCTION Geometry::R2ToI2{{{*/
957 I2 Geometry::R2ToI2(const R2 & P) const {
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)) );
968 }/*}}}*/
969 /*FUNCTION Geometry::UnMarkEdges{{{*/
970 void Geometry::UnMarkEdges() {
971 for (int i=0;i<nbe;i++) edges[i].SetUnMark();
972 }/*}}}*/
973}
Note: See TracBrowser for help on using the repository browser.