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

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

created 2 files Triangles and Geometry for corresponding classes
File size: 12.1 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}
Note: See TracBrowser for help on using the repository browser.