source: issm/trunk/src/c/Bamgx/Geometry.cpp@ 2809

Last change on this file since 2809 was 2809, checked in by Mathieu Morlighem, 15 years ago

removed Bamgx/MeshGeom.cpp
File size: 31.3 KB
Line 
1#include "../shared/shared.h"
2#include "../include/macros.h"
3#include "../toolkits/toolkits.h"
4
5#include <cstdio>
6#include <cstring>
7#include <cmath>
8#include <ctime>
9#include "Mesh2.h"
10#include "QuadTree.h"
11#include "SetOfE4.h"
12
13
14namespace bamg {
15
16 static const Direction NoDirOfSearch=Direction();
17
18 /*Constructors/Destructors*/
19
20 /*IO*/
21 /*FUNCTION Geometry::WriteGeometry{{{1*/
22 void Geometry::WriteGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts){
23
24 int verbose;
25 int nbreq=0;
26 int nbreqv=0;
27 int nbtan=0;
28 int nbcracked=0;
29 int i,count;
30
31 verbose=bamgopts->verbose;
32
33 //Vertices
34 if(verbose>3) printf(" writing Vertices\n");
35 bamggeom->NumVertices=nbv;
36 xfree((void**)&bamggeom->Vertices);
37 if (nbv){
38 bamggeom->Vertices=(double*)xmalloc(3*nbv*sizeof(double));
39 for (i=0;i<nbv;i++){
40 bamggeom->Vertices[i*3+0]=vertices[i].r.x;
41 bamggeom->Vertices[i*3+1]=vertices[i].r.y;
42 bamggeom->Vertices[i*3+2]=vertices[i].ref();
43
44 //update counters
45 if (vertices[i].Required()) nbreqv++;
46 }
47 }
48 else{
49 bamggeom->Vertices=NULL;
50 }
51
52 //Edges
53 if(verbose>3) printf(" writing Edges\n");
54 bamggeom->NumEdges=nbe;
55 xfree((void**)&bamggeom->Edges);
56 if (nbe){
57 bamggeom->Edges=(double*)xmalloc(3*nbe*sizeof(double));
58 for (i=0;i<nbe;i++){
59 bamggeom->Edges[i*3+0]=Number(edges[i][0])+1; //back to Matlab indexing
60 bamggeom->Edges[i*3+1]=Number(edges[i][1])+1; //back to Matlab indexing
61 //bamggeom->Edges[i*3+2]=(double)edges[i].ref(); //TEST does not compile???
62 bamggeom->Edges[i*3+2]=1; //TEST for now
63
64 //update counters
65 if (edges[i].Required()) nbreq++;
66 if (edges[i].Cracked()){
67 if (i<=Number(edges[i].link)) nbcracked++;
68 }
69 if (edges[i].TgA() && edges[i][0].Corner()) nbtan++;
70 if (edges[i].TgB() && edges[i][1].Corner()) nbtan++;
71 }
72 }
73 else{
74 bamggeom->Edges=NULL;
75 }
76
77 //CrackedEdges
78 if(verbose>3) printf(" writing CrackedEdges\n");
79 bamggeom->NumCrackedEdges=nbcracked;
80 xfree((void**)&bamggeom->CrackedEdges);
81 if (nbcracked){
82 bamggeom->CrackedEdges=(double*)xmalloc(2*nbcracked*sizeof(double));
83 count=0;
84 for (i=0;i<nbe;i++){
85 if (edges[i].Cracked()){
86 if (i<=Number(edges[i].link)){
87 bamggeom->CrackedEdges[count*2+0]=i+1; //back to Matlab indexing
88 bamggeom->CrackedEdges[count*2+1]=Number(edges[i].link)+1; //back to Matlab indexing
89 count=count+1;
90 }
91 }
92 }
93 }
94 else{
95 bamggeom->CrackedEdges=NULL;
96 }
97
98 //RequiredEdges
99 if(verbose>3) printf(" writing RequiredEdges\n");
100 bamggeom->NumRequiredEdges=nbreq;
101 xfree((void**)&bamggeom->RequiredEdges);
102 if (nbreq){
103 bamggeom->RequiredEdges=(double*)xmalloc(1*nbreq*sizeof(double));
104 count=0;
105 for (i=0;i<nbe;i++){
106 if (edges[i].Required()){
107 bamggeom->RequiredEdges[count]=i+1; //back to Matlab indexing
108 count=count+1;
109 }
110 }
111 }
112 else{
113 bamggeom->RequiredEdges=NULL;
114 }
115
116 //No corners
117
118 //RequiredVertices
119 if(verbose>3) printf(" writing RequiredVertices\n");
120 bamggeom->NumRequiredVertices=nbreqv;
121 xfree((void**)&bamggeom->RequiredVertices);
122 if (nbreqv){
123 bamggeom->RequiredVertices=(double*)xmalloc(1*nbreqv*sizeof(double));
124 count=0;
125 for (i=0;i<nbe;i++){
126 if (vertices[i].Required()){
127 bamggeom->RequiredVertices[count]=i+1; //back to Matlab indexing
128 count=count+1;
129 }
130 }
131 }
132 else{
133 bamggeom->RequiredVertices=NULL;
134 }
135
136 //SubDomains
137 if(verbose>3) printf(" writing SubDomains\n");
138 bamggeom->NumSubDomains=NbSubDomains;
139 xfree((void**)&bamggeom->SubDomains);
140 if (NbSubDomains){
141 bamggeom->SubDomains=(double*)xmalloc(4*NbSubDomains*sizeof(double));
142 for (i=0;i<NbSubDomains;i++){
143 bamggeom->SubDomains[4*i+0]=2;
144 bamggeom->SubDomains[4*i+1]=Number(subdomains[i].edge)+1; //back to Matlab indexing
145 bamggeom->SubDomains[4*i+2]=subdomains[i].sens;
146 bamggeom->SubDomains[4*i+3]=subdomains[i].ref;
147 }
148 }
149 else{
150 bamggeom->SubDomains=NULL;
151 }
152
153 //TangentAtEdges
154 if(verbose>3) printf(" writing TangentAtEdges\n");
155 bamggeom->NumTangentAtEdges=nbtan;
156 xfree((void**)&bamggeom->TangentAtEdges);
157 if (nbtan){
158 bamggeom->TangentAtEdges=(double*)xmalloc(4*nbtan*sizeof(double));
159 count=0;
160 for (i=0;i<nbe;i++){
161 if (edges[i].TgA() && edges[i][0].Corner()){
162 bamggeom->TangentAtEdges[4*i+0]=i+1; //back to Matlab indexing
163 bamggeom->TangentAtEdges[4*i+1]=1;
164 bamggeom->TangentAtEdges[4*i+2]=edges[i].tg[0].x;
165 bamggeom->TangentAtEdges[4*i+3]=edges[i].tg[0].y;
166 }
167 if (edges[i].TgB() && edges[i][1].Corner()){
168 bamggeom->TangentAtEdges[4*i+0]=i+1; //back to Matlab indexing
169 bamggeom->TangentAtEdges[4*i+1]=2;
170 bamggeom->TangentAtEdges[4*i+2]=edges[i].tg[1].x;
171 bamggeom->TangentAtEdges[4*i+3]=edges[i].tg[2].y;
172 }
173 count=count+1;
174 }
175 }
176 else{
177 bamggeom->TangentAtEdges=NULL;
178 }
179 }
180 /*}}}1*/
181 /*FUNCTION Geometry::ReadGeometry{{{1*/
182 void Geometry::ReadGeometry(BamgGeom* bamggeom,BamgOpts* bamgopts){
183
184 int verbose;
185 nbiv=nbv=nbvx=0;
186 nbe=nbt=nbtx=0;
187 NbOfCurves=0;
188
189 Real8 Hmin = HUGE_VAL;// the infinie value
190 Int4 hvertices=0;
191 int i,j,k,n;
192
193 //initialize some variables
194 int Version=1,dim=2;
195 nbv=bamggeom->NumVertices;
196 nbe=bamggeom->NumEdges;
197 nbvx = nbv;
198 nbiv = nbv;
199 verbose=bamgopts->verbose;
200
201 //some checks
202 if (bamggeom->NumVertices<=0 || bamggeom->Vertices==NULL){
203 throw ErrorException(__FUNCT__,exprintf("the domain provided does not contain any vertex"));
204 }
205 if (bamggeom->NumEdges<=0 || bamggeom->Edges==NULL){
206 throw ErrorException(__FUNCT__,exprintf("the domain provided does not contain any edge"));
207 }
208
209 //Vertices
210 if (bamggeom->Vertices){
211 if(verbose>3) printf(" processing Vertices\n");
212 vertices = new GeometricalVertex[nbvx];
213 for (i=0;i<nbv;i++) {
214 vertices[i].r.x=(double)bamggeom->Vertices[i*3+0];
215 vertices[i].r.y=(double)bamggeom->Vertices[i*3+1];
216 vertices[i].ReferenceNumber=(Int4)bamggeom->Vertices[i*3+2];
217 vertices[i].DirOfSearch=NoDirOfSearch;
218 vertices[i].color =0;
219 vertices[i].Set();
220 }
221 //find domain extrema for pmin,pmax
222 pmin = vertices[0].r;
223 pmax = vertices[0].r;
224 for (i=0;i<nbv;i++) {
225 pmin.x = Min(pmin.x,vertices[i].r.x);
226 pmin.y = Min(pmin.y,vertices[i].r.y);
227 pmax.x = Max(pmax.x,vertices[i].r.x);
228 pmax.y = Max(pmax.y,vertices[i].r.y);
229 }
230 R2 DD05 = (pmax-pmin)*0.05;
231 pmin -= DD05;
232 pmax += DD05;
233 coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
234 if(coefIcoor <=0){
235 throw ErrorException(__FUNCT__,exprintf("coefIcoor should be positive"));
236 }
237 }
238 else{
239 throw ErrorException(__FUNCT__,exprintf("No Vertex provided"));
240 }
241
242 //Edges
243 if (bamggeom->Edges){
244 int i1,i2;
245 R2 zero2(0,0);
246 Real4 *len =0;
247
248 if(verbose>3) printf(" processing Edges\n");
249 edges = new GeometricalEdge[nbe];
250
251 if (!hvertices) {
252 len = new Real4[nbv];
253 for(i=0;i<nbv;i++)
254 len[i]=0;
255 }
256
257 for (i=0;i<nbe;i++){
258 i1=(int)bamggeom->Edges[i*3+0]-1; //-1 for C indexing
259 i2=(int)bamggeom->Edges[i*3+1]-1; //-1 for C indexing
260 edges[i].ref=(Int4)bamggeom->Edges[i*3+2];
261 edges[i].v[0]= vertices + i1;
262 edges[i].v[1]= vertices + i2;
263 R2 x12 = vertices[i2].r-vertices[i1].r;
264 Real8 l12=sqrt((x12,x12));
265 edges[i].tg[0]=zero2;
266 edges[i].tg[1]=zero2;
267 edges[i].SensAdj[0] = edges[i].SensAdj[1] = -1;
268 edges[i].Adj[0] = edges[i].Adj[1] = 0;
269 edges[i].flag = 0;
270 if (!hvertices) {
271 vertices[i1].color++;
272 vertices[i2].color++;
273 len[i1] += l12;
274 len[i2] += l12;
275 }
276
277 Hmin = Min(Hmin,l12);
278 }
279
280 // definition the default of the given mesh size
281 if (!hvertices){
282 for (i=0;i<nbv;i++)
283 if (vertices[i].color > 0)
284 vertices[i].m= Metric(len[i] /(Real4) vertices[i].color);
285 else
286 vertices[i].m= Metric(Hmin);
287 delete [] len;
288 }
289 }
290 else{
291 throw ErrorException(__FUNCT__,exprintf("No edges provided"));
292 }
293
294 //hVertices
295 if(bamggeom->hVertices){
296 if(verbose>3) printf(" processing hVertices\n");
297 for (i=0;i< nbv;i++){
298 vertices[i].m=Metric((Real4)bamggeom->hVertices[i]);
299 }
300 }
301 else{
302 if(verbose>3) printf(" no hVertices found\n");
303 }
304
305 //MetricVertices
306 if(bamggeom->MetricVertices){
307 if(verbose>3) printf(" processing MetricVertices\n");
308 hvertices=1;
309 for (i=0;i< nbv;i++) {
310 vertices[i].m = Metric((Real4)bamggeom->MetricVertices[i*3+0],(Real4)bamggeom->MetricVertices[i*3+1],(Real4)bamggeom->MetricVertices[i*3+2]);
311 }
312 }
313 else{
314 if(verbose>3) printf(" no MetricVertices found\n");
315 }
316
317 //h1h2VpVertices
318 if(bamggeom->h1h2VpVertices){
319 if(verbose>3) printf(" processing h1h2VpVertices\n");
320 Real4 h1,h2,v1,v2;
321 hvertices =1;
322 for (i=0;i< nbv;i++) {
323 h1=(Real4)bamggeom->MetricVertices[i*4+0];
324 h2=(Real4)bamggeom->MetricVertices[i*4+1];
325 v1=(Real4)bamggeom->MetricVertices[i*4+2];
326 v2=(Real4)bamggeom->MetricVertices[i*4+3];
327 vertices[i].m = Metric(MatVVP2x2(1/(h1*h1),1/(h2*h2),D2(v1,v2)));
328 }
329 }
330 else{
331 if(verbose>3) printf(" no h1h2VpVertices found\n");
332 }
333
334 //MaximalAngleOfCorner
335 if (bamgopts->MaximalAngleOfCorner){
336 if(verbose>3) printf(" processing MaximalAngleOfCorner\n");
337 MaximalAngleOfCorner=bamgopts->MaximalAngleOfCorner*Pi/180;
338 }
339 else{
340 if(verbose>3) printf(" no MaximalAngleOfCorner found\n");
341 }
342
343 //TangentAtEdges
344 if (bamggeom->TangentAtEdges){
345 if(verbose>3) printf(" processing TangentAtEdges");
346 int n,i,j,k;
347 R2 tg;
348
349 n=bamggeom->NumTangentAtEdges;
350 for (k=0;k<n;k++) {
351 i=(int)bamggeom->TangentAtEdges[k*4+0]-1; //for C indexing
352 j=(int)bamggeom->TangentAtEdges[k*4+1]-1; //for C indexing
353 tg.x=bamggeom->TangentAtEdges[k*4+2];
354 tg.y=bamggeom->TangentAtEdges[k*4+3];
355 if (j!=0 && j!=1){
356 throw ErrorException(__FUNCT__,exprintf("TangentAtEdges second index should be 1 or 2 only"));
357 }
358 edges[i].tg[j] = tg;
359 }
360 }
361 else{
362 if(verbose>3) printf(" no TangentAtEdges found\n");
363 }
364
365 //Corners
366 if(bamggeom->Corners){
367 if(verbose>3) printf(" processing Corners");
368 n=bamggeom->NumCorners;
369 for (i=0;i<n;i++) {
370 j=(int)bamggeom->Corners[i]-1; //for C indexing
371 if (j>nbv-1 || j<0){
372 throw ErrorException(__FUNCT__,exprintf("Bad corner definition: should in [0 %i]",nbv));
373 }
374 vertices[j].SetCorner();
375 vertices[j].SetRequired(); }
376 }
377 else{
378 if(verbose>3) printf(" no Corners found\n");
379 }
380
381 //RequiredVertices
382 if(bamggeom->RequiredVertices){
383 if(verbose>3) printf(" processing RequiredVertices");
384 n=bamggeom->NumRequiredVertices;
385 for (i=0;i<n;i++) {
386 j=(int)bamggeom->RequiredVertices[i]-1; //for C indexing
387 if (j>nbv-1 || j<0){
388 throw ErrorException(__FUNCT__,exprintf("Bad RequiredVerticess definition: should in [0 %i]",nbv));
389 }
390 vertices[j].SetRequired(); }
391 }
392 else{
393 if(verbose>3) printf(" no RequiredVertices found\n");
394 }
395
396 //RequiredEdges
397 if(bamggeom->RequiredEdges){
398 if(verbose>3) printf(" processing RequiredEdges");
399 n=bamggeom->NumRequiredEdges;
400 for (i=0;i<n;i++) {
401 j=(int)bamggeom->RequiredEdges[i]-1; //for C indexing
402 if (j>nbe-1 || j<0){
403 throw ErrorException(__FUNCT__,exprintf("Bad RequiredEdges definition: should in [0 %i]",nbv));
404 }
405 edges[j].SetRequired(); }
406 }
407 else{
408 if(verbose>3) printf(" no RequiredEdges found\n");
409 }
410
411 //SubDomain
412 if(bamggeom->SubDomains){
413 Int4 i0,i1,i2,i3;
414 if(verbose>3) printf(" processing SubDomains\n");
415 NbSubDomains=bamggeom->NumSubDomains;
416 subdomains = new GeometricalSubDomain[NbSubDomains];
417 for (i=0;i<NbSubDomains;i++) {
418 i0=(int)bamggeom->SubDomains[i*4+0];
419 i1=(int)bamggeom->SubDomains[i*4+1];
420 i2=(int)bamggeom->SubDomains[i*4+2];
421 i3=(int)bamggeom->SubDomains[i*4+3];
422 if (i0!=2) throw ErrorException(__FUNCT__,exprintf("Bad Subdomain definition: first number should be 2 (for Edges)"));
423 if (i1>nbe || i1<=0) throw ErrorException(__FUNCT__,exprintf("Bad Subdomain definition: second number should in [1 %i] (edge number)",nbe));
424 subdomains[i].edge=edges + (i1-1);
425 subdomains[i].sens = (int) i2;
426 subdomains[i].ref = i3;
427 }
428 }
429 else{
430 if(verbose>3) printf(" no SubDomains found\n");
431 }
432 }
433 /*}}}1*/
434
435 /*Other*/
436
437void Geometry::EmptyGeometry() // empty geometry
438 {
439 OnDisk=0;
440 NbRef=0;
441 name =0;
442 quadtree=0;
443 curves=0;
444 // edgescomponante=0;
445 triangles=0;
446 edges=0;
447 vertices=0;
448 NbSubDomains=0;
449 // nbtf=0;
450// BeginOfCurve=0;
451 nbiv=nbv=nbvx=0;
452 nbe=nbt=nbtx=0;
453 NbOfCurves=0;
454// BeginOfCurve=0;
455 subdomains=0;
456 MaximalAngleOfCorner = 10*Pi/180;
457 }
458
459
460
461Geometry::Geometry(const Geometry & Gh)
462 { Int4 i;
463 *this = Gh;
464 NbRef =0;
465 quadtree=0;
466 name = new char[strlen(Gh.name)+4];
467 strcpy(name,"cp:");
468 strcat(name,Gh.name);
469 vertices = nbv ? new GeometricalVertex[nbv] : NULL;
470 triangles = nbt ? new Triangle[nbt]:NULL;
471 edges = nbe ? new GeometricalEdge[nbe]:NULL;
472 curves= NbOfCurves ? new Curve[NbOfCurves]:NULL;
473 subdomains = NbSubDomains ? new GeometricalSubDomain[NbSubDomains]:NULL;
474 for (i=0;i<nbv;i++)
475 vertices[i].Set(Gh.vertices[i],Gh,*this);
476 for (i=0;i<nbe;i++)
477 edges[i].Set(Gh.edges[i],Gh,*this);
478 for (i=0;i<NbOfCurves;i++)
479 curves[i].Set(Gh.curves[i],Gh,*this);
480 for (i=0;i<NbSubDomains;i++)
481 subdomains[i].Set(Gh.subdomains[i],Gh,*this);
482
483 // for (i=0;i<nbt;i++)
484 // triangles[i].Set(Gh.triangles[i],Gh,*this);
485 assert(!nbt);
486 }
487
488
489GeometricalEdge* Geometry::Contening(const R2 P, GeometricalEdge * start) const
490{
491 GeometricalEdge* on =start,* pon=0;
492 // walk with the cos on geometry
493 // cout << P ;
494 int k=0;
495 while(pon != on)
496 {
497 pon = on;
498 assert(k++<100);
499 R2 A= (*on)[0];
500 R2 B= (*on)[1];
501 R2 AB = B-A;
502 R2 AP = P-A;
503 R2 BP = P-B;
504 // cout << ":: " << on - edges << " " << AB*AP << " " << AB*BP << " " << A << B << endl;
505 if ( (AB,AP) < 0)
506 on = on->Adj[0];
507 else if ( (AB,BP) > 0)
508 on = on->Adj[1];
509 else
510 return on;
511 }
512 return on;
513}
514GeometricalEdge* Geometry::ProjectOnCurve(const Edge & e,Real8 s,Vertex &V,VertexOnGeom &GV ) const
515 {
516 Real8 save_s=s;
517 int NbTry=0;
518retry:
519 s=save_s;
520 GeometricalEdge * on = e.on;
521 assert(on);
522 assert( e[0].on && e[1].on);
523 const Vertex &v0=e[0],&v1=e[1];
524 V.m = Metric(1.0-s, v0,s, v1);
525#define MXE__LINE __LINE__+1
526 const int mxe =100;
527 GeometricalEdge *ge[mxe+1];
528 int sensge[mxe+1];
529 Real8 lge[mxe+1];
530 int bge=mxe/2,tge=bge;
531 ge[bge] = e.on;
532 sensge[bge]=1;
533
534 R2 V0 = v0,V1=v1,V01=V1-V0;
535 VertexOnGeom vg0= *v0.on, vg1=*v1.on;
536 if(NbTry) cout << "bug: s==== " << s << " e=" << V0 << " " << V1 << endl;
537
538 // GeometricalEdge * eg0 = e.on,* eg1 = e.on, *eg=NULL;
539 GeometricalEdge * eg0=on, *eg1=on;
540 R2 Ag=(R2) (*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag;
541 if(NbTry) cout <<" G edge= " << Ag << Bg << endl << " v edge" << V01 << " v geom " << AB << (V01,AB) <<endl;
542 int OppositeSens = (V01,AB) < 0;
543 int sens0=0,sens1=1;
544 if (OppositeSens)
545 s=1-s,Exchange(vg0,vg1),Exchange(V0,V1);
546 if(NbTry) cout << "bug: edge = " << v0.r << " -> " << v1.r << endl
547 << "sg 0 = " << vg0
548 << " on = " << Number(on) << ":" << Ag << Bg << "; "
549 << " sg 1= " << vg1
550 << "--------------------------------------------" << endl;
551 while (eg0 != (GeometricalEdge*) vg0 && (*eg0)(sens0) != (GeometricalVertex*) vg0)
552 {
553 if (bge<=0) {
554 // int kkk;
555 // if (NbTry) cout <<"Read (int) to Show Sioux window", cin >> kkk ;
556 if(NbTry)
557 {
558 cerr << " -- Fatal Error: on the class triangles before call Geometry::ProjectOnCurve" << endl;
559 cerr << " The mesh of the Geometry is to fine: ";
560 cerr << " 1) a mesh edge contening more than "<< mxe/2 << " geometrical edges." << endl;
561 cerr << " 2) code bug : be sure that we call Triangles::SetVertexFieldOn() before " << endl;
562 cerr << " To solve the problem do a coarsening of the geometrical mesh " << endl;
563 cerr << " or change the constant value of mxe in " << __FILE__ << " line " << MXE__LINE << "( dangerous way )" << endl;
564 MeshError(222);
565 }
566 NbTry++;
567 goto retry;}
568 GeometricalEdge* tmpge = eg0;
569 if(NbTry)
570 cout << "bug: --Edge @" << Number(tmpge) << " = "<< Number(eg0) << ":" <<Number(eg0->Adj[0]) << "," <<
571 Number(eg0->Adj[1]) <<"," ;
572 ge[--bge] =eg0 = eg0->Adj[sens0];
573 assert(bge>=0 && bge <= mxe);
574 sens0 = 1-( sensge[bge] = tmpge->SensAdj[sens0]);
575 if(NbTry)
576 cout << "bug: Edge " << Number(eg0) << " "<< 1-sens0 << " S "
577 << Number((*eg0)[1-sens0]) <<":" << Number(eg0->Adj[0]) << ","
578 << Number(eg0->Adj[1]) <<"," << endl
579 <<Number(eg0)<< (*eg0)[sens0].r << "v = " << Number((*eg1)(sens0)) << " e = " << eg0 << endl;
580 }
581 if(NbTry) cout << Number((GeometricalEdge*) vg1) << " " << Number((GeometricalVertex*) vg1) << endl;
582 while (eg1 != (GeometricalEdge*) vg1 && (*eg1)(sens1) != (GeometricalVertex*) vg1)
583 {
584 if(tge>=mxe ) {
585 cerr << " --Fatal Error: on the class triangles before call Geometry::ProjectOnCurve" << endl;
586 NbTry++;
587 if (NbTry<2) goto retry;
588 cerr << " The mesh of the Geometry is to fine:" ;
589 cerr << " 1) a mesh edge contening more than "<< mxe/2 << " geometrical edges." << endl;
590 cerr << " 2) code bug : be sure that we call Triangles::SetVertexFieldOn() before " << endl;
591 cerr << " To solve the problem do a coarsening of the geometrical mesh " << endl;
592 cerr << " or change the constant value of mxe in " << __FILE__ << " line " << MXE__LINE << "( dangerous way )" << endl;
593 MeshError(223);
594 }
595
596 GeometricalEdge* tmpge = eg1;
597 if(NbTry)
598 cout << "++Edge @" << tmpge << " = " << Number(eg1) <<"%" << Number(eg1->Adj[0]) << ","
599 << Number(eg1->Adj[1]) <<"," ;
600 ge[++tge] =eg1 = eg1->Adj[sens1];
601 sensge[tge]= sens1 = 1-tmpge->SensAdj[sens1];
602 assert(tge>=0 && tge <= mxe);
603 if(NbTry)
604 cout << " Edge " << Number(eg1) << " " << sens1 << " S "
605 <<Number((*eg1)[sens1]) <<"%"<< Number(eg1->Adj[0]) << "," << Number(eg1->Adj[1]) <<","
606 <<Number(eg1)<< (*eg1)[sens1].r << "v = " << Number((*eg1)(sens1)) << " e = " << Number(eg1) << endl;
607 }
608
609 if(NbTry) cout << endl;
610
611
612 if ( (*eg0)(sens0) == (GeometricalVertex*) vg0 )
613 vg0 = VertexOnGeom( *(Vertex *) vg0,*eg0,sens0);
614
615 if ( (*eg1)(sens1) == (GeometricalVertex*) vg1)
616 vg1 = VertexOnGeom( *(Vertex *) vg1,*eg1,sens1);
617
618 Real8 sg;
619 // cout << " " << Number(on) << " " << Number(eg0) << " " << Number(eg1) << " " ;
620 if (eg0 == eg1) {
621 register Real8 s0= vg0,s1=vg1;
622 sg = s0 * (1.0-s) + s * s1;
623 // cout <<" s0=" << s0 << " s1=" << s1
624 // << " s = " << s << " sens= " << OppositeSens << "\t\t sg = " << sg << endl ;
625 on=eg0;}
626 else {
627 R2 AA=V0,BB;
628 Real8 s0,s1;
629
630 //cout << endl << "s= " << s << Number(eg0) << " " << (Real8) vg0 << " "
631 // << Number(eg1) << " " << (Real8) vg1 << V0 << V1 << " Interpol = "
632 // << V0*(1-s)+V1*s << ";;; " << endl;
633 int i=bge;
634 Real8 ll=0;
635 for(i=bge;i<tge;i++)
636 {
637 assert( i>=0 && i <= mxe);
638 BB = (*ge[i])[sensge[i]];
639 lge[i]=ll += Norme2(AA-BB);
640 // cout << " ll " << i << BB << ll << " " <<sensge[i] <<" on = " <<
641 // Number(ge[i]) << " sens= " << sensge[i] ;
642 AA=BB ;}
643 lge[tge]=ll+=Norme2(AA-V1);
644 // cout << " ll " << tge << " " << ll << sensge[tge]
645 // <<" on = " << Number(ge[tge]) << " sens= " << sensge[tge] << endl;
646 // search the geometrical edge
647 assert(s <= 1.0);
648 Real8 ls= s*ll;
649 on =0;
650 s0 = vg0;
651 s1= sensge[bge];
652 Real8 l0=0,l1;
653 i=bge;
654 while ( (l1=lge[i]) < ls ) {
655 assert(i >= 0 && i <= mxe);
656 i++,s0=1-(s1=sensge[i]),l0=l1;}
657 on=ge[i];
658 if (i==tge)
659 s1=vg1;
660
661 s=(ls-l0)/(l1-l0);
662 // cout << "on =" << Number(on) << sens0 << sens1 << "s0 " << s0 << " s1 ="
663 // << s1 << " l0 =" << l0 << " ls= " << ls << " l1= " << l1 << " s= " << s;
664 sg = s0 * (1.0-s) + s * s1;
665 }
666 assert(on);
667 // assert(sg && sg-1);
668 V.r= on->F(sg);
669 // if (eg0 != eg1)
670 // cout << "----- sg = "<< sg << " Sens =" << OppositeSens << " Edge = "
671 // << Number(on) <<" V=" << V << endl;
672 GV=VertexOnGeom(V,*on,sg);
673 return on;
674 }
675
676void Geometry::AfterRead()
677 {// -----------------
678 long int verbosity=0;
679
680 if (verbosity>20)
681 cout << "Geometry::AfterRead()" << nbv << " " << nbe << endl;
682 Int4 i,k=0; ;
683 int jj; // jj in [0,1]
684 Int4 * hv = new Int4 [ nbv];
685 Int4 * ev = new Int4 [ 2 * nbe ];
686 float * eangle = new float[ nbe ];
687 {
688 double eps = 1e-20;
689 QuadTree quadtree; // to find same vertices
690 Vertex * v0 = vertices;
691 GeometricalVertex * v0g = (GeometricalVertex *) (void *) v0;
692 int k=0;
693 for (i=0;i<nbv;i++)
694 vertices[i].link = vertices +i;
695 for (i=0;i<nbv;i++)
696 {
697 vertices[i].i = toI2(vertices[i].r); // set integer coordinate
698 Vertex *v= quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y);
699 if( v && Norme1(v->r - vertices[i]) < eps )
700 { // link v & vertices[i]
701 // vieille ruse pour recuperer j
702 GeometricalVertex * vg = (GeometricalVertex *) (void *) v;
703 int j = vg-v0g;
704 assert( v == & (Vertex &) vertices[j]);
705 vertices[i].link = vertices + j;
706 k++;
707 }
708 else quadtree.Add(vertices[i]);
709 }
710 if (k) {
711 cout << " Number of distinte vertices " << nbv - k << " Over " << nbv << endl;
712 //if (verbosity>10)
713 {
714 cout << " The duplicate vertex " << endl;
715 for (i=0;i<nbv;i++)
716 if (!vertices[i].IsThe())
717 cout << " " << i << " and " << Number(vertices[i].The()) << endl;
718 MeshError(102);
719 //throw(ErrorExec("exit",1));
720 }
721 }
722
723 // verification of cracked edge
724 int err =0;
725 for (i=0;i<nbe;i++)
726 if (edges[i].Cracked() )
727 {
728 // verification of crack
729 GeometricalEdge & e1=edges[i];
730 GeometricalEdge & e2=*e1.link;
731 cerr << i << " " << e1[0].The() << " " << e2[0].The() << " " << e1[1].The() << " " << e2[1].The() << endl;
732 if ( e1[0].The() == e2[0].The() && e1[1].The() == e2[1].The() )
733 {
734 }
735 else
736 if ( e1[0].The() == e2[1].The() && e1[1].The() == e2[0].The() )
737 {
738 }
739 else
740 {
741 err++;
742 cerr << " Cracked edges with no same vertex " << &e1-edges << " " << &e2 -edges << endl;
743 }
744 }
745 else
746 {
747 // if (!edges[i][0].IsThe()) err++;
748 // if (!edges[i][1].IsThe()) err++;
749 }
750 if (err)
751 {
752 cerr << " Some vertex was not distint and not on cracked edge " << err<< endl;
753 MeshError(222);
754 }
755 }
756 if(verbosity>7)
757 for (i=0;i<nbv;i++)
758 if (vertices[i].Required())
759 cout << " The geo vertices " << i << " is required" << endl;
760
761 for (i=0;i<nbv;i++)
762 hv[i]=-1;// empty list
763
764 for (i=0;i<nbe;i++)
765 {
766 R2 v10 = edges[i].v[1]->r - edges[i].v[0]->r;
767 Real8 lv10 = Norme2(v10);
768 if(lv10 == 0) {
769 cerr << "The length of " <<i<< "th Egde is 0 " << endl ;
770 MeshError(1);}
771 eangle[i] = atan2(v10.y,v10.x) ; // angle in [ -Pi,Pi ]
772 if(verbosity>9)
773 cout << " angle edge " << i <<" " << eangle[i]*180/Pi<< v10<<endl;
774 for (jj=0;jj<2;jj++)
775 { // generation of list
776 Int4 v = Number(edges[i].v[jj]);
777 ev[k] = hv[v];
778 hv[v] = k++;
779 }
780 }
781 // bulle sort on the angle of edge
782 for (i=0;i<nbv;i++) {
783 int exch = 1,ord =0;
784 while (exch) {
785 exch = 0;
786 Int4 *p = hv + i, *po = p;
787 Int4 n = *p;
788 register float angleold = -1000 ; // angle = - infini
789 ord = 0;
790 while (n >=0)
791 {
792 ord++;
793 register Int4 i1= n /2;
794 register Int4 j1 = n % 2;
795 register Int4 *pn = ev + n;
796 float angle = j1 ? OppositeAngle(eangle[i1]): eangle[i1];
797 n = *pn;
798 if (angleold > angle) // exch to have : po -> pn -> p
799 exch=1,*pn = *po,*po=*p,*p=n,po = pn;
800 else // to have : po -> p -> pn
801 angleold = angle, po = p,p = pn;
802 }
803 } // end while (exch)
804
805 if (ord >= 1 )
806 { /*
807 Int4 n = hv[i];
808 while ( n >=0)
809 { Int4 i1 = n/2,j1 = n%2;
810 //float a = 180*(j1 ? OppositeAngle(eangle[i1]): eangle[i1])/Pi;
811 n = ev[n];
812 }
813 */
814 }
815 if(ord == 2) { // angulare test to find a corner
816 Int4 n1 = hv[i];
817 Int4 n2 = ev[n1];
818 Int4 i1 = n1 /2, i2 = n2/2; // edge number
819 Int4 j1 = n1 %2, j2 = n2%2; // vertex in the edge
820 float angle1= j1 ? OppositeAngle(eangle[i1]) : eangle[i1];
821 float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2];
822 float da12 = Abs(angle2-angle1);
823 if(verbosity>9)
824 cout <<" check angle " << i << " " << i1 << " " << i2 << " " << 180*(da12)/Pi
825 << " " << 180*MaximalAngleOfCorner/Pi << vertices[i] << endl;
826
827 if (( da12 >= MaximalAngleOfCorner )
828 && (da12 <= 2*Pi -MaximalAngleOfCorner)) {
829 vertices[i].SetCorner() ;
830 if(verbosity>7)
831 cout << " The vertex " << i << " is a corner (angle) "
832 << 180*(da12)/ Pi<< " " << 180*MaximalAngleOfCorner/Pi << endl;}
833 // if the ref a changing then is SetRequired();
834
835 if (edges[i1].flag != edges[i2].flag || edges[i1].Required())
836 {
837 vertices[i].SetRequired();
838 if(verbosity>7)
839 cout << " The vertex " << i << " is Required the flag change (crack or equi, or require)" << endl;}
840
841 if (edges[i1].ref != edges[i2].ref) {
842 vertices[i].SetRequired();
843 if(verbosity>7)
844 cout << " The vertex " << i << " is Required ref" << endl;}
845 } ;
846
847 if(ord != 2) {
848 vertices[i].SetCorner();
849 if(verbosity>7)
850 cout << " the vertex " << i << " is a corner ordre = " << ord << endl;
851 }
852 // close the liste around the vertex
853 { Int4 no=-1, ne = hv[i];
854 while ( ne >=0)
855 ne = ev[no=ne];
856 if(no>=0)
857 ev[no] = hv[i];
858 } // now the list around the vertex is circular
859
860 } // end for (i=0;i<nbv;i++)
861
862 k =0;
863 for (i=0;i<nbe;i++)
864 for (jj=0;jj<2;jj++){
865 Int4 n1 = ev[k++];
866 Int4 i1 = n1/2 ,j1=n1%2;
867 if( edges[i1].v[j1] != edges[i].v[jj])
868 { cerr << " Bug Adj edge " << i << " " << jj <<
869 " et " << i1 << " " << j1 << " k=" << k;
870 cerr << Number(edges[i].v[jj]) <<" <> "
871 << Number(edges[i1].v[j1]) <<endl;
872 cerr << "edge " << Number(edges[i].v[0]) << " "
873 << Number(edges[i].v[1]) << endl;
874 // cerr << "in file " <<filename <<endl;
875 MeshError(1);
876 }
877 edges[i1].Adj[j1] = edges + i;
878 edges[i1].SensAdj[j1] = jj;
879 if (verbosity>10)
880 cout << " edges. Adj " << i1 << " " << j1 << " <--- " << i << " " << jj << endl;
881 }
882
883 // generation of all the tangente
884 for (i=0;i<nbe;i++) {
885 R2 AB = edges[i].v[1]->r -edges[i].v[0]->r;
886 Real8 lAB = Norme2(AB); // length of current edge AB
887 Real8 ltg2[2];
888 ltg2[0]=0;ltg2[1]=0;
889 for (jj=0;jj<2;jj++) {
890 R2 tg = edges[i].tg[jj];
891 Real8 ltg = Norme2(tg); // length of tg
892 if(ltg == 0) {// no tg
893 if( ! edges[i].v[jj]->Corner()) { // not a Corner
894 tg = edges[i].v[1-jj]->r
895 - edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r;
896 ltg = Norme2(tg);
897 tg = tg *(lAB/ltg),ltg=lAB;
898 /*
899 if(edges[i].ref >=4)
900 cout << " tg " << tg.x << " "<< tg.y << " " << edges[i].v[1-jj]->r << edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r << " y-y = "
901 << edges[i].v[1-jj]->r.y -edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r.y << endl;
902 */
903 }
904
905 //else ;// a Corner with no tangent => nothing to do
906 } // a tg
907 else
908 tg = tg *(lAB/ltg),ltg=lAB;
909 ltg2[jj] = ltg;
910 if ( (tg,AB) < 0)
911 tg = -tg;
912 //if(edges[i].ref >=4) cout << " tg = " << tg << endl;
913 edges[i].tg[jj] = tg;
914 } // for (jj=0;jj<2;jj++)
915
916 if (ltg2[0]!=0) edges[i].SetTgA();
917 if (ltg2[1]!=0) edges[i].SetTgB();
918 } // for (i=0;i<nbe;i++)
919
920 if(verbosity>7)
921 for (i=0;i<nbv;i++)
922 if (vertices[i].Required())
923 cout << " The geo vertices " << i << " is required " << endl;
924
925 for (int step=0;step<2;step++)
926 {
927 for (i=0;i<nbe;i++)
928 edges[i].SetUnMark();
929
930 NbOfCurves = 0;
931 Int4 nbgem=0;
932 for (int level=0;level < 2 && nbgem != nbe;level++)
933 for (i=0;i<nbe;i++) {
934 GeometricalEdge & ei = edges[i];
935 for(jj=0;jj<2;jj++)
936 if (!ei.Mark() && (level || ei[jj].Required())) {
937 // warning ei.Mark() can be change in loop for(jj=0;jj<2;jj++)
938 int k0=jj,k1;
939 GeometricalEdge *e = & ei;
940 GeometricalVertex *a=(*e)(k0); // begin
941 if(curves) {
942 curves[NbOfCurves].be=e;
943 curves[NbOfCurves].kb=k0;
944 }
945 int nee=0;
946 for(;;) {
947 nee++;
948 k1 = 1-k0; // next vertex of the edge
949 e->SetMark();
950 nbgem++;
951 e->CurveNumber=NbOfCurves;
952 if(curves) {
953 curves[NbOfCurves].ee=e;
954 curves[NbOfCurves].ke=k1;
955 }
956
957 GeometricalVertex *b=(*e)(k1);
958 if (a == b || b->Required() ) break;
959 k0 = e->SensAdj[k1];// vertex in next edge
960 e = e->Adj[k1]; // next edge
961
962 }// for(;;)
963 if(verbosity>10 && curves==0) cout << NbOfCurves <<" curve : nb edges= "<< nee<< endl;
964 NbOfCurves++;
965 if(level) {
966 if(verbosity>4)
967 cout << " Warning: Curve "<< NbOfCurves << " without required vertex "
968 << "so the vertex " << Number(a) << " become required " <<endl;
969 a->SetRequired();
970 }
971
972 }}
973 assert(nbgem && nbe);
974
975 if(step==0) {
976 curves = new Curve[NbOfCurves];
977 }
978 }
979 for(int i=0;i<NbOfCurves ;i++)
980 {
981 GeometricalEdge * be=curves[i].be, *eqbe=be->link;
982 //GeometricalEdge * ee=curves[i].ee, *eqee=be->link;
983 curves[i].master=true;
984 if(be->Equi() || be->ReverseEqui() )
985 {
986 assert(eqbe);
987 int nc = eqbe->CurveNumber;
988 assert(i!=nc);
989 curves[i].next=curves[nc].next;
990 curves[i].master=false;
991 curves[nc].next=curves+i;
992 if(be->ReverseEqui())
993 curves[i].Reverse();
994 }
995 }
996
997 if(verbosity>3)
998 cout << " End ReadGeometry: Number of curves in geometry is " << NbOfCurves <<endl;
999 if(verbosity>4)
1000 for(int i=0;i<NbOfCurves ;i++)
1001 {
1002 cout << " Curve " << i << " begin e=" << Number(curves[i].be) << " k=" << curves[i].kb
1003 << " end e= " << Number(curves[i].ee) << " k=" << curves[i].ke << endl;
1004 }
1005 delete []ev;
1006 delete []hv;
1007 delete []eangle;
1008
1009}
1010Geometry::~Geometry()
1011{
1012 long int verbosity=0;
1013
1014 assert(NbRef<=0);
1015 if(verbosity>9)
1016 cout << "DELETE ~Geometry "<< this << endl;
1017 if(vertices) delete [] vertices;vertices=0;
1018 if(edges) delete [] edges;edges=0;
1019 // if(edgescomponante) delete [] edgescomponante; edgescomponante=0;
1020 if(triangles) delete [] triangles;triangles=0;
1021 if(quadtree) delete quadtree;quadtree=0;
1022 if(curves) delete []curves;curves=0;NbOfCurves=0;
1023 if(name) delete [] name;name=0;
1024 if(subdomains) delete [] subdomains;subdomains=0;
1025// if(ordre) delete [] ordre;
1026 EmptyGeometry();
1027
1028}
1029
1030
1031}
Note: See TracBrowser for help on using the repository browser.