Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
makeEllipsoidNektar.c
Go to the documentation of this file.
1 #include <math.h>
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <string.h>
5 #include <time.h>
6 #include <float.h>
7 #define PI 3.14159265358979323846
8 
9 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
10 //Dr. Laurent Risser (Imperial College London)
11 //l.risser@imperial.ac.uk
12 //25.08.2009
13 //
14 //Command line to compile: gcc -o makeEllipsoidNektar makeEllipsoidNektar.c -O -lm
15 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
16 
17 
18 // ---------------------------------------------------------------------------------------
19 // structures definition
20 // ---------------------------------------------------------------------------------------
21 
22 typedef struct
23 {
24 double X;
25 double Y; // coordinates of a vertex
26 double Z;
27 } Vertex;
28 
29 typedef struct
30 {
31 int V1; //id of the first vertex
32 int V2; //id of the second vertex
33 Vertex * SubiVertex;
34 int NbSubdiv;
35 } Edge;
36 
37 typedef struct
38 {
39 int E1; //id of the first edge
40 int E2; //id of the second edge
41 int E3; //id of the third edge
42 } Element;
43 
44 typedef struct
45 {
46 int NbVertexes;
47 Vertex * Vertexes;
48 int NbEdges;
49 Edge * Edges;
50 int NbElements;
51 Element * Elements;
52 } Mesh;
53 
54 
55 
56 
57 // ---------------------------------------------------------------------------------------
58 // functions
59 // ---------------------------------------------------------------------------------------
60 
61 
62 //create a simple cubic mesh
63 Mesh CreateBasicCubicMesh(){
64  Mesh LocMesh;
65  int i;
66 
67 
68  //allocations...
69  LocMesh.NbVertexes=8;
70  LocMesh.Vertexes=(Vertex *)malloc(LocMesh.NbVertexes*sizeof(Vertex));
71 
72  LocMesh.NbEdges=18;
73  LocMesh.Edges=(Edge *)malloc(LocMesh.NbEdges*sizeof(Edge));
74 
75  LocMesh.NbElements=12;
76  LocMesh.Elements=(Element *)malloc(LocMesh.NbElements*sizeof(Element));
77 
78  //fill the coordinates of each point
79  LocMesh.Vertexes[0].X=0; LocMesh.Vertexes[0].Y=0; LocMesh.Vertexes[0].Z=0;
80  LocMesh.Vertexes[1].X=1; LocMesh.Vertexes[1].Y=0; LocMesh.Vertexes[1].Z=0;
81  LocMesh.Vertexes[2].X=1; LocMesh.Vertexes[2].Y=1; LocMesh.Vertexes[2].Z=0;
82  LocMesh.Vertexes[3].X=0; LocMesh.Vertexes[3].Y=1; LocMesh.Vertexes[3].Z=0;
83  LocMesh.Vertexes[4].X=0; LocMesh.Vertexes[4].Y=0; LocMesh.Vertexes[4].Z=1;
84  LocMesh.Vertexes[5].X=1; LocMesh.Vertexes[5].Y=0; LocMesh.Vertexes[5].Z=1;
85  LocMesh.Vertexes[6].X=1; LocMesh.Vertexes[6].Y=1; LocMesh.Vertexes[6].Z=1;
86  LocMesh.Vertexes[7].X=0; LocMesh.Vertexes[7].Y=1; LocMesh.Vertexes[7].Z=1;
87 
88  //fill each edge
89  LocMesh.Edges[0].V1=0; LocMesh.Edges[0].V2=1;
90  LocMesh.Edges[1].V1=1; LocMesh.Edges[1].V2=2;
91  LocMesh.Edges[2].V1=2; LocMesh.Edges[2].V2=3;
92  LocMesh.Edges[3].V1=3; LocMesh.Edges[3].V2=0;
93  LocMesh.Edges[4].V1=4; LocMesh.Edges[4].V2=5;
94  LocMesh.Edges[5].V1=5; LocMesh.Edges[5].V2=6;
95  LocMesh.Edges[6].V1=6; LocMesh.Edges[6].V2=7;
96  LocMesh.Edges[7].V1=7; LocMesh.Edges[7].V2=4;
97  LocMesh.Edges[8].V1=0; LocMesh.Edges[8].V2=4;
98  LocMesh.Edges[9].V1=1; LocMesh.Edges[9].V2=5;
99  LocMesh.Edges[10].V1=2; LocMesh.Edges[10].V2=6;
100  LocMesh.Edges[11].V1=3; LocMesh.Edges[11].V2=7;
101  LocMesh.Edges[12].V1=1; LocMesh.Edges[12].V2=4;
102  LocMesh.Edges[13].V1=1; LocMesh.Edges[13].V2=3;
103  LocMesh.Edges[14].V1=1; LocMesh.Edges[14].V2=6;
104  LocMesh.Edges[15].V1=7; LocMesh.Edges[15].V2=0;
105  LocMesh.Edges[16].V1=7; LocMesh.Edges[16].V2=2;
106  LocMesh.Edges[17].V1=7; LocMesh.Edges[17].V2=5;
107 
108  for (i=0;i<LocMesh.NbEdges;i++) LocMesh.Edges[i].NbSubdiv=1;
109 
110  //fill each Element
111  LocMesh.Elements[0].E1=0; LocMesh.Elements[0].E2=8; LocMesh.Elements[0].E3=12;
112  LocMesh.Elements[1].E1=9; LocMesh.Elements[1].E2=4; LocMesh.Elements[1].E3=12;
113  LocMesh.Elements[2].E1=0; LocMesh.Elements[2].E2=3; LocMesh.Elements[2].E3=13;
114  LocMesh.Elements[3].E1=1; LocMesh.Elements[3].E2=2; LocMesh.Elements[3].E3=13;
115  LocMesh.Elements[4].E1=1; LocMesh.Elements[4].E2=10; LocMesh.Elements[4].E3=14;
116  LocMesh.Elements[5].E1=9; LocMesh.Elements[5].E2=5; LocMesh.Elements[5].E3=14;
117  LocMesh.Elements[6].E1=3; LocMesh.Elements[6].E2=11; LocMesh.Elements[6].E3=15;
118  LocMesh.Elements[7].E1=8; LocMesh.Elements[7].E2=7; LocMesh.Elements[7].E3=15;
119  LocMesh.Elements[8].E1=2; LocMesh.Elements[8].E2=11; LocMesh.Elements[8].E3=16;
120  LocMesh.Elements[9].E1=10; LocMesh.Elements[9].E2=6; LocMesh.Elements[9].E3=16;
121  LocMesh.Elements[10].E1=4; LocMesh.Elements[10].E2=7; LocMesh.Elements[10].E3=17;
122  LocMesh.Elements[11].E1=5; LocMesh.Elements[11].E2=6; LocMesh.Elements[11].E3=17;
123 
124  return LocMesh;
125 }
126 
127 
128 //subdivide into 2 edges the edge 'IdEdge' in the mesh 'LocMesh'
129 void CutEdge(Mesh * LocMesh,int IdEdge){
130  int i;
131  int IdNgbhVortex1;
132  int IdNgbhVortex2;
133  int IdNgbhVortex3;
134  int IdNgbhVortex4;
135  int IdNewVortex;
136  int IdSubdiviedEdge;
137  int IdNewEdge1;
138  int IdNewEdge2;
139  int IdNewEdge3;
140  int IdNgbhEdge1;
141  int IdNgbhEdge2;
142  int IdNgbhEdge3;
143  int IdNgbhEdge4;
144  int IdSubdiviedElement1;
145  int IdSubdiviedElement2;
146  int IdNewElement2;
147  int IdNewElement1;
148 
149  //reallocations in LocMesh
150  LocMesh->NbVertexes++;
151  LocMesh->NbEdges=LocMesh->NbEdges+3;
152  LocMesh->NbElements=LocMesh->NbElements+2;
153 
154  LocMesh->Vertexes=(Vertex *)realloc(LocMesh->Vertexes,LocMesh->NbVertexes*sizeof(Vertex));
155  LocMesh->Edges=(Edge *)realloc(LocMesh->Edges,LocMesh->NbEdges*sizeof(Edge));
156  LocMesh->Elements=(Element *)realloc(LocMesh->Elements,LocMesh->NbElements*sizeof(Element));
157 
158  //identifiers to take into account
159  IdNgbhVortex1=LocMesh->Edges[IdEdge].V1;
160  IdNgbhVortex2=LocMesh->Edges[IdEdge].V2;
161  IdNewVortex=LocMesh->NbVertexes-1;
162 
163  IdSubdiviedEdge=IdEdge;
164  IdNewEdge1=LocMesh->NbEdges-3;
165  IdNewEdge2=LocMesh->NbEdges-2;
166  IdNewEdge3=LocMesh->NbEdges-1;
167 
168  IdSubdiviedElement1=-1;
169  IdSubdiviedElement2=-1;
170  for (i=0;i<LocMesh->NbElements-2;i++) if ((LocMesh->Elements[i].E1==IdEdge)||(LocMesh->Elements[i].E2==IdEdge)||(LocMesh->Elements[i].E3==IdEdge)){
171  if (IdSubdiviedElement1==-1){
172  IdSubdiviedElement1=i;
173  IdNewElement1=LocMesh->NbElements-2;
174  if ((LocMesh->Elements[i].E1!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E1].V1==IdNgbhVortex1)||(LocMesh->Edges[LocMesh->Elements[i].E1].V2==IdNgbhVortex1))){
175  IdNgbhEdge1=LocMesh->Elements[i].E1;
176  if (LocMesh->Edges[LocMesh->Elements[i].E1].V1==IdNgbhVortex1) IdNgbhVortex3=LocMesh->Edges[LocMesh->Elements[i].E1].V2;
177  else IdNgbhVortex3=LocMesh->Edges[LocMesh->Elements[i].E1].V1;
178  }
179  if ((LocMesh->Elements[i].E2!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E2].V1==IdNgbhVortex1)||(LocMesh->Edges[LocMesh->Elements[i].E2].V2==IdNgbhVortex1))){
180  IdNgbhEdge1=LocMesh->Elements[i].E2;
181  if (LocMesh->Edges[LocMesh->Elements[i].E2].V1==IdNgbhVortex1) IdNgbhVortex3=LocMesh->Edges[LocMesh->Elements[i].E2].V2;
182  else IdNgbhVortex3=LocMesh->Edges[LocMesh->Elements[i].E2].V1;
183  }
184  if ((LocMesh->Elements[i].E3!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E3].V1==IdNgbhVortex1)||(LocMesh->Edges[LocMesh->Elements[i].E3].V2==IdNgbhVortex1))){
185  IdNgbhEdge1=LocMesh->Elements[i].E3;
186  if (LocMesh->Edges[LocMesh->Elements[i].E3].V1==IdNgbhVortex1) IdNgbhVortex3=LocMesh->Edges[LocMesh->Elements[i].E3].V2;
187  else IdNgbhVortex3=LocMesh->Edges[LocMesh->Elements[i].E3].V1;
188  }
189  if ((LocMesh->Elements[i].E1!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E1].V1==IdNgbhVortex2)||(LocMesh->Edges[LocMesh->Elements[i].E1].V2==IdNgbhVortex2))){
190  IdNgbhEdge3=LocMesh->Elements[i].E1;
191  }
192  if ((LocMesh->Elements[i].E2!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E2].V1==IdNgbhVortex2)||(LocMesh->Edges[LocMesh->Elements[i].E2].V2==IdNgbhVortex2))){
193  IdNgbhEdge3=LocMesh->Elements[i].E2;
194  }
195  if ((LocMesh->Elements[i].E3!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E3].V1==IdNgbhVortex2)||(LocMesh->Edges[LocMesh->Elements[i].E3].V2==IdNgbhVortex2))){
196  IdNgbhEdge3=LocMesh->Elements[i].E3;
197  }
198  }
199  else{
200  IdSubdiviedElement2=i;
201  IdNewElement2=LocMesh->NbElements-1;
202  if ((LocMesh->Elements[i].E1!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E1].V1==IdNgbhVortex1)||(LocMesh->Edges[LocMesh->Elements[i].E1].V2==IdNgbhVortex1))){
203  IdNgbhEdge2=LocMesh->Elements[i].E1;
204  if (LocMesh->Edges[LocMesh->Elements[i].E1].V1==IdNgbhVortex1) IdNgbhVortex4=LocMesh->Edges[LocMesh->Elements[i].E1].V2;
205  else IdNgbhVortex4=LocMesh->Edges[LocMesh->Elements[i].E1].V1;
206  }
207  if ((LocMesh->Elements[i].E2!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E2].V1==IdNgbhVortex1)||(LocMesh->Edges[LocMesh->Elements[i].E2].V2==IdNgbhVortex1))){
208  IdNgbhEdge2=LocMesh->Elements[i].E2;
209  if (LocMesh->Edges[LocMesh->Elements[i].E2].V1==IdNgbhVortex1) IdNgbhVortex4=LocMesh->Edges[LocMesh->Elements[i].E2].V2;
210  else IdNgbhVortex4=LocMesh->Edges[LocMesh->Elements[i].E2].V1;
211  }
212  if ((LocMesh->Elements[i].E3!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E3].V1==IdNgbhVortex1)||(LocMesh->Edges[LocMesh->Elements[i].E3].V2==IdNgbhVortex1))){
213  IdNgbhEdge2=LocMesh->Elements[i].E3;
214  if (LocMesh->Edges[LocMesh->Elements[i].E3].V1==IdNgbhVortex1) IdNgbhVortex4=LocMesh->Edges[LocMesh->Elements[i].E3].V2;
215  else IdNgbhVortex4=LocMesh->Edges[LocMesh->Elements[i].E3].V1;
216  }
217  if ((LocMesh->Elements[i].E1!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E1].V1==IdNgbhVortex2)||(LocMesh->Edges[LocMesh->Elements[i].E1].V2==IdNgbhVortex2))){
218  IdNgbhEdge4=LocMesh->Elements[i].E1;
219  }
220  if ((LocMesh->Elements[i].E2!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E2].V1==IdNgbhVortex2)||(LocMesh->Edges[LocMesh->Elements[i].E2].V2==IdNgbhVortex2))){
221  IdNgbhEdge4=LocMesh->Elements[i].E2;
222  }
223  if ((LocMesh->Elements[i].E3!=IdEdge)&&((LocMesh->Edges[LocMesh->Elements[i].E3].V1==IdNgbhVortex2)||(LocMesh->Edges[LocMesh->Elements[i].E3].V2==IdNgbhVortex2))){
224  IdNgbhEdge4=LocMesh->Elements[i].E3;
225  }
226  }
227  }
228  IdNewElement1=LocMesh->NbElements-2;
229  IdNewElement2=LocMesh->NbElements-1;
230 
231  //Coordinates of the new vertex
232  LocMesh->Vertexes[IdNewVortex].X=(LocMesh->Vertexes[IdNgbhVortex1].X+LocMesh->Vertexes[IdNgbhVortex2].X)/2.;
233  LocMesh->Vertexes[IdNewVortex].Y=(LocMesh->Vertexes[IdNgbhVortex1].Y+LocMesh->Vertexes[IdNgbhVortex2].Y)/2.;
234  LocMesh->Vertexes[IdNewVortex].Z=(LocMesh->Vertexes[IdNgbhVortex1].Z+LocMesh->Vertexes[IdNgbhVortex2].Z)/2.;
235 
236  //update the vertex identifiers of the edge ends
237  LocMesh->Edges[IdSubdiviedEdge].V1=IdNgbhVortex1;
238  LocMesh->Edges[IdSubdiviedEdge].V2=IdNewVortex;
239 
240  LocMesh->Edges[IdNewEdge1].V1=IdNewVortex;
241  LocMesh->Edges[IdNewEdge1].V2=IdNgbhVortex2;
242  LocMesh->Edges[IdNewEdge1].NbSubdiv=1;
243 
244  LocMesh->Edges[IdNewEdge2].V1=IdNewVortex;
245  LocMesh->Edges[IdNewEdge2].V2=IdNgbhVortex4;
246  LocMesh->Edges[IdNewEdge2].NbSubdiv=1;
247 
248  LocMesh->Edges[IdNewEdge3].V1=IdNewVortex;
249  LocMesh->Edges[IdNewEdge3].V2=IdNgbhVortex3;
250  LocMesh->Edges[IdNewEdge3].NbSubdiv=1;
251 
252 
253  //update the edges identifiers of the elements
254  LocMesh->Elements[IdSubdiviedElement1].E1=IdSubdiviedEdge;
255  LocMesh->Elements[IdSubdiviedElement1].E2=IdNewEdge3;
256  LocMesh->Elements[IdSubdiviedElement1].E3=IdNgbhEdge1;
257 
258  LocMesh->Elements[IdSubdiviedElement2].E1=IdNgbhEdge2;
259  LocMesh->Elements[IdSubdiviedElement2].E2=IdNewEdge2;
260  LocMesh->Elements[IdSubdiviedElement2].E3=IdSubdiviedEdge;
261 
262  LocMesh->Elements[IdNewElement1].E1=IdNewEdge1;
263  LocMesh->Elements[IdNewElement1].E2=IdNgbhEdge3;
264  LocMesh->Elements[IdNewElement1].E3=IdNewEdge3;
265 
266  LocMesh->Elements[IdNewElement2].E1=IdNgbhEdge4;
267  LocMesh->Elements[IdNewElement2].E2=IdNewEdge1;
268  LocMesh->Elements[IdNewElement2].E3=IdNewEdge2;
269 }
270 
271 //Project de coordinates of a mesh on a ellispoid of radius 'Rx,Ry,Rz'
272 void ProjectMeshEllipsoid(Mesh * LocMesh, double Rx, double Ry, double Rz){
273  int i;
274  double Xmax,Xmin,Ymax,Ymin,Zmax,Zmin,Xmean,Ymean,Zmean;
275  double LocNorm;
276 
277  //translate the mesh centre to the origin
278  Xmax=LocMesh->Vertexes[0].X; Xmin=LocMesh->Vertexes[0].X;
279  Ymax=LocMesh->Vertexes[0].Y; Ymin=LocMesh->Vertexes[0].Y;
280  Zmax=LocMesh->Vertexes[0].Z; Zmin=LocMesh->Vertexes[0].Z;
281 
282  for (i=1;i<LocMesh->NbVertexes;i++){
283  if (LocMesh->Vertexes[i].X>Xmax) Xmax=LocMesh->Vertexes[i].X;
284  if (LocMesh->Vertexes[i].X<Xmin) Xmin=LocMesh->Vertexes[i].X;
285  if (LocMesh->Vertexes[i].Y>Ymax) Ymax=LocMesh->Vertexes[i].Y;
286  if (LocMesh->Vertexes[i].Y<Ymin) Ymin=LocMesh->Vertexes[i].Y;
287  if (LocMesh->Vertexes[i].Z>Zmax) Zmax=LocMesh->Vertexes[i].Z;
288  if (LocMesh->Vertexes[i].Z<Zmin) Zmin=LocMesh->Vertexes[i].Z;
289  }
290 
291  Xmean=(Xmax+Xmin)/2;
292  Ymean=(Ymax+Ymin)/2;
293  Zmean=(Zmax+Zmin)/2;
294 
295  for (i=0;i<LocMesh->NbVertexes;i++){
296  LocMesh->Vertexes[i].X-=Xmean;
297  LocMesh->Vertexes[i].Y-=Ymean;
298  LocMesh->Vertexes[i].Z-=Zmean;
299  }
300 
301  //project each vertex on a ellipsoid of radius 'Rx,Ry,Rz'
302  for (i=0;i<LocMesh->NbVertexes;i++){
303  LocNorm=sqrt((pow(LocMesh->Vertexes[i].X,2.0)/pow(Rx,2.0))+
304  (pow(LocMesh->Vertexes[i].Y,2.0)/pow(Ry,2.0))+
305  (pow(LocMesh->Vertexes[i].Z,2.0)/pow(Rz,2.0)));
306  LocMesh->Vertexes[i].X=LocMesh->Vertexes[i].X/LocNorm;
307  LocMesh->Vertexes[i].Y=LocMesh->Vertexes[i].Y/LocNorm;
308  LocMesh->Vertexes[i].Z=LocMesh->Vertexes[i].Z/LocNorm;
309  }
310 }
311 
312 
313 //order the edges of all element so that this order is clockwise if we see it from outside the shape
314 void MakeClockwiseOrder(Mesh * LocMesh){
315  int i,j,jswap;
316  int Ve1,Ed1,Ve2,Ed2,Ve3,Ed3; //Order: Vertex 1 -> Edge1 -> Vertex 2 -> Edge 2 -> Vertex 3 -> Edge 3 -> Vertex 1 ...
317  int temp;
318  double x1,x2,x3,y1,y2,y3,z1,z2,z3; //coordinates
319  double u1,u2,u3,v1,v2,v3; //vectors
320  double vp1,vp2,vp3; //cross product
321  double tempDbl;
322 
323  //test all elements of the mesh
324  for (i=0;i<LocMesh->NbElements;i++){
325  //order the vertexes and edges
326  Ed1=LocMesh->Elements[i].E1;
327  Ed2=LocMesh->Elements[i].E2;
328  Ed3=LocMesh->Elements[i].E3;
329 
330  if ((LocMesh->Edges[Ed1].V2!=LocMesh->Edges[Ed2].V1)&&(LocMesh->Edges[Ed1].V2!=LocMesh->Edges[Ed2].V2)){ //swap edges
331  temp=Ed2; Ed2=Ed3; Ed3=temp;
332  }
333 
334  if (LocMesh->Edges[Ed1].V2!=LocMesh->Edges[Ed2].V1){ //swap vertexes
335  temp=LocMesh->Edges[Ed2].V1; LocMesh->Edges[Ed2].V1=LocMesh->Edges[Ed2].V2; LocMesh->Edges[Ed2].V2=temp;
336  if (LocMesh->Edges[Ed2].NbSubdiv>1)
337  for (j=0;j<=LocMesh->Edges[Ed2].NbSubdiv/2;j++){
338  jswap=LocMesh->Edges[Ed2].NbSubdiv-j;
339  tempDbl=LocMesh->Edges[Ed2].SubiVertex[j].X; LocMesh->Edges[Ed2].SubiVertex[j].X=LocMesh->Edges[Ed2].SubiVertex[jswap].X; LocMesh->Edges[Ed2].SubiVertex[jswap].X=tempDbl;
340  tempDbl=LocMesh->Edges[Ed2].SubiVertex[j].Y; LocMesh->Edges[Ed2].SubiVertex[j].Y=LocMesh->Edges[Ed2].SubiVertex[jswap].Y; LocMesh->Edges[Ed2].SubiVertex[jswap].Y=tempDbl;
341  tempDbl=LocMesh->Edges[Ed2].SubiVertex[j].Z; LocMesh->Edges[Ed2].SubiVertex[j].Z=LocMesh->Edges[Ed2].SubiVertex[jswap].Z; LocMesh->Edges[Ed2].SubiVertex[jswap].Z=tempDbl;
342  }
343  }
344 
345  if (LocMesh->Edges[Ed2].V2!=LocMesh->Edges[Ed3].V1){ //swap vertexes
346  temp=LocMesh->Edges[Ed3].V1; LocMesh->Edges[Ed3].V1=LocMesh->Edges[Ed3].V2; LocMesh->Edges[Ed3].V2=temp;
347  if (LocMesh->Edges[Ed3].NbSubdiv>1)
348  for (j=0;j<=LocMesh->Edges[Ed3].NbSubdiv/2;j++){
349  jswap=LocMesh->Edges[Ed3].NbSubdiv-j;
350  tempDbl=LocMesh->Edges[Ed3].SubiVertex[j].X; LocMesh->Edges[Ed3].SubiVertex[j].X=LocMesh->Edges[Ed3].SubiVertex[jswap].X; LocMesh->Edges[Ed3].SubiVertex[jswap].X=tempDbl;
351  tempDbl=LocMesh->Edges[Ed3].SubiVertex[j].Y; LocMesh->Edges[Ed3].SubiVertex[j].Y=LocMesh->Edges[Ed3].SubiVertex[jswap].Y; LocMesh->Edges[Ed3].SubiVertex[jswap].Y=tempDbl;
352  tempDbl=LocMesh->Edges[Ed3].SubiVertex[j].Z; LocMesh->Edges[Ed3].SubiVertex[j].Z=LocMesh->Edges[Ed3].SubiVertex[jswap].Z; LocMesh->Edges[Ed3].SubiVertex[jswap].Z=tempDbl;
353  }
354  }
355 
356  if (LocMesh->Edges[Ed3].V2!=LocMesh->Edges[Ed1].V1){
357  printf("This is not good!!!\n");
358  }
359 
360  //invert the order if necessary
361  x1=LocMesh->Vertexes[LocMesh->Edges[Ed1].V1].X; y1=LocMesh->Vertexes[LocMesh->Edges[Ed1].V1].Y; z1=LocMesh->Vertexes[LocMesh->Edges[Ed1].V1].Z;
362  x2=LocMesh->Vertexes[LocMesh->Edges[Ed1].V2].X; y2=LocMesh->Vertexes[LocMesh->Edges[Ed1].V2].Y; z2=LocMesh->Vertexes[LocMesh->Edges[Ed1].V2].Z;
363  x3=LocMesh->Vertexes[LocMesh->Edges[Ed3].V1].X; y3=LocMesh->Vertexes[LocMesh->Edges[Ed3].V1].Y; z3=LocMesh->Vertexes[LocMesh->Edges[Ed3].V1].Z;
364 
365  u1=x2-x1; u2=y2-y1; u3=z2-z1;
366  v1=x3-x1; v2=y3-y1; v3=z3-z1;
367 
368  vp1=u2*v3-u3*v2;
369  vp2=u3*v1-u1*v3;
370  vp3=u1*v2-u2*v1;
371 
372  //printf("%3.2lf %3.2lf %3.2lf | %3.2lf %3.2lf %3.2lf | %3.2lf %3.2lf %3.2lf || %3.2lf %3.2lf %3.2lf || %3.2lf %3.2lf %3.2lf | %lf\n",x1,y1,z1,x2,y2,z2,x3,y3,z3,u1,u2,u3,v1,v2,v3,x1*vp1+y1*vp2+z1*vp3);
373  //printf("%3.2lf %3.2lf %3.2lf | %3.2lf %3.2lf %3.2lf | %3.2lf\n",x1,y1,z1,vp1,vp2,vp3,x1*vp1+y1*vp2+z1*vp3);
374 
375  if (x1*vp1+y1*vp2+z1*vp3<0){//we consider that the origin is within the spheric shape so the vector (x1,y1,z1) points out of the shape
376  //swap edges
377  temp=Ed3; Ed3=Ed2; Ed2=temp;
378  //swap vertexes
379  temp=LocMesh->Edges[Ed1].V1; LocMesh->Edges[Ed1].V1=LocMesh->Edges[Ed1].V2; LocMesh->Edges[Ed1].V2=temp;
380  if (LocMesh->Edges[Ed1].NbSubdiv>1)
381  for (j=0;j<=LocMesh->Edges[Ed1].NbSubdiv/2;j++){
382  jswap=LocMesh->Edges[Ed1].NbSubdiv-j;
383  tempDbl=LocMesh->Edges[Ed1].SubiVertex[j].X; LocMesh->Edges[Ed1].SubiVertex[j].X=LocMesh->Edges[Ed1].SubiVertex[jswap].X; LocMesh->Edges[Ed1].SubiVertex[jswap].X=tempDbl;
384  tempDbl=LocMesh->Edges[Ed1].SubiVertex[j].Y; LocMesh->Edges[Ed1].SubiVertex[j].Y=LocMesh->Edges[Ed1].SubiVertex[jswap].Y; LocMesh->Edges[Ed1].SubiVertex[jswap].Y=tempDbl;
385  tempDbl=LocMesh->Edges[Ed1].SubiVertex[j].Z; LocMesh->Edges[Ed1].SubiVertex[j].Z=LocMesh->Edges[Ed1].SubiVertex[jswap].Z; LocMesh->Edges[Ed1].SubiVertex[jswap].Z=tempDbl;
386  }
387  }
388  //new order of the elements
389  LocMesh->Elements[i].E1=Ed1;
390  LocMesh->Elements[i].E2=Ed2;
391  LocMesh->Elements[i].E3=Ed3;
392  }
393 
394 }
395 
396 //subdivide into 2 edges all edges in the mesh 'LocMesh'
397 void RefineMesh(Mesh * LocMesh){
398  int i,j,k;
399  int * order;
400  double SqLengthI,SqLengthJ;
401  int OrigNbEdges;
402 
403  order=(int*)malloc(LocMesh->NbEdges*sizeof(int));
404 
405  for (i=0;i<LocMesh->NbEdges;i++) order[i]=i;
406 
407  for (i=0;i<LocMesh->NbEdges-1;i++)
408  for (j=i+1;j<LocMesh->NbEdges;j++){
409  SqLengthI=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].X-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].X,2.0);
410  SqLengthI+=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].Y-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].Y,2.0);
411  SqLengthI+=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].Z-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].Z,2.0);
412  SqLengthJ=pow(LocMesh->Vertexes[LocMesh->Edges[order[j]].V1].X-LocMesh->Vertexes[LocMesh->Edges[order[j]].V2].X,2.0);
413  SqLengthJ+=pow(LocMesh->Vertexes[LocMesh->Edges[order[j]].V1].Y-LocMesh->Vertexes[LocMesh->Edges[order[j]].V2].Y,2.0);
414  SqLengthJ+=pow(LocMesh->Vertexes[LocMesh->Edges[order[j]].V1].Z-LocMesh->Vertexes[LocMesh->Edges[order[j]].V2].Z,2.0);
415 
416  if (SqLengthI<SqLengthJ){
417  k=order[i];
418  order[i]=order[j];
419  order[j]=k;
420  }
421  }
422 
423  OrigNbEdges=LocMesh->NbEdges;
424  for (i=0;i<OrigNbEdges;i++) CutEdge(LocMesh,order[i]);
425 }
426 
427 
428 //Project de coordinates of a mesh on an ellipsoid of radius 'Rx,Ry,Rz' and refine the mesh until
429 //the longest edge has a length smaller than 'MaxEdgeLength'.
430 //If NbSubdiv>1, each edge is curved into 'NbSubdiv' parts in the end
431 void ProjectMeshEllipsoidAndRefine(Mesh * LocMesh, double Rx,double Ry,double Rz,double MaxEdgeLength,int NbSubdiv){
432  int i,j,k;
433  int * order;
434  double SqLengthI,SqLengthJ;
435  int OrigNbEdges;
436  int OK;
437  double LocNorm;
438  double dj,dNbSubdiv;
439  int jm1;
440  double EdgeV1_X,EdgeV1_Y,EdgeV1_Z,EdgeV2_X,EdgeV2_Y,EdgeV2_Z;
441 
442  //initialize
443  order=(int*)malloc(LocMesh->NbEdges*sizeof(int));
444  for (i=0;i<LocMesh->NbEdges;i++) order[i]=i;
445  OK=1;
446 
447  //mesh projection on an ellipsoid
448  ProjectMeshEllipsoid(LocMesh,Rx,Ry,Rz);
449 
450  //Big loop to subdivide the mesh until each edge has a length smaller than 'MaxEdgeLength'
451  while (OK==1){
452  OK=0;
453 
454  order=(int*)realloc(order,LocMesh->NbEdges*sizeof(int));
455  for (i=0;i<LocMesh->NbEdges;i++) order[i]=i;
456 
457  //order the edges as a function of their length
458  for (i=0;i<LocMesh->NbEdges-1;i++) for (j=i+1;j<LocMesh->NbEdges;j++){
459  SqLengthI=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].X-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].X,2.0);
460  SqLengthI+=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].Y-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].Y,2.0);
461  SqLengthI+=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].Z-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].Z,2.0);
462  SqLengthJ=pow(LocMesh->Vertexes[LocMesh->Edges[order[j]].V1].X-LocMesh->Vertexes[LocMesh->Edges[order[j]].V2].X,2.0);
463  SqLengthJ+=pow(LocMesh->Vertexes[LocMesh->Edges[order[j]].V1].Y-LocMesh->Vertexes[LocMesh->Edges[order[j]].V2].Y,2.0);
464  SqLengthJ+=pow(LocMesh->Vertexes[LocMesh->Edges[order[j]].V1].Z-LocMesh->Vertexes[LocMesh->Edges[order[j]].V2].Z,2.0);
465 
466  if (SqLengthI<SqLengthJ){ k=order[i]; order[i]=order[j]; order[j]=k;}
467  }
468 
469  //subdivide the edges if their length is larger than 'MaxEdgeLength'
470  OrigNbEdges=LocMesh->NbEdges;
471  for (i=0;i<OrigNbEdges;i++){
472  SqLengthI=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].X-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].X,2.0);
473  SqLengthI+=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].Y-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].Y,2.0);
474  SqLengthI+=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].Z-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].Z,2.0);
475 
476  if (sqrt(SqLengthI)>MaxEdgeLength){
477  //subdvision
478  OK=1;
479  CutEdge(LocMesh,order[i]); //rem: 'LocMesh' is changed here. In particular LocMesh->NbEdges is larger.
480  }
481  }
482 
483  //projection of the new edges on an ellipsoid
484  ProjectMeshEllipsoid(LocMesh,Rx,Ry,Rz);
485  }
486 
487 
488 
489  //subdivide the edges end project them on the ellipsoid if required
490  if (NbSubdiv>1){
491  for (i=0;i<LocMesh->NbEdges;i++){
492  //init
493  LocMesh->Edges[i].NbSubdiv=NbSubdiv;
494  LocMesh->Edges[i].SubiVertex=(Vertex *)malloc((LocMesh->Edges[i].NbSubdiv+1)*sizeof(Vertex));
495 
496  //coordinates of the subdivised edge
497  EdgeV1_X=LocMesh->Vertexes[LocMesh->Edges[i].V1].X;
498  EdgeV1_Y=LocMesh->Vertexes[LocMesh->Edges[i].V1].Y;
499  EdgeV1_Z=LocMesh->Vertexes[LocMesh->Edges[i].V1].Z;
500  EdgeV2_X=LocMesh->Vertexes[LocMesh->Edges[i].V2].X;
501  EdgeV2_Y=LocMesh->Vertexes[LocMesh->Edges[i].V2].Y;
502  EdgeV2_Z=LocMesh->Vertexes[LocMesh->Edges[i].V2].Z;
503  dNbSubdiv=(double)LocMesh->Edges[i].NbSubdiv;
504 
505  for (j=0;j<LocMesh->Edges[i].NbSubdiv+1;j++){
506  dj=(double)j;
507 
508  //interpolation
509  LocMesh->Edges[i].SubiVertex[j].X=EdgeV1_X*(dNbSubdiv-dj)/dNbSubdiv+EdgeV2_X*dj/dNbSubdiv;
510  LocMesh->Edges[i].SubiVertex[j].Y=EdgeV1_Y*(dNbSubdiv-dj)/dNbSubdiv+EdgeV2_Y*dj/dNbSubdiv;
511  LocMesh->Edges[i].SubiVertex[j].Z=EdgeV1_Z*(dNbSubdiv-dj)/dNbSubdiv+EdgeV2_Z*dj/dNbSubdiv;
512 
513  //projection on an ellipsoid
514  LocNorm=sqrt((pow(LocMesh->Edges[i].SubiVertex[j].X,2.0)/pow(Rx,2.0))+
515  (pow(LocMesh->Edges[i].SubiVertex[j].Y,2.0)/pow(Ry,2.0))+
516  (pow(LocMesh->Edges[i].SubiVertex[j].Z,2.0)/pow(Rz,2.0)));
517  LocMesh->Edges[i].SubiVertex[j].X=LocMesh->Edges[i].SubiVertex[j].X/LocNorm;
518  LocMesh->Edges[i].SubiVertex[j].Y=LocMesh->Edges[i].SubiVertex[j].Y/LocNorm;
519  LocMesh->Edges[i].SubiVertex[j].Z=LocMesh->Edges[i].SubiVertex[j].Z/LocNorm;
520  }
521  }
522  }
523 }
524 
525 
526 //Extract a submesh of 'LocMesh' in the bounding box [Xmin,Xmax],[Ymin,Ymax],[Zmin,Zmax] (NEW)
527 void ExtractSubmesh(Mesh * LocMesh,double Xmin,double Xmax,double Ymin,double Ymax,double Zmin,double Zmax){
528  int i,j,NbCurvedEdge;
529  int * VertexesToKeep;
530  int * EdgesToKeep;
531  int * ElementsToKeep;
532  int NewValue;
533  int NewI;
534  double epsilon;
535  Mesh TempMesh;
536 
537  //1) allocations and initializations
538  epsilon=0.001;
539 
540  VertexesToKeep=(int*)malloc(LocMesh->NbVertexes*sizeof(int));
541  for (i=0;i<LocMesh->NbVertexes;i++) VertexesToKeep[i]=0;
542 
543  EdgesToKeep=(int*)malloc(LocMesh->NbEdges*sizeof(int));
544  for (i=0;i<LocMesh->NbEdges;i++) EdgesToKeep[i]=0;
545 
546  ElementsToKeep=(int*)malloc(LocMesh->NbElements*sizeof(int));
547  for (i=0;i<LocMesh->NbElements;i++) ElementsToKeep[i]=0;
548 
549  //2) find the vertexes to keep
550  for (i=0;i<LocMesh->NbVertexes;i++)
551  if ((LocMesh->Vertexes[i].X>Xmin-epsilon)&&(LocMesh->Vertexes[i].X<Xmax+epsilon)&&
552  (LocMesh->Vertexes[i].Y>Ymin-epsilon)&&(LocMesh->Vertexes[i].Y<Ymax+epsilon)&&
553  (LocMesh->Vertexes[i].Z>Zmin-epsilon)&&(LocMesh->Vertexes[i].Z<Zmax+epsilon)) VertexesToKeep[i]=1;
554 
555  //3) find the edges to keep
556  for (i=0;i<LocMesh->NbEdges;i++) if ((VertexesToKeep[LocMesh->Edges[i].V1]==1)&&(VertexesToKeep[LocMesh->Edges[i].V2]==1))
557  EdgesToKeep[i]=1;
558 
559  //4) find the elements to keep
560  for (i=0;i<LocMesh->NbElements;i++) if ((EdgesToKeep[LocMesh->Elements[i].E1]==1)&&(EdgesToKeep[LocMesh->Elements[i].E2]==1)&&(EdgesToKeep[LocMesh->Elements[i].E3]==1))
561  ElementsToKeep[i]=1;
562 
563 
564  //5) extraction of the submesh
565 
566  //5.1) new IDs (ID+1 actually)
567  NewValue=1;
568  for (i=0;i<LocMesh->NbVertexes;i++) if (VertexesToKeep[i]==1){
569  VertexesToKeep[i]=NewValue;
570  NewValue++;
571  }
572 
573  NewValue=1;
574  for (i=0;i<LocMesh->NbEdges;i++) if (EdgesToKeep[i]==1){
575  EdgesToKeep[i]=NewValue;
576  NewValue++;
577  }
578 
579  NewValue=1;
580  for (i=0;i<LocMesh->NbElements;i++) if (ElementsToKeep[i]==1){
581  ElementsToKeep[i]=NewValue;
582  NewValue++;
583  }
584 
585  //5.2) copy the mesh in a temporary mesh
586  TempMesh.NbVertexes=LocMesh->NbVertexes;
587  TempMesh.Vertexes=(Vertex *)malloc(TempMesh.NbVertexes*sizeof(Vertex));
588  for (i=0;i<TempMesh.NbVertexes;i++){
589  TempMesh.Vertexes[i].X=LocMesh->Vertexes[i].X;
590  TempMesh.Vertexes[i].Y=LocMesh->Vertexes[i].Y;
591  TempMesh.Vertexes[i].Z=LocMesh->Vertexes[i].Z;
592  }
593 
594  TempMesh.NbEdges=LocMesh->NbEdges;
595  TempMesh.Edges=(Edge *)malloc(TempMesh.NbEdges*sizeof(Edge));
596  for (i=0;i<TempMesh.NbEdges;i++){
597  TempMesh.Edges[i].V1=LocMesh->Edges[i].V1;
598  TempMesh.Edges[i].V2=LocMesh->Edges[i].V2;
599  }
600  for (i=0;i<LocMesh->NbEdges;i++) if (LocMesh->Edges[i].NbSubdiv>1){
601  TempMesh.Edges[i].NbSubdiv=LocMesh->Edges[i].NbSubdiv;
602  TempMesh.Edges[i].SubiVertex=(Vertex *)malloc((TempMesh.Edges[i].NbSubdiv+1)*sizeof(Vertex));
603  for (j=0;j<TempMesh.Edges[i].NbSubdiv+1;j++){
604  TempMesh.Edges[i].SubiVertex[j].X=LocMesh->Edges[i].SubiVertex[j].X;
605  TempMesh.Edges[i].SubiVertex[j].Y=LocMesh->Edges[i].SubiVertex[j].Y;
606  TempMesh.Edges[i].SubiVertex[j].Z=LocMesh->Edges[i].SubiVertex[j].Z;
607  }
608  }
609 
610  TempMesh.NbElements=LocMesh->NbElements;
611  TempMesh.Elements=(Element *)malloc(TempMesh.NbElements*sizeof(Element));
612  for (i=0;i<TempMesh.NbElements;i++){
613  TempMesh.Elements[i].E1=LocMesh->Elements[i].E1;
614  TempMesh.Elements[i].E2=LocMesh->Elements[i].E2;
615  TempMesh.Elements[i].E3=LocMesh->Elements[i].E3;
616  }
617 
618  //5.3) change the mesh
619  NewI=0;
620  for (i=0;i<TempMesh.NbVertexes;i++) if (VertexesToKeep[i]>0){
621  LocMesh->Vertexes[NewI].X=TempMesh.Vertexes[i].X;
622  LocMesh->Vertexes[NewI].Y=TempMesh.Vertexes[i].Y;
623  LocMesh->Vertexes[NewI].Z=TempMesh.Vertexes[i].Z;
624  //printf("Vertex %d - %f %f %f\n",NewI,LocMesh->Vertexes[NewI].X,LocMesh->Vertexes[NewI].Y,LocMesh->Vertexes[NewI].Z);
625  NewI++;
626  }
627  LocMesh->NbVertexes=NewI;
628 
629  NewI=0;
630  for (i=0;i<TempMesh.NbEdges;i++) if (EdgesToKeep[i]>0){
631  LocMesh->Edges[NewI].V1=VertexesToKeep[TempMesh.Edges[i].V1]-1;
632  LocMesh->Edges[NewI].V2=VertexesToKeep[TempMesh.Edges[i].V2]-1;
633  if (TempMesh.Edges[i].NbSubdiv>1) for (j=0;j<TempMesh.Edges[i].NbSubdiv+1;j++){
634  LocMesh->Edges[NewI].SubiVertex[j].X=TempMesh.Edges[i].SubiVertex[j].X;
635  LocMesh->Edges[NewI].SubiVertex[j].Y=TempMesh.Edges[i].SubiVertex[j].Y;
636  LocMesh->Edges[NewI].SubiVertex[j].Z=TempMesh.Edges[i].SubiVertex[j].Z;
637  } //we suppose that all edges have the same number of subdivisions
638  //printf("Edge %d - %d %d\n",NewI,LocMesh->Edges[NewI].V1,LocMesh->Edges[NewI].V2);
639  NewI++;
640  }
641  LocMesh->NbEdges=NewI;
642 
643  NewI=0;
644  for (i=0;i<TempMesh.NbElements;i++) if (ElementsToKeep[i]>0){
645  LocMesh->Elements[NewI].E1=EdgesToKeep[TempMesh.Elements[i].E1]-1;
646  LocMesh->Elements[NewI].E2=EdgesToKeep[TempMesh.Elements[i].E2]-1;
647  LocMesh->Elements[NewI].E3=EdgesToKeep[TempMesh.Elements[i].E3]-1;
648  //printf("Element %d - %d %d %d \n",NewI,LocMesh->Elements[NewI].E1,LocMesh->Elements[NewI].E2,LocMesh->Elements[NewI].E3);
649  NewI++;
650  }
651  LocMesh->NbElements=NewI;
652  }
653 
654 
655 
656 //write the geometry of mesh contained in a Mesh structure in a XML-Nektar file (NEW STUFFS THERE)
657 void WriteMesh(Mesh * LocMesh,char FileName[256]){
658  FILE * XmlMeshGeomFile;
659  int i,j,NbCurvedEdge;
660  int * BoundaryEdgesList;
661  int FirstBoundary;
662 
663  //to have a proper order of the edges within the elements
664  MakeClockwiseOrder(LocMesh);
665 
666  //open file
667  XmlMeshGeomFile = fopen(FileName,"w");
668 
669  //write header
670  fprintf(XmlMeshGeomFile,"<?xml version=\"1.0\" encoding=\"utf-8\"?>\n");
671  fprintf(XmlMeshGeomFile,"\n");
672  fprintf(XmlMeshGeomFile,"<NEKTAR>\n");
673  fprintf(XmlMeshGeomFile,"<!-- Embed a 2-dimensional object in a 2-dimensional space -->\n");
674  fprintf(XmlMeshGeomFile,"<!-- DIM <= SPACE -->\n");
675  fprintf(XmlMeshGeomFile,"<!-- This provides a method of optimizing code for a 1-D curve embedded in 3-space. -->\n");
676  fprintf(XmlMeshGeomFile,"<GEOMETRY DIM=\"2\" SPACE=\"3\">\n");
677  fprintf(XmlMeshGeomFile,"\n");
678 
679  //write vertex
680  fprintf(XmlMeshGeomFile," <VERTEX>\n");
681  fprintf(XmlMeshGeomFile," <!-- Always must have four values per entry. -->\n");
682  for (i=0;i<LocMesh->NbVertexes;i++)
683  fprintf(XmlMeshGeomFile," <V ID=\"%d\"> %f %f %f </V>\n",i,LocMesh->Vertexes[i].X,LocMesh->Vertexes[i].Y,LocMesh->Vertexes[i].Z);
684  fprintf(XmlMeshGeomFile," </VERTEX>\n");
685  fprintf(XmlMeshGeomFile," \n");
686 
687  //write edges
688  fprintf(XmlMeshGeomFile," <EDGE>\n");
689  fprintf(XmlMeshGeomFile," <!--Edges are vertex pairs -->\n");
690  for (i=0;i<LocMesh->NbEdges;i++)
691  fprintf(XmlMeshGeomFile," <E ID=\"%d\"> %d %d </E>\n",i,LocMesh->Edges[i].V1,LocMesh->Edges[i].V2);
692  fprintf(XmlMeshGeomFile," </EDGE>\n");
693  fprintf(XmlMeshGeomFile," \n");
694 
695  //write Element
696  fprintf(XmlMeshGeomFile," <ELEMENT>\n");
697  for (i=0;i<LocMesh->NbElements;i++)
698  fprintf(XmlMeshGeomFile," <T ID=\"%d\"> %d %d %d </T>\n",i,LocMesh->Elements[i].E1,LocMesh->Elements[i].E2,LocMesh->Elements[i].E3);
699  fprintf(XmlMeshGeomFile," </ELEMENT>\n");
700  fprintf(XmlMeshGeomFile," \n");
701 
702 
703 
704  //write curved elements
705  fprintf(XmlMeshGeomFile," <CURVED>\n");
706  NbCurvedEdge=0;
707  for (i=0;i<LocMesh->NbEdges;i++) if (LocMesh->Edges[i].NbSubdiv>1){
708  fprintf(XmlMeshGeomFile," <E ID=\"%d\" EDGEID=\"%d\" TYPE=\"PolyEvenlySpaced\" NUMPOINTS=\"%d\">\n",NbCurvedEdge,i,LocMesh->Edges[i].NbSubdiv+1);
709  for (j=0;j<LocMesh->Edges[i].NbSubdiv+1;j++){
710  fprintf(XmlMeshGeomFile," %lf %lf %lf\n",LocMesh->Edges[i].SubiVertex[j].X,LocMesh->Edges[i].SubiVertex[j].Y,LocMesh->Edges[i].SubiVertex[j].Z);
711  }
712  fprintf(XmlMeshGeomFile," </E>\n\n");
713  NbCurvedEdge++;
714  }
715  fprintf(XmlMeshGeomFile," </CURVED>\n");
716 
717 
718  //find the boundaries (points at 0 in BoundaryEdgesList)
719  BoundaryEdgesList=(int *)malloc(LocMesh->NbEdges*sizeof(int));
720 
721  for (i=0;i<LocMesh->NbEdges;i++) BoundaryEdgesList[i]=-1;
722 
723  for (i=0;i<LocMesh->NbElements;i++){
724  if (BoundaryEdgesList[LocMesh->Elements[i].E1]==1) printf("Problem between the elements and the edges\n");
725  if (BoundaryEdgesList[LocMesh->Elements[i].E2]==1) printf("Problem between the elements and the edges\n");
726  if (BoundaryEdgesList[LocMesh->Elements[i].E3]==1) printf("Problem between the elements and the edges\n");
727  if (BoundaryEdgesList[LocMesh->Elements[i].E1]==0) BoundaryEdgesList[LocMesh->Elements[i].E1]=1;
728  if (BoundaryEdgesList[LocMesh->Elements[i].E2]==0) BoundaryEdgesList[LocMesh->Elements[i].E2]=1;
729  if (BoundaryEdgesList[LocMesh->Elements[i].E3]==0) BoundaryEdgesList[LocMesh->Elements[i].E3]=1;
730  if (BoundaryEdgesList[LocMesh->Elements[i].E1]==-1) BoundaryEdgesList[LocMesh->Elements[i].E1]=0;
731  if (BoundaryEdgesList[LocMesh->Elements[i].E2]==-1) BoundaryEdgesList[LocMesh->Elements[i].E2]=0;
732  if (BoundaryEdgesList[LocMesh->Elements[i].E3]==-1) BoundaryEdgesList[LocMesh->Elements[i].E3]=0;
733  }
734 
735  //write the end (could be more evolved)
736  fprintf(XmlMeshGeomFile,"<!-- V - vertex, E - edge, F - face, L - element -->\n");
737  fprintf(XmlMeshGeomFile," <COMPOSITE>\n");
738  fprintf(XmlMeshGeomFile," <C ID=\"0\"> T[0-%d]\n",LocMesh->NbElements-1);
739  fprintf(XmlMeshGeomFile," </C>\n");
740  fprintf(XmlMeshGeomFile," <C ID=\"1\"> E[");
741  FirstBoundary=1;
742  for (i=0;i<LocMesh->NbEdges;i++) if (BoundaryEdgesList[i]==0){
743  if (FirstBoundary==1) FirstBoundary=0;
744  else fprintf(XmlMeshGeomFile,",");
745  fprintf(XmlMeshGeomFile,"%d",i);
746  }
747  fprintf(XmlMeshGeomFile,"]\n");
748  fprintf(XmlMeshGeomFile," </C>\n");
749  fprintf(XmlMeshGeomFile," </COMPOSITE>\n");
750  fprintf(XmlMeshGeomFile," \n");
751  fprintf(XmlMeshGeomFile," <DOMAIN> C[0] </DOMAIN>\n");
752  fprintf(XmlMeshGeomFile," \n");
753 
754  //EOF
755  fprintf(XmlMeshGeomFile,"</GEOMETRY>\n");
756 
757  //close file
758  fclose(XmlMeshGeomFile);
759 }
760 
761 
762 
763 
764 
765 // ---------------------------------------------------------------------------------------
766 // main function
767 // ---------------------------------------------------------------------------------------
768 
769 
770 void usage()
771 {
772  printf("Usage: makeEllipsoidNektar <options>\n");
773  printf("Where <options> are one or more of the following:\n");
774  printf("\t<-NameFile n> Name of the output file containg the ellipsoid (default=\"toto.xml\")\n");
775  printf("\t<-MaxLength n> Maximum length allowed for an edge (default=0.5)\n");
776  printf("\t<-NbSubdiv n> Number of subdivisions of each edge to make it curved (default=1)\n");
777  printf("\t<-RadiusX n> Radius of the ellipsoid in the X direction (default=1.)\n");
778  printf("\t<-RadiusY n> Radius of the ellipsoid in the Y direction (default=1.)\n");
779  printf("\t<-RadiusZ n> Radius of the ellipsoid in the Z direction (default=1.)\n");
780  printf("\t<-MinX n> Bounding box - min X (default=-1000.)\n");
781  printf("\t<-MaxX n> Bounding box - max X (default=1000.)\n");
782  printf("\t<-MinY n> Bounding box - min Y (default=-1000.)\n");
783  printf("\t<-MaxY n> Bounding box - max Y (default=1000.)\n");
784  printf("\t<-MinZ n> Bounding box - min Z (default=-1000.)\n");
785  printf("\t<-MaxZ n> Bounding box - max Z (default=1000.)\n");
786  exit(1);
787 }
788 
789 int main(int argc, char **argv)
790 {
791  int ok,i;
792  char *output_name = "toto.xml";
793  double RadiusX=1.;
794  double RadiusY=1.;
795  double RadiusZ=1.;
796  double MinX=-1000;
797  double MaxX=1000;
798  double MinY=-1000;
799  double MaxY=1000;
800  double MinZ=-1000;
801  double MaxZ=1000;
802  double Prec=0.5;
803  int NbSubdiv=1;
804  Mesh LocMesh;
805 
806 
807  //1) Parse parameters
808  while (argc > 1) {
809  ok = 0;
810  if ((ok == 0) && (strcmp(argv[1], "help") == 0)) {
811  usage();
812  }
813  if ((ok == 0) && (strcmp(argv[1], "-NameFile") == 0)) {
814  argc--;
815  argv++;
816  output_name = argv[1];
817  argc--;
818  argv++;
819  ok = 1;
820  }
821  if ((ok == 0) && (strcmp(argv[1], "-RadiusX") == 0)) {
822  argc--;
823  argv++;
824  RadiusX = atof(argv[1]);
825  if (RadiusX<0) RadiusX=1.;
826  argc--;
827  argv++;
828  ok = 1;
829  }
830  if ((ok == 0) && (strcmp(argv[1], "-RadiusY") == 0)) {
831  argc--;
832  argv++;
833  RadiusY = atof(argv[1]);
834  if (RadiusY<0) RadiusY=1.;
835  argc--;
836  argv++;
837  ok = 1;
838  }
839  if ((ok == 0) && (strcmp(argv[1], "-RadiusZ") == 0)) {
840  argc--;
841  argv++;
842  RadiusZ = atof(argv[1]);
843  if (RadiusZ<0) RadiusZ=1.;
844  argc--;
845  argv++;
846  ok = 1;
847  }
848  if ((ok == 0) && (strcmp(argv[1], "-MinX") == 0)) {
849  argc--;
850  argv++;
851  MinX = atof(argv[1]);
852  argc--;
853  argv++;
854  ok = 1;
855  }
856  if ((ok == 0) && (strcmp(argv[1], "-MaxX") == 0)) {
857  argc--;
858  argv++;
859  MaxX = atof(argv[1]);
860  argc--;
861  argv++;
862  ok = 1;
863  }
864  if ((ok == 0) && (strcmp(argv[1], "-MinY") == 0)) {
865  argc--;
866  argv++;
867  MinY = atof(argv[1]);
868  argc--;
869  argv++;
870  ok = 1;
871  }
872  if ((ok == 0) && (strcmp(argv[1], "-MaxY") == 0)) {
873  argc--;
874  argv++;
875  MaxY = atof(argv[1]);
876  argc--;
877  argv++;
878  ok = 1;
879  }
880  if ((ok == 0) && (strcmp(argv[1], "-MinZ") == 0)) {
881  argc--;
882  argv++;
883  MinZ = atof(argv[1]);
884  argc--;
885  argv++;
886  ok = 1;
887  }
888  if ((ok == 0) && (strcmp(argv[1], "-MaxZ") == 0)) {
889  argc--;
890  argv++;
891  MaxZ = atof(argv[1]);
892  argc--;
893  argv++;
894  ok = 1;
895  }
896  if ((ok == 0) && (strcmp(argv[1], "-MaxLength") == 0)) {
897  argc--;
898  argv++;
899  Prec = atof(argv[1]);
900  if ((Prec<(RadiusX+RadiusY+RadiusZ)*0.001/3.)) Prec=(RadiusX+RadiusY+RadiusZ)*0.001/3;
901  argc--;
902  argv++;
903  ok = 1;
904  }
905  if ((ok == 0) && (strcmp(argv[1], "-NbSubdiv") == 0)) {
906  argc--;
907  argv++;
908  NbSubdiv = atoi(argv[1]);
909  if ((NbSubdiv<1)) NbSubdiv=1;
910  argc--;
911  argv++;
912  ok = 1;
913  }
914  if (ok == 0) {
915  printf("Unknown option: \n");
916  usage();
917  }
918  }
919 
920  printf("Output file: %s / Radius: %lf / Radius Y: %lf / Radius Z: %lf / Precision: %lf / Subdivisions: %d\n",output_name,RadiusX,RadiusY,RadiusZ,Prec,NbSubdiv);
921 
922 
923  //2) run functions
924  LocMesh=CreateBasicCubicMesh();
925  ProjectMeshEllipsoidAndRefine(&LocMesh,RadiusX,RadiusY,RadiusZ,Prec,NbSubdiv);
926  ExtractSubmesh(&LocMesh,MinX,MaxX,MinY,MaxY,MinZ,MaxZ);
927  WriteMesh(&LocMesh,output_name);
928 
929  return 0;
930 }
Generated on Sun Mar 15 2015 00:11:32 for Nektar++ by   doxygen 1.8.1.2