/[escript]/branches/domexper/dudley/src/Mesh_tet4.c
ViewVC logotype

Annotation of /branches/domexper/dudley/src/Mesh_tet4.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3221 - (hide annotations)
Wed Sep 29 01:00:21 2010 UTC (8 years, 7 months ago) by jfenwick
File MIME type: text/plain
File size: 20886 byte(s)
Comment stripping

1 jgs 150
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13 ksteube 1312
14 ksteube 1811
15 jgs 82 /**************************************************************/
16    
17 jfenwick 3086 /* Dudley: generates rectangular meshes */
18 jgs 82
19     /* Generates a numElements[0] x numElements[1] x numElements[2] mesh with first order elements (Hex8) in the brick */
20     /* [0,Length[0]] x [0,Length[1]] x [0,Length[2]]. order is the desired accuracy of the */
21     /* integration scheme. */
22    
23     /**************************************************************/
24    
25 jfenwick 3114 #include "TriangularMesh.h"
26 jgs 82
27 jfenwick 3114
28 jfenwick 3118 /* Be careful reading this function. The X? and NStride? are 1,2,3 but the loop vars are 0,1,2 */
29 jfenwick 3114 Dudley_Mesh* Dudley_TriangularMesh_Tet4(dim_t* numElements,
30 ksteube 1312 double* Length,
31 jfenwick 3114 index_t order,
32 ksteube 1312 index_t reduced_order,
33     bool_t optimize)
34 bcumming 751 {
35 ksteube 1312 #define N_PER_E 1
36     #define DIM 3
37 gross 2748 dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,Nstride0=0, Nstride1=0,Nstride2=0, local_NE0, local_NE1, local_NE2, local_N0=0, local_N1=0, local_N2=0;
38     dim_t totalNECount,faceNECount,NDOF0=0,NDOF1=0,NDOF2=0,NFaceElements=0, NN;
39     index_t node0, myRank, e_offset2, e_offset1, e_offset0=0, offset1=0, offset2=0, offset0=0, global_i0, global_i1, global_i2;
40 jfenwick 3086 Dudley_Mesh* out;
41 ksteube 1312 Paso_MPIInfo *mpi_info = NULL;
42     char name[50];
43 jfenwick 3086 #ifdef Dudley_TRACE
44     double time0=Dudley_timer();
45 jfenwick 1981 #endif
46 bcumming 751
47 jfenwick 3114 const int LEFTTAG=1; /* boundary x1=0 */
48     const int RIGHTTAG=2; /* boundary x1=1 */
49     const int BOTTOMTAG=100; /* boundary x3=1 */
50     const int TOPTAG=200; /* boundary x3=0 */
51     const int FRONTTAG=10; /* boundary x2=0 */
52     const int BACKTAG=20; /* boundary x2=1 */
53    
54    
55    
56    
57 ksteube 1312 /* get MPI information */
58     mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
59 jfenwick 3086 if (! Dudley_noError()) {
60 ksteube 1312 return NULL;
61 bcumming 751 }
62 ksteube 1312 myRank=mpi_info->rank;
63 bcumming 751
64 ksteube 1312 /* set up the global dimensions of the mesh */
65 bcumming 751
66 jgs 82 NE0=MAX(1,numElements[0]);
67     NE1=MAX(1,numElements[1]);
68     NE2=MAX(1,numElements[2]);
69 ksteube 1312 N0=N_PER_E*NE0+1;
70     N1=N_PER_E*NE1+1;
71     N2=N_PER_E*NE2+1;
72 jgs 82
73 ksteube 1312 /* allocate mesh: */
74 jfenwick 3114 sprintf(name,"Triangular %d x %d x %d (x 5) mesh",N0,N1,N2);
75 jfenwick 3086 out=Dudley_Mesh_alloc(name,DIM, mpi_info);
76     if (! Dudley_noError()) {
77 ksteube 1312 Paso_MPIInfo_free( mpi_info );
78     return NULL;
79 jgs 153 }
80 jfenwick 3086 if ( Dudley_noError()) {
81 gross 2748
82 jfenwick 3216 Dudley_Mesh_setPoints(out,Dudley_ElementFile_alloc(Point1, mpi_info));
83     Dudley_Mesh_setFaceElements(out,Dudley_ElementFile_alloc(Tri3, mpi_info));
84     Dudley_Mesh_setElements(out,Dudley_ElementFile_alloc(Tet4, mpi_info));
85 ksteube 1312
86 gross 2748 /* work out the largest dimension */
87     if (N2==MAX3(N0,N1,N2)) {
88     Nstride0=1;
89     Nstride1=N0;
90     Nstride2=N0*N1;
91     local_NE0=NE0;
92     e_offset0=0;
93     local_NE1=NE1;
94     e_offset1=0;
95     Paso_MPIInfo_Split(mpi_info,NE2,&local_NE2,&e_offset2);
96     } else if (N1==MAX3(N0,N1,N2)) {
97     Nstride0=N2;
98     Nstride1=N0*N2;
99     Nstride2=1;
100     local_NE0=NE0;
101     e_offset0=0;
102     Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
103     local_NE2=NE2;
104     e_offset2=0;
105     } else {
106     Nstride0=N1*N2;
107     Nstride1=1;
108     Nstride2=N1;
109     Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
110     local_NE1=NE1;
111     e_offset1=0;
112     local_NE2=NE2;
113     e_offset2=0;
114     }
115     offset0=e_offset0*N_PER_E;
116     offset1=e_offset1*N_PER_E;
117     offset2=e_offset2*N_PER_E;
118     local_N0=local_NE0>0 ? local_NE0*N_PER_E+1 : 0;
119     local_N1=local_NE1>0 ? local_NE1*N_PER_E+1 : 0;
120     local_N2=local_NE2>0 ? local_NE2*N_PER_E+1 : 0;
121 bcumming 751
122 gross 2748 /* get the number of surface elements */
123 jgs 82
124 gross 2748 NFaceElements=0;
125 jfenwick 3114 if (local_NE2>0) {
126 gross 2748 NDOF2=N2;
127 jfenwick 3114 if (offset2==0) NFaceElements+=2*local_NE1*local_NE0; // each face is split
128     if (local_NE2+e_offset2 == NE2) NFaceElements+=2*local_NE1*local_NE0;
129 gross 2748 } else {
130     NDOF2=N2-1;
131     }
132 ksteube 1312
133 jfenwick 3114 if (local_NE0>0) {
134 gross 2748 NDOF0=N0;
135 jfenwick 3114 if (e_offset0 == 0) NFaceElements+=2*local_NE1*local_NE2;
136     if (local_NE0+e_offset0 == NE0) NFaceElements+=2*local_NE1*local_NE2;
137 gross 2748 } else {
138     NDOF0=N0-1;
139     }
140 jfenwick 3114
141     if (local_NE1>0) {
142     NDOF1=N1;
143     if (e_offset1 == 0) NFaceElements+=2*local_NE0*local_NE2;
144     if (local_NE1+e_offset1 == NE1) NFaceElements+=2*local_NE0*local_NE2;
145     } else {
146     NDOF1=N1-1;
147     }
148 jgs 82 }
149 gross 2748
150    
151 ksteube 1312 /* allocate tables: */
152 jfenwick 3086 if (Dudley_noError()) {
153 jgs 82
154 jfenwick 3086 Dudley_NodeFile_allocTable(out->Nodes,local_N0*local_N1*local_N2);
155 jfenwick 3114 /* we split the rectangular prism this code used to produce into 5 tetrahedrons */
156     Dudley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1*local_NE2*5);
157     /* each border face will be split in half */
158 jfenwick 3086 Dudley_ElementFile_allocTable(out->FaceElements,NFaceElements);
159 gross 2748 }
160    
161 jfenwick 3086 if (Dudley_noError()) {
162 ksteube 1312 /* create nodes */
163 jgs 82
164 ksteube 1312 #pragma omp parallel for private(i0,i1,i2,k,global_i0,global_i1,global_i2)
165     for (i2=0;i2<local_N2;i2++) {
166     for (i1=0;i1<local_N1;i1++) {
167     for (i0=0;i0<local_N0;i0++) {
168     k=i0+local_N0*i1+local_N0*local_N1*i2;
169     global_i0=i0+offset0;
170     global_i1=i1+offset1;
171     global_i2=i2+offset2;
172     out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
173     out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
174     out->Nodes->Coordinates[INDEX2(2,k,DIM)]=DBLE(global_i2)/DBLE(N2-1)*Length[2];
175     out->Nodes->Id[k]=Nstride0*global_i0+Nstride1*global_i1+Nstride2*global_i2;
176     out->Nodes->Tag[k]=0;
177     out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
178     +Nstride1*(global_i1%NDOF1)
179     +Nstride2*(global_i2%NDOF2);
180 jgs 82 }
181     }
182     }
183 ksteube 1312 /* set the elements: */
184 jfenwick 3114
185 jfenwick 3122 int global_adjustment=(offset0+offset1+offset2)%2; // If we are not the only rank we may need to shift our pattern to match neighbours
186 jfenwick 3114
187 ksteube 1312 NN=out->Elements->numNodes;
188     #pragma omp parallel for private(i0,i1,i2,k,node0)
189     for (i2=0;i2<local_NE2;i2++) {
190     for (i1=0;i1<local_NE1;i1++) {
191     for (i0=0;i0<local_NE0;i0++) {
192    
193 jfenwick 3114 k=5*(i0+local_NE0*i1+local_NE0*local_NE1*i2);
194 ksteube 1312 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
195 jgs 82
196 jfenwick 3114 index_t res=5*((i0+e_offset0)+NE0*(i1+e_offset1)+NE0*NE1*(i2+e_offset2));
197     for (int j=0;j<5;++j)
198     {
199     out->Elements->Id[k+j]=res+j;
200     out->Elements->Tag[k+j]=0;
201     out->Elements->Owner[k+j]=myRank;
202     }
203    
204     index_t v[8];
205     /* in non-rotated orientation the points are numbered as follows:
206 jfenwick 3117 The bottom face (anticlockwise= 0,1,3,2), top face (anticlockwise 4,5,7,6)*/
207 jfenwick 3114 if ((global_adjustment+i0+i1+i2)%2==0)
208     {
209     v[0]=node0;
210     v[1]=node0+Nstride0;
211     v[2]=node0+Nstride1;
212     v[3]=node0+Nstride1+Nstride0;
213     v[4]=node0+Nstride2;
214     v[5]=node0+Nstride0+Nstride2;
215     v[6]=node0+Nstride1+Nstride2;
216     v[7]=node0+Nstride2+Nstride1+Nstride0;
217     }
218     else
219     {
220     // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
221    
222 jfenwick 3117 v[0]=node0+Nstride1; // node 0 ends up in position 2
223     v[2]=node0+Nstride1+Nstride2; // node 2 ends up in position 6
224     v[6]=node0+Nstride2; // node 6 ends up in position 4
225     v[4]=node0; // node 4 ends up in position 0
226 jfenwick 3114
227 jfenwick 3117 v[1]=node0+Nstride0+Nstride1; // node 1 -> pos 3
228     v[3]=node0+Nstride2+Nstride1+Nstride0; // node 3-> pos 7
229     v[7]=node0+Nstride0+Nstride2; // node 7 -> pos 5
230     v[5]=node0+Nstride0; // node 5 -> pos 1
231 jfenwick 3114 }
232 jfenwick 3117
233 jfenwick 3114 // elements nodes are numbered: centre, x, y, z
234    
235     out->Elements->Nodes[INDEX2(0,k,NN)] =v[4];
236     out->Elements->Nodes[INDEX2(1,k,NN)] =v[5];
237     out->Elements->Nodes[INDEX2(2,k,NN)] =v[6];
238     out->Elements->Nodes[INDEX2(3,k,NN)] =v[0];
239    
240     out->Elements->Nodes[INDEX2(0,k+1,NN)] =v[7];
241     out->Elements->Nodes[INDEX2(1,k+1,NN)] =v[6];
242     out->Elements->Nodes[INDEX2(2,k+1,NN)] =v[5];
243     out->Elements->Nodes[INDEX2(3,k+1,NN)] =v[3];
244    
245     out->Elements->Nodes[INDEX2(0,k+2,NN)] =v[2];
246     out->Elements->Nodes[INDEX2(1,k+2,NN)] =v[3];
247     out->Elements->Nodes[INDEX2(2,k+2,NN)] =v[0];
248     out->Elements->Nodes[INDEX2(3,k+2,NN)] =v[6];
249    
250 jfenwick 3117
251 jfenwick 3114 out->Elements->Nodes[INDEX2(0,k+3,NN)] =v[1];
252     out->Elements->Nodes[INDEX2(1,k+3,NN)] =v[0];
253     out->Elements->Nodes[INDEX2(2,k+3,NN)] =v[3];
254     out->Elements->Nodes[INDEX2(3,k+3,NN)] =v[5];
255    
256     // I can't work out where the center is for this one
257     out->Elements->Nodes[INDEX2(0,k+4,NN)] =v[5];
258     out->Elements->Nodes[INDEX2(1,k+4,NN)] =v[0];
259     out->Elements->Nodes[INDEX2(2,k+4,NN)] =v[6];
260     out->Elements->Nodes[INDEX2(3,k+4,NN)] =v[3];
261 jgs 82 }
262 jfenwick 3114
263 jgs 82 }
264 jfenwick 3114 } /* end for */
265 ksteube 1312 /* face elements */
266     NN=out->FaceElements->numNodes;
267 jfenwick 3114 totalNECount=5*NE0*NE1*NE2;
268 ksteube 1312 faceNECount=0;
269     /* these are the quadrilateral elements on boundary 1 (x3=0): */
270 jfenwick 3114 if (local_NE2>0)
271     {
272 ksteube 1312 /* ** elements on boundary 100 (x3=0): */
273 jfenwick 3114 if (e_offset2==0)
274     {
275 ksteube 1312 #pragma omp parallel for private(i0,i1,k,node0)
276 jfenwick 3114 for (i1=0;i1<local_NE1;i1++)
277     {
278     for (i0=0;i0<local_NE0;i0++)
279     {
280     k=2*(i0+local_NE0*i1)+faceNECount;
281 ksteube 1312 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
282 jfenwick 3114 index_t res=2*((i0+e_offset0)+NE0*(i1+e_offset1))+totalNECount;
283     out->FaceElements->Id[k]=res;
284 jfenwick 3118 out->FaceElements->Tag[k]=BOTTOMTAG;
285 ksteube 1312 out->FaceElements->Owner[k]=myRank;
286 jfenwick 3114 out->FaceElements->Id[k+1]=res+1;
287 jfenwick 3118 out->FaceElements->Tag[k+1]=BOTTOMTAG;
288 jfenwick 3114 out->FaceElements->Owner[k+1]=myRank;
289    
290 jfenwick 3118 index_t n0=node0;
291     index_t n1=node0+Nstride0;
292     index_t n2=node0+Nstride1;
293     index_t n3=node0+Nstride0+Nstride1;
294 jfenwick 3117
295 jfenwick 3114 if ((global_adjustment+i0+i1)%2==0)
296     {
297 jfenwick 3118 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0;
298     out->FaceElements->Nodes[INDEX2(1,k,NN)]=n3;
299     out->FaceElements->Nodes[INDEX2(2,k,NN)]=n1;
300 jfenwick 3114
301 jfenwick 3118 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n0;
302     out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n2;
303     out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n3;
304 jfenwick 3114
305     }
306     else
307     {
308 jfenwick 3118 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0;
309     out->FaceElements->Nodes[INDEX2(1,k,NN)]=n2;
310     out->FaceElements->Nodes[INDEX2(2,k,NN)]=n1;
311 jfenwick 3114
312 jfenwick 3118 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n1;
313     out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n2;
314     out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n3;
315 jfenwick 3114
316 jfenwick 3117
317    
318 jfenwick 3114 }
319 ksteube 1312 }
320     }
321 jfenwick 3114 faceNECount+=2*local_NE1*local_NE0;
322 jgs 82 }
323 jfenwick 3114 totalNECount+=2*NE1*NE0;
324 jfenwick 3221 /* ** elements on boundary 200 (x3=1) - Top*/
325 ksteube 1312 if (local_NE2+e_offset2 == NE2) {
326     #pragma omp parallel for private(i0,i1,k,node0)
327     for (i1=0;i1<local_NE1;i1++) {
328     for (i0=0;i0<local_NE0;i0++) {
329    
330 jfenwick 3114 k=2*(i0+local_NE0*i1)+faceNECount;
331 ksteube 1312 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(NE2-1);
332 jfenwick 3114
333     index_t res=2*((i0+e_offset0)+NE0*(i1+e_offset1))+totalNECount;
334     out->FaceElements->Id[k]=res;
335 jfenwick 3118 out->FaceElements->Tag[k]=TOPTAG;
336 ksteube 1312 out->FaceElements->Owner[k]=myRank;
337 jfenwick 3114 out->FaceElements->Id[k+1]=res+1;
338 jfenwick 3118 out->FaceElements->Tag[k+1]=TOPTAG;
339 jfenwick 3114 out->FaceElements->Owner[k+1]=myRank;
340    
341 jfenwick 3118 index_t n4=node0+Nstride2;
342     index_t n5=node0+Nstride0+Nstride2;
343     index_t n6=node0+Nstride1+Nstride2;
344     index_t n7=node0+Nstride1+Nstride0+Nstride2;
345 jfenwick 3117
346 jfenwick 3114 if ((global_adjustment+i0+i1+local_NE2-1)%2==0)
347     {
348 jfenwick 3118 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n4;
349     out->FaceElements->Nodes[INDEX2(1,k,NN)]=n5;
350     out->FaceElements->Nodes[INDEX2(2,k,NN)]=n6;
351 jfenwick 3114
352 jfenwick 3118 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n5;
353     out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n7;
354     out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n6;
355 jfenwick 3114 }
356     else
357     {
358    
359 jfenwick 3118 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n4;
360     out->FaceElements->Nodes[INDEX2(1,k,NN)]=n5;
361     out->FaceElements->Nodes[INDEX2(2,k,NN)]=n7;
362 jfenwick 3114
363 jfenwick 3118 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n4;
364     out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n7;
365     out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n6;
366 jfenwick 3114 }
367 ksteube 1312 }
368     }
369 jfenwick 3114 faceNECount+=2*local_NE1*local_NE0;
370 jgs 82 }
371 jfenwick 3114 totalNECount+=2*NE1*NE0;
372 jgs 82 }
373 jfenwick 3114 if (local_NE0>0) {
374 jfenwick 3221 /* ** elements on boundary 001 (x1=0): - Left*/
375 ksteube 1312
376     if (e_offset0 == 0) {
377     #pragma omp parallel for private(i1,i2,k,node0)
378     for (i2=0;i2<local_NE2;i2++) {
379     for (i1=0;i1<local_NE1;i1++) {
380    
381 jfenwick 3114 k=2*(i1+local_NE1*i2)+faceNECount;
382 ksteube 1312 node0=Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
383 jfenwick 3114 index_t res=2*((i1+e_offset1)+NE1*(i2+e_offset2))+totalNECount;
384     out->FaceElements->Id[k]=res;
385     out->FaceElements->Tag[k]=LEFTTAG;
386 ksteube 1312 out->FaceElements->Owner[k]=myRank;
387 jfenwick 3114 out->FaceElements->Id[k+1]=res+1;
388     out->FaceElements->Tag[k+1]=LEFTTAG;
389     out->FaceElements->Owner[k+1]=myRank;
390 ksteube 1312
391 jfenwick 3117 index_t n0=node0;
392     index_t n2=node0+Nstride1;
393     index_t n4=node0+Nstride2;
394     index_t n6=node0+Nstride1+Nstride2;
395    
396 jfenwick 3114 if ((global_adjustment+0+i1+i2)%2==0)
397     {
398 jfenwick 3117 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0;
399     out->FaceElements->Nodes[INDEX2(1,k,NN)]=n4;
400     out->FaceElements->Nodes[INDEX2(2,k,NN)]=n6;
401 jfenwick 3114
402 jfenwick 3117 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n0;
403     out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n6;
404     out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n2;
405 jfenwick 3114 }
406     else
407     {
408     // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
409    
410 jfenwick 3117 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0;
411     out->FaceElements->Nodes[INDEX2(1,k,NN)]=n4;
412     out->FaceElements->Nodes[INDEX2(2,k,NN)]=n2;
413 jfenwick 3114
414 jfenwick 3117 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n4;
415     out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n6;
416     out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n2;
417 jfenwick 3114 }
418 ksteube 1312 }
419     }
420 jfenwick 3114 faceNECount+=2*local_NE1*local_NE2;
421 ksteube 1312 }
422 jfenwick 3114 totalNECount+=2*NE1*NE2;
423 jfenwick 3221 /* ** elements on boundary 002 (x1=1): - Right*/
424 ksteube 1312 if (local_NE0+e_offset0 == NE0) {
425     #pragma omp parallel for private(i1,i2,k,node0)
426     for (i2=0;i2<local_NE2;i2++) {
427     for (i1=0;i1<local_NE1;i1++) {
428 jfenwick 3114 k=2*(i1+local_NE1*i2)+faceNECount;
429    
430 ksteube 1312 node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
431 jfenwick 3114 index_t res=2*((i1+e_offset1)+NE1*(i2+e_offset2))+totalNECount;
432     out->FaceElements->Id[k]=res;
433     out->FaceElements->Tag[k]=RIGHTTAG;
434 ksteube 1312 out->FaceElements->Owner[k]=myRank;
435 jfenwick 3114 out->FaceElements->Id[k+1]=res+1;
436     out->FaceElements->Tag[k+1]=RIGHTTAG;
437     out->FaceElements->Owner[k+1]=myRank;
438    
439 jfenwick 3117 index_t n1=node0+Nstride0;
440     index_t n3=node0+Nstride0+Nstride1;
441     index_t n5=node0+Nstride0+Nstride2;
442     index_t n7=node0+Nstride0+Nstride1+Nstride2;
443 jfenwick 3114 if ((global_adjustment+local_NE0-1+i1+i2)%2==0)
444     {
445 jfenwick 3117 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n1;
446     out->FaceElements->Nodes[INDEX2(1,k,NN)]=n3;
447     out->FaceElements->Nodes[INDEX2(2,k,NN)]=n5;
448 jfenwick 3114
449 jfenwick 3117 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n3;
450     out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n7;
451     out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n5;
452 jfenwick 3114 }
453     else
454     {
455     // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
456    
457 jfenwick 3117 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n1;
458     out->FaceElements->Nodes[INDEX2(1,k,NN)]=n7;
459     out->FaceElements->Nodes[INDEX2(2,k,NN)]=n5;
460 jfenwick 3114
461 jfenwick 3117 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n1;
462     out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n3;
463     out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n7;
464 jfenwick 3114 }
465 ksteube 1312 }
466     }
467 jfenwick 3114 faceNECount+=2*local_NE1*local_NE2;
468 jgs 82 }
469 jfenwick 3114 totalNECount+=2*NE1*NE2;
470 jgs 82 }
471 jfenwick 3114 if (local_NE1>0) {
472 jfenwick 3221 /* ** elements on boundary 010 (x2=0): -Front*/
473 ksteube 1312 if (e_offset1 == 0) {
474     #pragma omp parallel for private(i0,i2,k,node0)
475     for (i2=0;i2<local_NE2;i2++) {
476     for (i0=0;i0<local_NE0;i0++) {
477 jfenwick 3114 k=2*(i0+local_NE0*i2)+faceNECount;
478 ksteube 1312 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride2*N_PER_E*(i2+e_offset2);
479 jfenwick 3114 index_t res=2*((i2+e_offset2)+NE2*(e_offset0+i0))+totalNECount;
480     out->FaceElements->Id[k]=res;
481     out->FaceElements->Tag[k]=FRONTTAG;
482 ksteube 1312 out->FaceElements->Owner[k]=myRank;
483 jfenwick 3114 out->FaceElements->Id[k+1]=res+1;
484     out->FaceElements->Tag[k+1]=FRONTTAG;
485     out->FaceElements->Owner[k+1]=myRank;
486    
487 jfenwick 3117 index_t n0=node0;
488     index_t n1=node0+Nstride0;
489     index_t n4=node0+Nstride2;
490     index_t n5=node0+Nstride0+Nstride2;
491    
492 jfenwick 3114 if ((global_adjustment+i0+0+i2)%2==0)
493     {
494 jfenwick 3117 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0;
495     out->FaceElements->Nodes[INDEX2(1,k,NN)]=n1;
496     out->FaceElements->Nodes[INDEX2(2,k,NN)]=n5;
497 jfenwick 3114
498 jfenwick 3117 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n0;
499     out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n5;
500     out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n4;
501 jfenwick 3114
502    
503     }
504     else
505     {
506     // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
507    
508 jfenwick 3117 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0;
509     out->FaceElements->Nodes[INDEX2(1,k,NN)]=n1;
510     out->FaceElements->Nodes[INDEX2(2,k,NN)]=n4;
511 jfenwick 3114
512 jfenwick 3117 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n1;
513     out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n5;
514     out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n4;
515 jfenwick 3114
516     }
517 ksteube 1312 }
518     }
519 jfenwick 3114 faceNECount+=2*local_NE0*local_NE2;
520 ksteube 1312 }
521 jfenwick 3114 totalNECount+=2*NE0*NE2;
522 jfenwick 3221 /* ** elements on boundary 020 (x2=1): - Back*/
523 ksteube 1312 if (local_NE1+e_offset1 == NE1) {
524     #pragma omp parallel for private(i0,i2,k,node0)
525     for (i2=0;i2<local_NE2;i2++) {
526     for (i0=0;i0<local_NE0;i0++) {
527 jfenwick 3114 k=2*(i0+local_NE0*i2)+faceNECount;
528 ksteube 1312 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1)+Nstride2*N_PER_E*(i2+e_offset2);
529 jfenwick 3114 index_t res=2*((i2+e_offset2)+NE2*(i0+e_offset0))+totalNECount;
530     out->FaceElements->Id[k]=res;
531     out->FaceElements->Tag[k]=BACKTAG;
532 ksteube 1312 out->FaceElements->Owner[k]=myRank;
533 jfenwick 3114 out->FaceElements->Id[k+1]=res+1;
534     out->FaceElements->Tag[k+1]=BACKTAG;
535     out->FaceElements->Owner[k+1]=myRank;
536    
537 jfenwick 3117 index_t n2=node0+Nstride1;
538     index_t n6=node0+Nstride1+Nstride2;
539     index_t n7=node0+Nstride0+Nstride1+Nstride2;
540     index_t n3=node0+Nstride0+Nstride1;
541    
542 jfenwick 3114 if ((global_adjustment+i0+local_NE1-1+i2)%2==0)
543     {
544 jfenwick 3117 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n2;
545     out->FaceElements->Nodes[INDEX2(1,k,NN)]=n6;
546     out->FaceElements->Nodes[INDEX2(2,k,NN)]=n3;
547 jfenwick 3114
548 jfenwick 3117 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n6;
549     out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n7;
550     out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n3;
551 jfenwick 3114
552    
553     }
554     else
555     {
556     // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
557 jfenwick 3117 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n2;
558     out->FaceElements->Nodes[INDEX2(1,k,NN)]=n6;
559     out->FaceElements->Nodes[INDEX2(2,k,NN)]=n7;
560 jfenwick 3114
561 jfenwick 3117 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n2;
562     out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n7;
563     out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n3;
564 jfenwick 3114
565     }
566 ksteube 1312 }
567     }
568 jfenwick 3114 faceNECount+=2*local_NE0*local_NE2;
569 ksteube 1312 }
570 jfenwick 3114 totalNECount+=2*NE0*NE2;
571 bcumming 751 }
572 gross 2748 }
573 jfenwick 3086 if (Dudley_noError()) {
574 ksteube 1312 /* add tag names */
575 jfenwick 3114 Dudley_Mesh_addTagMap(out,"top", TOPTAG);
576     Dudley_Mesh_addTagMap(out,"bottom", BOTTOMTAG);
577     Dudley_Mesh_addTagMap(out,"left", LEFTTAG);
578     Dudley_Mesh_addTagMap(out,"right", RIGHTTAG);
579     Dudley_Mesh_addTagMap(out,"front", FRONTTAG);
580     Dudley_Mesh_addTagMap(out,"back", BACKTAG);
581 gross 2748 }
582     /* prepare mesh for further calculatuions:*/
583 jfenwick 3086 if (Dudley_noError()) {
584     Dudley_Mesh_resolveNodeIds(out);
585 gross 2748 }
586 jfenwick 3086 if (Dudley_noError()) {
587     Dudley_Mesh_prepare(out, optimize);
588 bcumming 751 }
589    
590 jfenwick 3086 if (!Dudley_noError()) {
591     Dudley_Mesh_free(out);
592 bcumming 751 }
593 gross 2748 /* free up memory */
594 ksteube 1312 Paso_MPIInfo_free( mpi_info );
595    
596     return out;
597 bcumming 751 }
598 jfenwick 3114
599    

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo

  ViewVC Help
Powered by ViewVC 1.1.26