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

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

merged trunk-jpl and trunk for revision 13974
File size: 30.1 KB
Line 
1#include <cstdio>
2#include <cstring>
3#include <cmath>
4#include <ctime>
5
6#include "../objects/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,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 long *head_v = new long[nbv];
503 long *next_p = new long[2 *nbe];
504 float *eangle = new float[nbe];
505 double eps = 1e-20;
506 BamgQuadtree quadtree; // build quadtree to find duplicates
507 BamgVertex *v0 = vertices;
508
509 k=0;
510
511 //build quadtree for this geometry
512 for (i=0;i<nbv;i++){
513
514 /*build integer coordinates (non unique)
515 these coordinates are used by the quadtree to group
516 the vertices by groups of 5:
517 All the coordinates are transformed to ]0,1[^2
518 then, the integer coordinates are computed using
519 the transformation ]0,1[^2 -> [0 2^30-1[^2 for a quadtree of depth 30*/
520 vertices[i].i=R2ToI2(vertices[i].r);
521
522 /*find nearest vertex already present in the quadtree (NULL if empty)*/
523 BamgVertex* v=quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y);
524
525 /*if there is a vertex found that is to close to vertices[i] -> error*/
526 if( v && Norme1(v->r - vertices[i].r) < eps ){
527 _printLine_("reference numbers: " << v->ReferenceNumber << " " << vertices[i].ReferenceNumber);
528 _printLine_("Id: " << i+1);
529 _printLine_("Coords: ["<<v->r.x<<" "<<v->r.y<<"] ["<<vertices[i].r.x<<" "<<vertices[i].r.y<<"]");
530
531 delete [] next_p;
532 delete [] head_v;
533 delete [] eangle;
534 _error_("two points of the geometry are very closed to each other (see reference numbers above)");
535 }
536
537 /*Add vertices[i] to the quadtree*/
538 quadtree.Add(vertices[i]);
539 }
540
541 /* Here we use a powerful chaining algorithm
542 *
543 * 1. What is a chaining algorithm?
544 *
545 * If F is a function that goes from i in [0 n] to j in [0 m]
546 * and we want to compute the reciprocal function F-1 of F
547 * (what are the antecedents of a given j in Im(F) )
548 * We use 2 lists:
549 * head_F[j] that holds the head of lists
550 * next_F[i] that holds the list of elements that have the same image
551 *
552 * Example:
553 * i1, i2, ..., ip in [0,n] are all antecedents of a given j in [0 m]
554 * head_F[j] = ip
555 * next_F[ip]= ip-1
556 * ....
557 * next_F[i2]= i1
558 * next_F[i1]= -1 //end of the list
559 *
560 * Algorithm:
561 * for(j=0;j<m;j++) head_F[j] = -1 //initialization
562 * for(i=0;i<n;i++){
563 * j=F[i];
564 * next_F[i]= head_F[j];
565 * head_F[j]=i;
566 * }
567 *
568 * Then, we can go through all the elements that have for image j:
569 * for(i=head_F[j]; i!=-1; i=next_F[i])
570 * initialization of i by i=head_F[j]
571 * stop the loop when i=-1 (end of the chain)
572 * iterate using i=next_F[i] (next element that have for image j)
573 *
574 * 2. How to use this algorithm here?
575 *
576 * Here F is a function that associates two vertices v0 and v1 for a given edge E
577 * We want to build the reciprocal function: what are the edges that contains
578 * a vertex v?
579 * To do so, we use the same chaining algorithm but there is a difficulty
580 * coming from the fact that for F we have a couple of vertices and not one
581 * vertex.
582 * To overcome this difficulty, we use a global indexing exactly like in
583 * C/C++ so that
584 * a member of a 2-column-table can be described by one index p=i*2+j
585 * i=(int)p/2 line number of p
586 * j=p%2 column number of p
587 *
588 * Algorithm:
589 * for(i=0;i<nbv;i++) head_v[i] = -1 //initialization
590 * for(i=0;i<nbe;i++){
591 * for(j=0;j<2;j++){
592 * p=2*i+j;
593 * v=edges(i,j);
594 * next_p[p]= head_v[v];
595 * head_v[v]=p;
596 * }
597 * }
598 */
599
600 //initialize head_v as -1
601 for (i=0;i<nbv;i++) head_v[i]=-1;
602 k=0;
603 for (i=0;i<nbe;i++) {
604 //compute vector of edge i that goes from vertex 0 to vertex 1
605 R2 v10=edges[i].v[1]->r - edges[i].v[0]->r;
606 double lv10=Norme2(v10);
607 //check that its length is not 0
608 if(lv10==0){
609 delete [] next_p;
610 delete [] head_v;
611 delete [] eangle;
612 _error_("Length of edge " << i << " is 0");
613 }
614 //compute angle in [-Pi Pi]
615 eangle[i] = atan2(v10.y,v10.x);
616 //build chains head_v and next_p
617 for (j=0;j<2;j++){
618 long v=GetId(edges[i].v[j]);
619 next_p[k]=head_v[v];
620 head_v[v]=k++; //post increment: head_v[v]=k; and then k=k+1;
621 }
622 }
623
624 //sort head_v by order of increasing edges angle
625 for (i=0;i<nbv;i++) {
626 int exch=1,ord=0;
627
628 //exchange vertices position in head_v and next_p till tey are sorted
629 while (exch){
630 long *p=head_v+i;
631 long *po=p;
632 long n=*p;
633 register float angleold=-1000 ; // angle = - infinity
634 ord=0; exch=0;
635
636 // loop over the edges that contain the vertex i (till -1)
637 while (n >=0){
638 ord++;
639 long i1=n/2; // i1 = floor (n/2) -> Edge number
640 long j1=n%2; // j1 = 1 if n is odd -> Vertex index for this edge (0 or 1)
641 long* pn=next_p+n;
642
643 //Next vertex index
644 n=*pn;
645
646 //compute angle between horizontal axis and v0->v1
647 float angle = j1 ? OppositeAngle(eangle[i1]): eangle[i1];
648
649 //exchange if the current edge angle is smaller than the previous one
650 if (angleold > angle){
651 exch=1;
652 *pn=*po; // next_p[n] = n + 1
653 *po=*p; //
654 *p=n; // next_p[n+1] = n
655 po=pn; // po now point toward pn (invert next and current)
656 }
657
658 //else, continue going to the next edge position
659 else{ // to have : po -> p -> pn
660 angleold=angle; // update maximum angle
661 po=p; // po now point toward p (current position)
662 p=pn; // p now point toward pn (next position)
663 }
664 }
665 }
666
667 // angular test on current vertex to guess whether it is a corner (ord = number of edges holding i)
668 if(ord==2) {
669 long n1 = head_v[i];
670 long n2 = next_p[n1];
671 long i1 = n1/2, i2 = n2/2; // edge number
672 long j1 = n1%2, j2 = n2%2; // vertex in the edge
673 float angle1= j1 ? OppositeAngle(eangle[i1]) : eangle[i1];
674 float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2];
675 float da12 = Abs(angle2-angle1);
676 if (( da12 >= MaxCornerAngle ) && (da12 <= 2*Pi -MaxCornerAngle)) {
677 vertices[i].SetCorner() ;
678 }
679 // if the edge type/referencenumber a changing then is SetRequired();
680 if (edges[i1].type != edges[i2].type || edges[i1].Required()){
681 vertices[i].SetRequired();
682 }
683 if (edges[i1].ReferenceNumber != edges[i2].ReferenceNumber) {
684 vertices[i].SetRequired();
685 }
686 }
687 if(ord != 2) {
688 vertices[i].SetCorner();
689 }
690
691 /*close the list around the vertex to have a circular loop*/
692 long no=-1, ne = head_v[i];
693 while (ne >=0) ne = next_p[no=ne];
694 if(no>=0) next_p[no] = head_v[i];
695 }
696
697 /*Check that the list makes sense (we have all the time the same vertex)
698 * and build adjacent edges*/
699 k =0;
700 for (i=0;i<nbe;i++){
701 for (j=0;j<2;j++){
702
703 long n1 = next_p[k++];
704 long i1 = n1/2 ,j1=n1%2;
705
706 if( edges[i1].v[j1] != edges[i].v[j]) _error_("Problem while processing edges: check the edge list");
707
708 edges[i1].Adj[j1] = edges + i;
709 edges[i1].AdjVertexIndex[j1] = j;
710 }
711 }
712
713 /* generation of all the tangents*/
714 for (i=0;i<nbe;i++) {
715 R2 AB =edges[i].v[1]->r -edges[i].v[0]->r;// AB = vertex0 -> vertex1
716 double lAB=Norme2(AB); // Get length of AB
717 double ltg2[2]={0.0}; // initialize tangent
718
719 //loop over the 2 vertices of the edge
720 for (j=0;j<2;j++) {
721 R2 tg =edges[i].tg[j];
722 double ltg=Norme2(tg);
723
724 //by default, tangent=[0 0]
725 if(ltg==0){
726 //if the current vertex of the edge is not a corner
727 if(!edges[i].v[j]->Corner()){
728 /*The tangent is set as the vector between the
729 * previous and next vertices connected to current vertex
730 * normed by the edge length*/
731 tg = edges[i].v[1-j]->r - edges[i].Adj[j]->v[1-edges[i].AdjVertexIndex[j]]->r;
732 ltg= Norme2(tg);
733 tg = tg *(lAB/ltg);
734 ltg= lAB;
735 }
736 //else: a Corner no tangent => nothing to do
737 }
738 else{
739 //tangent has already been computed
740 tg = tg*(lAB/ltg),ltg=lAB;
741 }
742
743 ltg2[j] = ltg;
744
745 if ((tg,AB)<0) tg = -tg;
746
747 edges[i].tg[j]=tg;
748 }
749 if (ltg2[0]!=0) edges[i].SetTgA();
750 if (ltg2[1]!=0) edges[i].SetTgB();
751 }
752
753 /* generation of all curves (from corner to corner)*/
754 /*We proceed in 2 steps: first allocate, second build*/
755 for (int step=0;step<2;step++){
756
757 //unmark all edges
758 for (i=0;i<nbe;i++) edges[i].SetUnMark();
759 long nb_marked_edges=0;
760
761 //initialize number of curves
762 nbcurves = 0;
763
764 for (int level=0;level<2 && nb_marked_edges!=nbe;level++){
765 for (i=0;i<nbe;i++){
766
767 GeomEdge & ei=edges[i];
768 for(j=0;j<2;j++){
769 /*If current edge ei is unmarked and (level=1 or vertex i is required (corner)):
770 * we do have the first edge of a new curve*/
771 if (!ei.Mark() && (level || ei[j].Required())) {
772 int k0=j,k1;
773 GeomEdge *e=&ei;
774 GeomVertex *a=(*e)(k0); // begin
775 if(curves){
776 curves[nbcurves].FirstEdge=e;
777 curves[nbcurves].FirstVertexIndex=k0;
778 }
779 int nee=0;
780 for(;;){
781 nee++;
782 k1 = 1-k0; // next vertex of the edge
783 e->SetMark();
784 nb_marked_edges++;
785 e->CurveNumber=nbcurves;
786 GeomVertex *b=(*e)(k1);
787
788 //break if we have reached the other end of the curve
789 if (a==b || b->Required()){
790 if(curves){
791 curves[nbcurves].LastEdge=e;
792 curves[nbcurves].LastVertexIndex=k1;
793 }
794 break;
795 }
796 //else: go to next edge (adjacent)
797 else{
798 k0 = e->AdjVertexIndex[k1];// vertex in next edge
799 e = e->Adj[k1]; // next edge
800 }
801 }
802 nbcurves++;
803 if(level) a->SetRequired();
804 }
805 }
806 }
807 }
808 _assert_(nb_marked_edges && nbe);
809 //allocate if first step
810 if(step==0) curves=new Curve[nbcurves];
811 }
812
813 /*clean up*/
814 delete [] next_p;
815 delete [] head_v;
816 delete [] eangle;
817
818 }
819 /*}}}*/
820 /*FUNCTION Geometry::ProjectOnCurve {{{*/
821 GeomEdge* Geometry::ProjectOnCurve(const Edge &e,double s,BamgVertex &V,VertexOnGeom &GV) const {
822 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/ProjectOnCurve)*/
823 /*Add a vertex on an existing geometrical edge according to the metrics of the two vertices constituting the edge*/
824
825 double save_s=s;
826 int NbTry=0;
827
828retry:
829
830 s=save_s;
831 GeomEdge* on=e.GeomEdgeHook;
832 if (!on){
833 _error_("ProjectOnCurve error message: edge provided should be on geometry");
834 }
835 if (!e[0].GeomEdgeHook || !e[1].GeomEdgeHook){
836 _error_("ProjectOnCurve error message: at least one of the vertex of the edge provided is not on geometry");
837 }
838
839 //Get the two vertices of the edge
840 const BamgVertex &v0=e[0];
841 const BamgVertex &v1=e[1];
842
843 //Get position of V0, V1 and vector v0->v1
844 R2 V0=v0,V1=v1,V01=V1-V0;
845
846 //Get geometrical vertices corresponding to v0 and v1
847 VertexOnGeom vg0=*v0.GeomEdgeHook, vg1=*v1.GeomEdgeHook;
848
849 //build two pointers towrad current geometrical edge
850 GeomEdge *eg0=on, *eg1=on;
851
852 //Get edge direction and swap v0 and v1 if necessary
853 R2 Ag=(R2)(*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag;
854 int OppositeSens = (V01,AB)<0;
855 int direction0=0,direction1=1;
856 if (OppositeSens) s=1-s,Exchange(vg0,vg1),Exchange(V0,V1);
857
858 //Compute metric of new vertex as a linear interpolation of the two others
859 V.m=Metric(1.0-s,v0,s,v1);
860
861 const int mxe=100;
862 GeomEdge* ge[mxe+1];
863 int directionge[mxe+1];
864 double lge[mxe+1];
865 int bge=mxe/2,tge=bge;
866 ge[bge] = e.GeomEdgeHook;
867 directionge[bge]=1;
868
869 while(eg0!=(GeomEdge*)vg0 && (*eg0)(direction0)!=(GeomVertex*)vg0){
870 if (bge<=0) {
871 if(NbTry) {
872 _printLine_("Fatal Error: on the class Mesh before call Geometry::ProjectOnCurve");
873 _printLine_("That bug might come from:");
874 _printLine_(" 1) a mesh edge containing more than " << mxe/2 << " geometrical edges");
875 _printLine_(" 2) code bug : be sure that we call Mesh::SetVertexFieldOn() before");
876 _printLine_("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)");
877 _error_("see above");
878 }
879 NbTry++;
880 goto retry;
881 }
882 GeomEdge* tmpge = eg0;
883 ge[--bge] =eg0 = eg0->Adj[direction0];
884 _assert_(bge>=0 && bge<=mxe);
885 direction0 = 1-( directionge[bge] = tmpge->AdjVertexIndex[direction0]);
886 }
887 while (eg1 != (GeomEdge*) vg1 && (*eg1)(direction1) != (GeomVertex*) vg1) {
888 if(tge>=mxe ) {
889 _printLine_("WARNING: on the class Mesh before call Geometry::ProjectOnCurve is having issues (isn't it Eric?)");
890 NbTry++;
891 if (NbTry<2) goto retry;
892 _printLine_("Fatal Error: on the class Mesh before call Geometry::ProjectOnCurve");
893 _printLine_("That bug might come from:");
894 _printLine_(" 1) a mesh edge contening more than " << mxe/2 << " geometrical edges");
895 _printLine_(" 2) code bug : be sure that we call Mesh::SetVertexFieldOn() before");
896 _printLine_("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)");
897 _error_("see above");
898 }
899 GeomEdge* tmpge = eg1;
900 ge[++tge] =eg1 = eg1->Adj[direction1];
901 directionge[tge]= direction1 = 1-tmpge->AdjVertexIndex[direction1];
902 _assert_(tge>=0 && tge<=mxe);
903 }
904
905 if ((*eg0)(direction0)==(GeomVertex*)vg0)
906 vg0=VertexOnGeom(*(BamgVertex*) vg0,*eg0,direction0); //vg0 = absisce
907
908 if ((*eg1)(direction1)==(GeomVertex*)vg1)
909 vg1=VertexOnGeom(*(BamgVertex*) vg1,*eg1,direction1);
910
911 double sg;
912 if (eg0 == eg1) {
913 register double s0=vg0,s1=vg1;
914 sg = s0*(1.0-s) + s*s1;
915 on=eg0;
916 }
917 else {
918 R2 AA=V0,BB;
919 double s0,s1;
920 int i=bge;
921 double ll=0;
922 for(i=bge;i<tge;i++){
923 _assert_(i>=0 && i<=mxe);
924 BB = (*ge[i])[directionge[i]];
925 lge[i]=ll += Norme2(AA-BB);
926 AA=BB ;}
927 lge[tge]=ll+=Norme2(AA-V1);
928 // search the geometrical edge
929 _assert_(s<=1.0);
930 double ls= s*ll;
931 on =0;
932 s0 = vg0;
933 s1= directionge[bge];
934 double l0=0,l1;
935 i=bge;
936 while ( (l1=lge[i]) < ls ) {
937 _assert_(i>=0 && i<=mxe);
938 i++,s0=1-(s1=directionge[i]),l0=l1;
939 }
940 on=ge[i];
941 if (i==tge)
942 s1=vg1;
943
944 s =(ls-l0)/(l1-l0);
945 sg =s0*(1.0-s)+s*s1;
946 }
947 _assert_(on);
948 V.r= on->F(sg);
949 GV=VertexOnGeom(V,*on,sg);
950 return on;
951 }
952 /*}}}*/
953 /*FUNCTION Geometry::R2ToI2{{{*/
954 I2 Geometry::R2ToI2(const R2 & P) const {
955 /*coefIcoor is the coefficient used for integer coordinates:
956 * (x-pmin.x)
957 * Icoor x = (2^30 -1) ------------
958 * D
959 * where D is the longest side of the domain (direction x or y)
960 * so that (x-pmin.x)/D is in ]0 1[
961 *
962 * coefIcoor = (2^30 -1)/D
963 */
964 return I2( (Icoor1) (coefIcoor*(P.x-pmin.x)) ,(Icoor1) (coefIcoor*(P.y-pmin.y)) );
965 }/*}}}*/
966 /*FUNCTION Geometry::UnMarkEdges{{{*/
967 void Geometry::UnMarkEdges() {
968 for (int i=0;i<nbe;i++) edges[i].SetUnMark();
969 }/*}}}*/
970}
Note: See TracBrowser for help on using the repository browser.