Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MeshSphericalshell.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 //12.06.2009
13 //
14 //Command line to compile: gcc -o MeshSphericalshell MeshSphericalshell.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 sphere of radius 'R'
272 void ProjectMeshSphere(Mesh * LocMesh, double R){
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 sphere of radius 'R'
302  for (i=0;i<LocMesh->NbVertexes;i++){
303  LocNorm=sqrt(pow(LocMesh->Vertexes[i].X,2.0)+pow(LocMesh->Vertexes[i].Y,2.0)+pow(LocMesh->Vertexes[i].Z,2.0));
304  LocMesh->Vertexes[i].X=LocMesh->Vertexes[i].X*R/LocNorm;
305  LocMesh->Vertexes[i].Y=LocMesh->Vertexes[i].Y*R/LocNorm;
306  LocMesh->Vertexes[i].Z=LocMesh->Vertexes[i].Z*R/LocNorm;
307  }
308 }
309 
310 //order the edges of all element so that this order is clockwise if we see it from outside the shape
311 void MakeClockwiseOrder(Mesh * LocMesh){
312  int i,j,jswap;
313  int Ve1,Ed1,Ve2,Ed2,Ve3,Ed3; //Order: Vertex 1 -> Edge1 -> Vertex 2 -> Edge 2 -> Vertex 3 -> Edge 3 -> Vertex 1 ...
314  int temp;
315  double x1,x2,x3,y1,y2,y3,z1,z2,z3; //coordinates
316  double u1,u2,u3,v1,v2,v3; //vectors
317  double vp1,vp2,vp3; //cross product
318  double tempDbl;
319 
320  //test all elements of the mesh
321  for (i=0;i<LocMesh->NbElements;i++){
322  //order the vertexes and edges
323  Ed1=LocMesh->Elements[i].E1;
324  Ed2=LocMesh->Elements[i].E2;
325  Ed3=LocMesh->Elements[i].E3;
326 
327  if ((LocMesh->Edges[Ed1].V2!=LocMesh->Edges[Ed2].V1)&&(LocMesh->Edges[Ed1].V2!=LocMesh->Edges[Ed2].V2)){ //swap edges
328  temp=Ed2; Ed2=Ed3; Ed3=temp;
329  }
330 
331  if (LocMesh->Edges[Ed1].V2!=LocMesh->Edges[Ed2].V1){ //swap vertexes
332  temp=LocMesh->Edges[Ed2].V1; LocMesh->Edges[Ed2].V1=LocMesh->Edges[Ed2].V2; LocMesh->Edges[Ed2].V2=temp;
333  if (LocMesh->Edges[Ed2].NbSubdiv>1)
334  for (j=0;j<=LocMesh->Edges[Ed2].NbSubdiv/2;j++){
335  jswap=LocMesh->Edges[Ed2].NbSubdiv-j;
336  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;
337  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;
338  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;
339  }
340  }
341 
342  if (LocMesh->Edges[Ed2].V2!=LocMesh->Edges[Ed3].V1){ //swap vertexes
343  temp=LocMesh->Edges[Ed3].V1; LocMesh->Edges[Ed3].V1=LocMesh->Edges[Ed3].V2; LocMesh->Edges[Ed3].V2=temp;
344  if (LocMesh->Edges[Ed3].NbSubdiv>1)
345  for (j=0;j<=LocMesh->Edges[Ed3].NbSubdiv/2;j++){
346  jswap=LocMesh->Edges[Ed3].NbSubdiv-j;
347  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;
348  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;
349  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;
350  }
351  }
352 
353  if (LocMesh->Edges[Ed3].V2!=LocMesh->Edges[Ed1].V1){
354  printf("This is not good!!!\n");
355  }
356 
357  //invert the order if necessary
358  x1=LocMesh->Vertexes[LocMesh->Edges[Ed1].V1].X; y1=LocMesh->Vertexes[LocMesh->Edges[Ed1].V1].Y; z1=LocMesh->Vertexes[LocMesh->Edges[Ed1].V1].Z;
359  x2=LocMesh->Vertexes[LocMesh->Edges[Ed1].V2].X; y2=LocMesh->Vertexes[LocMesh->Edges[Ed1].V2].Y; z2=LocMesh->Vertexes[LocMesh->Edges[Ed1].V2].Z;
360  x3=LocMesh->Vertexes[LocMesh->Edges[Ed3].V1].X; y3=LocMesh->Vertexes[LocMesh->Edges[Ed3].V1].Y; z3=LocMesh->Vertexes[LocMesh->Edges[Ed3].V1].Z;
361 
362  u1=x2-x1; u2=y2-y1; u3=z2-z1;
363  v1=x3-x1; v2=y3-y1; v3=z3-z1;
364 
365  vp1=u2*v3-u3*v2;
366  vp2=u3*v1-u1*v3;
367  vp3=u1*v2-u2*v1;
368 
369  //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);
370  //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);
371 
372  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
373  //swap edges
374  temp=Ed3; Ed3=Ed2; Ed2=temp;
375  //swap vertexes
376  temp=LocMesh->Edges[Ed1].V1; LocMesh->Edges[Ed1].V1=LocMesh->Edges[Ed1].V2; LocMesh->Edges[Ed1].V2=temp;
377  if (LocMesh->Edges[Ed1].NbSubdiv>1)
378  for (j=0;j<=LocMesh->Edges[Ed1].NbSubdiv/2;j++){
379  jswap=LocMesh->Edges[Ed1].NbSubdiv-j;
380  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;
381  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;
382  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;
383  }
384  }
385  //new order of the elements
386  LocMesh->Elements[i].E1=Ed1;
387  LocMesh->Elements[i].E2=Ed2;
388  LocMesh->Elements[i].E3=Ed3;
389  }
390 
391 }
392 
393 //subdivide into 2 edges all edges in the mesh 'LocMesh'
394 void RefineMesh(Mesh * LocMesh){
395  int i,j,k;
396  int * order;
397  double SqLengthI,SqLengthJ;
398  int OrigNbEdges;
399 
400  order=(int*)malloc(LocMesh->NbEdges*sizeof(int));
401 
402  for (i=0;i<LocMesh->NbEdges;i++) order[i]=i;
403 
404  for (i=0;i<LocMesh->NbEdges-1;i++)
405  for (j=i+1;j<LocMesh->NbEdges;j++){
406  SqLengthI=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].X-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].X,2.0);
407  SqLengthI+=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].Y-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].Y,2.0);
408  SqLengthI+=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].Z-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].Z,2.0);
409  SqLengthJ=pow(LocMesh->Vertexes[LocMesh->Edges[order[j]].V1].X-LocMesh->Vertexes[LocMesh->Edges[order[j]].V2].X,2.0);
410  SqLengthJ+=pow(LocMesh->Vertexes[LocMesh->Edges[order[j]].V1].Y-LocMesh->Vertexes[LocMesh->Edges[order[j]].V2].Y,2.0);
411  SqLengthJ+=pow(LocMesh->Vertexes[LocMesh->Edges[order[j]].V1].Z-LocMesh->Vertexes[LocMesh->Edges[order[j]].V2].Z,2.0);
412 
413  if (SqLengthI<SqLengthJ){
414  k=order[i];
415  order[i]=order[j];
416  order[j]=k;
417  }
418  }
419 
420  OrigNbEdges=LocMesh->NbEdges;
421  for (i=0;i<OrigNbEdges;i++) CutEdge(LocMesh,order[i]);
422 }
423 
424 
425 //Project de coordinates of a mesh on a sphere of radius 'R' and refine the mesh until
426 //the longest edge has a length smaller than 'MaxEdgeLength'.
427 //If NbSubdiv>1, each edge is curved into 'NbSubdiv' parts in the end
428 void ProjectMeshSphereAndRefine(Mesh * LocMesh, double R,double MaxEdgeLength,int NbSubdiv){
429  int i,j,k;
430  int * order;
431  double SqLengthI,SqLengthJ;
432  int OrigNbEdges;
433  int OK;
434  double LocNorm;
435  double dj,dNbSubdiv;
436  int jm1;
437  double EdgeV1_X,EdgeV1_Y,EdgeV1_Z,EdgeV2_X,EdgeV2_Y,EdgeV2_Z;
438 
439  //initialize
440  order=(int*)malloc(LocMesh->NbEdges*sizeof(int));
441  for (i=0;i<LocMesh->NbEdges;i++) order[i]=i;
442  OK=1;
443 
444  //mesh projection on a sphere
445  ProjectMeshSphere(LocMesh,R);
446 
447  //Big loop to subdivide the mesh until each edge has a length smaller than 'MaxEdgeLength'
448  while (OK==1){
449  OK=0;
450 
451  order=(int*)realloc(order,LocMesh->NbEdges*sizeof(int));
452  for (i=0;i<LocMesh->NbEdges;i++) order[i]=i;
453 
454  //order the edges as a functino of their length
455  for (i=0;i<LocMesh->NbEdges-1;i++) for (j=i+1;j<LocMesh->NbEdges;j++){
456  SqLengthI=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].X-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].X,2.0);
457  SqLengthI+=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].Y-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].Y,2.0);
458  SqLengthI+=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].Z-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].Z,2.0);
459  SqLengthJ=pow(LocMesh->Vertexes[LocMesh->Edges[order[j]].V1].X-LocMesh->Vertexes[LocMesh->Edges[order[j]].V2].X,2.0);
460  SqLengthJ+=pow(LocMesh->Vertexes[LocMesh->Edges[order[j]].V1].Y-LocMesh->Vertexes[LocMesh->Edges[order[j]].V2].Y,2.0);
461  SqLengthJ+=pow(LocMesh->Vertexes[LocMesh->Edges[order[j]].V1].Z-LocMesh->Vertexes[LocMesh->Edges[order[j]].V2].Z,2.0);
462 
463  if (SqLengthI<SqLengthJ){ k=order[i]; order[i]=order[j]; order[j]=k;}
464  }
465 
466  //subdivide the edges if their length is larger than 'MaxEdgeLength'
467  OrigNbEdges=LocMesh->NbEdges;
468  for (i=0;i<OrigNbEdges;i++){
469  SqLengthI=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].X-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].X,2.0);
470  SqLengthI+=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].Y-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].Y,2.0);
471  SqLengthI+=pow(LocMesh->Vertexes[LocMesh->Edges[order[i]].V1].Z-LocMesh->Vertexes[LocMesh->Edges[order[i]].V2].Z,2.0);
472 
473  if (sqrt(SqLengthI)>MaxEdgeLength){
474  //subdvision
475  OK=1;
476  CutEdge(LocMesh,order[i]); //rem: 'LocMesh' is changed here. In particular LocMesh->NbEdges is larger.
477  }
478  }
479 
480  //projection of the new edges on a sphere
481  ProjectMeshSphere(LocMesh,R);
482  }
483 
484 
485 
486  //subdivide the edges end project them on the sphere if required
487  if (NbSubdiv>1){
488  for (i=0;i<LocMesh->NbEdges;i++){
489  //init
490  LocMesh->Edges[i].NbSubdiv=NbSubdiv;
491  LocMesh->Edges[i].SubiVertex=(Vertex *)malloc((LocMesh->Edges[i].NbSubdiv+1)*sizeof(Vertex));
492 
493  //coordinates of the subdivised edge
494  EdgeV1_X=LocMesh->Vertexes[LocMesh->Edges[i].V1].X;
495  EdgeV1_Y=LocMesh->Vertexes[LocMesh->Edges[i].V1].Y;
496  EdgeV1_Z=LocMesh->Vertexes[LocMesh->Edges[i].V1].Z;
497  EdgeV2_X=LocMesh->Vertexes[LocMesh->Edges[i].V2].X;
498  EdgeV2_Y=LocMesh->Vertexes[LocMesh->Edges[i].V2].Y;
499  EdgeV2_Z=LocMesh->Vertexes[LocMesh->Edges[i].V2].Z;
500  dNbSubdiv=(double)LocMesh->Edges[i].NbSubdiv;
501 
502  for (j=0;j<LocMesh->Edges[i].NbSubdiv+1;j++){
503  dj=(double)j;
504 
505  //interpolation
506  LocMesh->Edges[i].SubiVertex[j].X=EdgeV1_X*(dNbSubdiv-dj)/dNbSubdiv+EdgeV2_X*dj/dNbSubdiv;
507  LocMesh->Edges[i].SubiVertex[j].Y=EdgeV1_Y*(dNbSubdiv-dj)/dNbSubdiv+EdgeV2_Y*dj/dNbSubdiv;
508  LocMesh->Edges[i].SubiVertex[j].Z=EdgeV1_Z*(dNbSubdiv-dj)/dNbSubdiv+EdgeV2_Z*dj/dNbSubdiv;
509 
510  //projection on a sphere
511  LocNorm=sqrt(pow(LocMesh->Edges[i].SubiVertex[j].X,2.0)+pow(LocMesh->Edges[i].SubiVertex[j].Y,2.0)+pow(LocMesh->Edges[i].SubiVertex[j].Z,2.0));
512  LocMesh->Edges[i].SubiVertex[j].X=LocMesh->Edges[i].SubiVertex[j].X*R/LocNorm;
513  LocMesh->Edges[i].SubiVertex[j].Y=LocMesh->Edges[i].SubiVertex[j].Y*R/LocNorm;
514  LocMesh->Edges[i].SubiVertex[j].Z=LocMesh->Edges[i].SubiVertex[j].Z*R/LocNorm;
515  }
516  }
517  }
518 }
519 
520 
521 
522 
523 //write the geometry of mesh contained in a Mesh structure in a XML-Nektar file
524 void WriteMesh(Mesh * LocMesh,char FileName[256]){
525  FILE * XmlMeshGeomFile;
526  int i,j,NbCurvedEdge;
527 
528  //to have a proper order of the edges within the elements
529  MakeClockwiseOrder(LocMesh);
530 
531  //open file
532  XmlMeshGeomFile = fopen(FileName,"w");
533 
534  //write header
535  fprintf(XmlMeshGeomFile,"<?xml version=\"1.0\" encoding=\"utf-8\"?>\n");
536  fprintf(XmlMeshGeomFile,"\n");
537  fprintf(XmlMeshGeomFile,"<NEKTAR>\n");
538  fprintf(XmlMeshGeomFile,"<!-- Embed a 2-dimensional object in a 2-dimensional space -->\n");
539  fprintf(XmlMeshGeomFile,"<!-- DIM <= SPACE -->\n");
540  fprintf(XmlMeshGeomFile,"<!-- This provides a method of optimizing code for a 1-D curve embedded in 3-space. -->\n");
541  fprintf(XmlMeshGeomFile,"<GEOMETRY DIM=\"2\" SPACE=\"3\">\n");
542  fprintf(XmlMeshGeomFile,"\n");
543 
544  //write vertex
545  fprintf(XmlMeshGeomFile," <VERTEX>\n");
546  fprintf(XmlMeshGeomFile," <!-- Always must have four values per entry. -->\n");
547  for (i=0;i<LocMesh->NbVertexes;i++)
548  fprintf(XmlMeshGeomFile," <V ID=\"%d\"> %f %f %f </V>\n",i,LocMesh->Vertexes[i].X,LocMesh->Vertexes[i].Y,LocMesh->Vertexes[i].Z);
549  fprintf(XmlMeshGeomFile," </VERTEX>\n");
550  fprintf(XmlMeshGeomFile," \n");
551 
552  //write edges
553  fprintf(XmlMeshGeomFile," <EDGE>\n");
554  fprintf(XmlMeshGeomFile," <!--Edges are vertex pairs -->\n");
555  for (i=0;i<LocMesh->NbEdges;i++)
556  fprintf(XmlMeshGeomFile," <E ID=\"%d\"> %d %d </E>\n",i,LocMesh->Edges[i].V1,LocMesh->Edges[i].V2);
557  fprintf(XmlMeshGeomFile," </EDGE>\n");
558  fprintf(XmlMeshGeomFile," \n");
559 
560  //write Element
561  fprintf(XmlMeshGeomFile," <ELEMENT>\n");
562  for (i=0;i<LocMesh->NbElements;i++)
563  fprintf(XmlMeshGeomFile," <T ID=\"%d\"> %d %d %d </T>\n",i,LocMesh->Elements[i].E1,LocMesh->Elements[i].E2,LocMesh->Elements[i].E3);
564  fprintf(XmlMeshGeomFile," </ELEMENT>\n");
565  fprintf(XmlMeshGeomFile," \n");
566 
567 
568 
569  //write curved elements
570  fprintf(XmlMeshGeomFile," <CURVED>\n");
571  NbCurvedEdge=0;
572  for (i=0;i<LocMesh->NbEdges;i++) if (LocMesh->Edges[i].NbSubdiv>1){
573  fprintf(XmlMeshGeomFile," <E ID=\"%d\" EDGEID=\"%d\" TYPE=\"PolyEvenlySpaced\" NUMPOINTS=\"%d\">\n",NbCurvedEdge,i,LocMesh->Edges[i].NbSubdiv+1);
574  for (j=0;j<LocMesh->Edges[i].NbSubdiv+1;j++){
575  fprintf(XmlMeshGeomFile," %lf %lf %lf\n",LocMesh->Edges[i].SubiVertex[j].X,LocMesh->Edges[i].SubiVertex[j].Y,LocMesh->Edges[i].SubiVertex[j].Z);
576  }
577  fprintf(XmlMeshGeomFile," </E>\n\n");
578  NbCurvedEdge++;
579  }
580  fprintf(XmlMeshGeomFile," </CURVED>\n");
581 
582 
583  //write the end (could be more evolved)
584  fprintf(XmlMeshGeomFile,"<!-- V - vertex, E - edge, F - face, L - element -->\n");
585  fprintf(XmlMeshGeomFile," <COMPOSITE>\n");
586  fprintf(XmlMeshGeomFile," <C ID=\"0\"> T[0-%d]\n",LocMesh->NbElements-1);
587  fprintf(XmlMeshGeomFile," </C>\n");
588  fprintf(XmlMeshGeomFile," <C ID=\"1\"> E[2,3,4,5,7,8]\n");
589  fprintf(XmlMeshGeomFile," </C>\n");
590  fprintf(XmlMeshGeomFile," </COMPOSITE>\n");
591  fprintf(XmlMeshGeomFile," \n");
592  fprintf(XmlMeshGeomFile," <DOMAIN> C[0] </DOMAIN>\n");
593  fprintf(XmlMeshGeomFile," \n");
594 
595  //EOF
596  fprintf(XmlMeshGeomFile,"</GEOMETRY>\n");
597 
598  //close file
599  fclose(XmlMeshGeomFile);
600 }
601 
602 
603 
604 
605 
606 // ---------------------------------------------------------------------------------------
607 // main function
608 // ---------------------------------------------------------------------------------------
609 
610 
611 void usage()
612 {
613  printf("Usage: makeSphereNektar <options>\n");
614  printf("Where <options> are one or more of the following:\n");
615  printf("\t<-NameFile n> Name of the output file containg the sphere (default=\"toto.xml\")\n");
616  printf("\t<-Radius n> Radius of the sphere (default=1.)\n");
617  printf("\t<-MaxLength n> Maximum length allowed for an edge (default=0.5)\n");
618  printf("\t<-NbSubdiv n> Number of subdivisions of each edge to make it curved (default=1)\n");
619  exit(1);
620 }
621 
622 int main(int argc, char **argv)
623 {
624  int ok,i;
625  char *output_name = "toto.xml";
626  double Radius=1.;
627  double Prec=0.5;
628  int NbSubdiv=1;
629  Mesh LocMesh;
630 
631 
632  //1) Parse parameters
633  while (argc > 1) {
634  ok = 0;
635  if ((ok == 0) && (strcmp(argv[1], "help") == 0)) {
636  usage();
637  }
638  if ((ok == 0) && (strcmp(argv[1], "-NameFile") == 0)) {
639  argc--;
640  argv++;
641  output_name = argv[1];
642  argc--;
643  argv++;
644  ok = 1;
645  }
646  if ((ok == 0) && (strcmp(argv[1], "-Radius") == 0)) {
647  argc--;
648  argv++;
649  Radius = atof(argv[1]);
650  if (Radius<0) Radius=1.;
651  argc--;
652  argv++;
653  ok = 1;
654  }
655  if ((ok == 0) && (strcmp(argv[1], "-MaxLength") == 0)) {
656  argc--;
657  argv++;
658  Prec = atof(argv[1]);
659  if ((Prec<Radius*0.001)) Prec=Radius*0.001;
660  argc--;
661  argv++;
662  ok = 1;
663  }
664  if ((ok == 0) && (strcmp(argv[1], "-NbSubdiv") == 0)) {
665  argc--;
666  argv++;
667  NbSubdiv = atoi(argv[1]);
668  if ((NbSubdiv<1)) NbSubdiv=1;
669  argc--;
670  argv++;
671  ok = 1;
672  }
673  if (ok == 0) {
674  printf("Unknown option: \n");
675  usage();
676  }
677  }
678 
679  printf("Output file: %s / Radius: %lf / Precision: %lf / Subdivisions: %d\n",output_name,Radius,Prec,NbSubdiv);
680 
681 
682  //2) run functions
683  LocMesh=CreateBasicCubicMesh();
684  ProjectMeshSphereAndRefine(&LocMesh,Radius,Prec,NbSubdiv);
685  WriteMesh(&LocMesh,output_name);
686 
687  return 0;
688 }
Generated on Sun Mar 15 2015 00:11:32 for Nektar++ by   doxygen 1.8.1.2