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 |
|
---|
10 | namespace 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 |
|
---|
828 | retry:
|
---|
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 | }
|
---|