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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1733 - (hide annotations)
Thu Aug 28 01:46:07 2008 UTC (11 years, 5 months ago) by gross
Original Path: trunk/finley/src/Mesh_hex20.c
File MIME type: text/plain
File size: 34952 byte(s)
bug fix: if the numer of elements to be generated suppose to be zero
still face elements and nodes are generated. this is fixed now.


1 jgs 150
2 ksteube 1312 /* $Id$ */
3    
4     /*******************************************************
5     *
6     * Copyright 2003-2007 by ACceSS MNRF
7     * Copyright 2007 by University of Queensland
8     *
9     * http://esscc.uq.edu.au
10     * Primary Business: Queensland, Australia
11     * Licensed under the Open Software License version 3.0
12     * http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16 jgs 82 /**************************************************************/
17    
18     /* Finley: generates rectangular meshes */
19    
20     /* Generates a numElements[0] x numElements[1] x numElements[2] mesh with second order elements (Hex20) in the brick */
21     /* [0,Length[0]] x [0,Length[1]] x [0,Length[2]]. order is the desired accuracy of the */
22     /* integration scheme. */
23    
24    
25     /**************************************************************/
26    
27     #include "RectangularMesh.h"
28    
29 ksteube 1312 Finley_Mesh* Finley_RectangularMesh_Hex20(dim_t* numElements,
30     double* Length,
31     bool_t* periodic,
32     index_t order,
33     index_t reduced_order,
34     bool_t useElementsOnFace,
35     bool_t useFullElementOrder,
36     bool_t optimize)
37 bcumming 782 {
38 ksteube 1312 #define N_PER_E 2
39     #define DIM 3
40     dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,Nstride0,Nstride1,Nstride2, local_NE0, local_NE1, local_NE2;
41     dim_t totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements, local_N0, local_N1, local_N2, NN;
42     index_t node0, myRank, e_offset0, e_offset1, e_offset2, offset0, offset1, offset2, global_i0, global_i1, global_i2;
43     Finley_Mesh* out;
44     Paso_MPIInfo *mpi_info = NULL;
45     char name[50];
46     double time0=Finley_timer();
47 bcumming 782
48 ksteube 1312 /* get MPI information */
49     mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
50     if (! Finley_noError()) {
51     return NULL;
52 bcumming 782 }
53 ksteube 1312 myRank=mpi_info->rank;
54 bcumming 782
55 ksteube 1312 /* set up the global dimensions of the mesh */
56 bcumming 782
57 jgs 82 NE0=MAX(1,numElements[0]);
58     NE1=MAX(1,numElements[1]);
59     NE2=MAX(1,numElements[2]);
60 ksteube 1312 N0=N_PER_E*NE0+1;
61     N1=N_PER_E*NE1+1;
62     N2=N_PER_E*NE2+1;
63 jgs 82
64 ksteube 1312 /* allocate mesh: */
65     sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
66     out=Finley_Mesh_alloc(name,DIM,order, reduced_order, mpi_info);
67     if (! Finley_noError()) {
68     Paso_MPIInfo_free( mpi_info );
69     return NULL;
70     }
71    
72     if (useFullElementOrder) {
73 gross 1414 /* Finley_setError(SYSTEM_ERROR,"full element order for Hex elements is not supported yet."); */
74 ksteube 1312 Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Hex27,
75     out->order,
76     out->reduced_order,
77 gross 1414 mpi_info));
78 ksteube 1312 if (useElementsOnFace) {
79     Finley_setError(SYSTEM_ERROR,"rich elements for Hex27 elements is not supported yet.");
80 jgs 153 } else {
81 ksteube 1312 Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Rec9,
82     out->order,
83     out->reduced_order,
84     mpi_info));
85     Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Rec9_Contact,
86     out->order,
87     out->reduced_order,
88     mpi_info));
89 jgs 153 }
90 ksteube 1312
91     } else {
92     Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Hex20,out->order,out->reduced_order,mpi_info));
93     if (useElementsOnFace) {
94     Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Hex20Face,
95     out->order,
96     out->reduced_order,
97     mpi_info));
98     Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Hex20Face_Contact,
99     out->order,
100     out->reduced_order,
101     mpi_info));
102 jgs 153 } else {
103 ksteube 1312 Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Rec8,
104     out->order,
105     out->reduced_order,
106     mpi_info));
107     Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Rec8_Contact,
108     out->order,
109     out->reduced_order,
110     mpi_info));
111 jgs 153 }
112     }
113 ksteube 1312 Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(Point1,
114     out->order,
115     out->reduced_order,
116     mpi_info));
117 bcumming 782 if (! Finley_noError()) {
118 ksteube 1312 Paso_MPIInfo_free( mpi_info );
119     Finley_Mesh_free(out);
120 bcumming 782 return NULL;
121     }
122    
123 ksteube 1312 /* work out the largest dimension */
124     if (N2==MAX3(N0,N1,N2)) {
125     Nstride0=1;
126     Nstride1=N0;
127     Nstride2=N0*N1;
128     local_NE0=NE0;
129     e_offset0=0;
130     local_NE1=NE1;
131     e_offset1=0;
132     Paso_MPIInfo_Split(mpi_info,NE2,&local_NE2,&e_offset2);
133     } else if (N1==MAX3(N0,N1,N2)) {
134     Nstride0=N2;
135     Nstride1=N0*N2;
136     Nstride2=1;
137     local_NE0=NE0;
138     e_offset0=0;
139     Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
140     local_NE2=NE2;
141     e_offset2=0;
142 jgs 82 } else {
143 ksteube 1312 Nstride0=N1*N2;
144     Nstride1=1;
145     Nstride2=N1;
146     Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
147     local_NE1=NE1;
148     e_offset1=0;
149     local_NE2=NE2;
150     e_offset2=0;
151 jgs 82 }
152 ksteube 1312 offset0=e_offset0*N_PER_E;
153     offset1=e_offset1*N_PER_E;
154     offset2=e_offset2*N_PER_E;
155 gross 1733 local_N0=local_NE0>0 ? local_NE0*N_PER_E+1 : 0;
156     local_N1=local_NE1>0 ? local_NE1*N_PER_E+1 : 0;
157     local_N2=local_NE0>0 ? local_NE2*N_PER_E+1 : 0;
158 jgs 82
159 ksteube 1312 /* get the number of surface elements */
160 jgs 153
161 ksteube 1312 NFaceElements=0;
162 gross 1733 if (!periodic[2] && (local_NE2>0) ) {
163 ksteube 1312 NDOF2=N2;
164     if (offset2==0) NFaceElements+=local_NE1*local_NE0;
165     if (local_NE2+e_offset2 == NE2) NFaceElements+=local_NE1*local_NE0;
166     } else {
167     NDOF2=N2-1;
168 jgs 82 }
169 ksteube 1312
170 gross 1733 if (!periodic[0] && (local_NE1>0) ) {
171 ksteube 1312 NDOF0=N0;
172     if (e_offset0 == 0) NFaceElements+=local_NE1*local_NE2;
173     if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1*local_NE2;
174     } else {
175     NDOF0=N0-1;
176 jgs 82 }
177 gross 1733 if (!periodic[1] && (local_NE1>0) ) {
178 ksteube 1312 NDOF1=N1;
179     if (e_offset1 == 0) NFaceElements+=local_NE0*local_NE2;
180     if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0*local_NE2;
181 jgs 82 } else {
182 ksteube 1312 NDOF1=N1-1;
183 jgs 82 }
184 gross 1733 printf("HEX20 %d %d %d\n",local_NE0,local_NE1,local_NE2);
185 jgs 82
186 ksteube 1312 /* allocate tables: */
187 jgs 82
188 ksteube 1312 Finley_NodeFile_allocTable(out->Nodes,local_N0*local_N1*local_N2);
189     Finley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1*local_NE2);
190     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
191 jgs 82
192 ksteube 1312 if (Finley_noError()) {
193     /* create nodes */
194 jgs 82
195 ksteube 1312 #pragma omp parallel for private(i0,i1,i2,k,global_i0,global_i1,global_i2)
196     for (i2=0;i2<local_N2;i2++) {
197     for (i1=0;i1<local_N1;i1++) {
198     for (i0=0;i0<local_N0;i0++) {
199     k=i0+local_N0*i1+local_N0*local_N1*i2;
200     global_i0=i0+offset0;
201     global_i1=i1+offset1;
202     global_i2=i2+offset2;
203     out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
204     out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
205     out->Nodes->Coordinates[INDEX2(2,k,DIM)]=DBLE(global_i2)/DBLE(N2-1)*Length[2];
206     out->Nodes->Id[k]=Nstride0*global_i0+Nstride1*global_i1+Nstride2*global_i2;
207     out->Nodes->Tag[k]=0;
208     out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
209     +Nstride1*(global_i1%NDOF1)
210     +Nstride2*(global_i2%NDOF2);
211 jgs 82 }
212     }
213     }
214 ksteube 1312 /* set the elements: */
215     NN=out->Elements->numNodes;
216     #pragma omp parallel for private(i0,i1,i2,k,node0)
217     for (i2=0;i2<local_NE2;i2++) {
218     for (i1=0;i1<local_NE1;i1++) {
219     for (i0=0;i0<local_NE0;i0++) {
220    
221     k=i0+local_NE0*i1+local_NE0*local_NE1*i2;
222     node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
223 jgs 82
224 ksteube 1312 out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+NE0*NE1*(i2+e_offset2);
225     out->Elements->Tag[k]=0;
226     out->Elements->Owner[k]=myRank;
227 jgs 82
228 ksteube 1312 out->Elements->Nodes[INDEX2(0,k,NN)] =node0 ;
229     out->Elements->Nodes[INDEX2(1,k,NN)] =node0+ 2*Nstride0;
230     out->Elements->Nodes[INDEX2(2,k,NN)] =node0+ 2*Nstride1+2*Nstride0;
231     out->Elements->Nodes[INDEX2(3,k,NN)] =node0+ 2*Nstride1;
232     out->Elements->Nodes[INDEX2(4,k,NN)] =node0+2*Nstride2 ;
233     out->Elements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride2 +2*Nstride0;
234     out->Elements->Nodes[INDEX2(6,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
235     out->Elements->Nodes[INDEX2(7,k,NN)] =node0+2*Nstride2+2*Nstride1 ;
236     out->Elements->Nodes[INDEX2(8,k,NN)] =node0+ 1*Nstride0;
237     out->Elements->Nodes[INDEX2(9,k,NN)] =node0+ 1*Nstride1+2*Nstride0;
238     out->Elements->Nodes[INDEX2(10,k,NN)]=node0+ 2*Nstride1+1*Nstride0;
239     out->Elements->Nodes[INDEX2(11,k,NN)]=node0+ 1*Nstride1 ;
240     out->Elements->Nodes[INDEX2(12,k,NN)]=node0+1*Nstride2 ;
241     out->Elements->Nodes[INDEX2(13,k,NN)]=node0+1*Nstride2 +2*Nstride0;
242     out->Elements->Nodes[INDEX2(14,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
243     out->Elements->Nodes[INDEX2(15,k,NN)]=node0+1*Nstride2+2*Nstride1 ;
244     out->Elements->Nodes[INDEX2(16,k,NN)]=node0+2*Nstride2 +1*Nstride0;
245     out->Elements->Nodes[INDEX2(17,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
246     out->Elements->Nodes[INDEX2(18,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
247     out->Elements->Nodes[INDEX2(19,k,NN)]=node0+2*Nstride2+1*Nstride1 ;
248     if (useFullElementOrder) {
249     out->Elements->Nodes[INDEX2(20,k,NN)]=node0+ 1*Nstride1+1*Nstride0;
250     out->Elements->Nodes[INDEX2(21,k,NN)]=node0+1*Nstride2 +1*Nstride0;
251     out->Elements->Nodes[INDEX2(22,k,NN)]=node0+1*Nstride2+1*Nstride1+2*Nstride0;
252     out->Elements->Nodes[INDEX2(23,k,NN)]=node0+1*Nstride2+2*Nstride1+1*Nstride0;
253     out->Elements->Nodes[INDEX2(24,k,NN)]=node0+1*Nstride2+1*Nstride1 ;
254     out->Elements->Nodes[INDEX2(25,k,NN)]=node0+2*Nstride2+1*Nstride1+1*Nstride0;
255     out->Elements->Nodes[INDEX2(26,k,NN)]=node0+1*Nstride2+1*Nstride1+1*Nstride0;
256     }
257 jgs 82 }
258     }
259     }
260 ksteube 1312 /* face elements */
261     NN=out->FaceElements->numNodes;
262     totalNECount=NE0*NE1*NE2;
263     faceNECount=0;
264     /* these are the quadrilateral elements on boundary 1 (x3=0): */
265 gross 1733 if (!periodic[2] && (local_NE2>0)) {
266 ksteube 1312 /* ** elements on boundary 100 (x3=0): */
267     if (offset2==0) {
268     #pragma omp parallel for private(i0,i1,k,node0)
269     for (i1=0;i1<local_NE1;i1++) {
270     for (i0=0;i0<local_NE0;i0++) {
271    
272     k=i0+local_NE0*i1+faceNECount;
273     node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
274    
275     out->FaceElements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+totalNECount;
276     out->FaceElements->Tag[k]=100;
277     out->FaceElements->Owner[k]=myRank;
278    
279     if (useElementsOnFace) {
280     out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0 ;
281     out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0 +2*Nstride1 ;
282     out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0 +2*Nstride1+2*Nstride0;
283     out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+ 2*Nstride0 ;
284     out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+2*Nstride2 ;
285     out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride2+2*Nstride1 ;
286     out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
287     out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+2*Nstride2 +2*Nstride0;
288     out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+ 1*Nstride1 ;
289     out->FaceElements->Nodes[INDEX2(9,k,NN)] =node0+ 2*Nstride1+1*Nstride0;
290     out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+ 1*Nstride1+2*Nstride0;
291     out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+ 1*Nstride0;
292     out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+1*Nstride2 ;
293     out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+1*Nstride2+2*Nstride1 ;
294     out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
295     out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+1*Nstride2 +2*Nstride0;
296     out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+2*Nstride2+1*Nstride1;
297     out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
298     out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
299     out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+2*Nstride2 +1*Nstride0;
300     } else {
301     out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0 ;
302     out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+ 2*Nstride1 ;
303     out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+ 2*Nstride1+2*Nstride0;
304     out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+ 2*Nstride0;
305     out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+ 1*Nstride1 ;
306     out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+ 2*Nstride1+1*Nstride0;
307     out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+ 1*Nstride1+2*Nstride0;
308     out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+ 1*Nstride0;
309     if (useFullElementOrder){
310     out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+ 1*Nstride1+1*Nstride0;
311     }
312     }
313     }
314     }
315     faceNECount+=local_NE1*local_NE0;
316 jgs 82 }
317 ksteube 1312 totalNECount+=NE1*NE0;
318     /* ** elements on boundary 200 (x3=1) */
319     if (local_NE2+e_offset2 == NE2) {
320     #pragma omp parallel for private(i0,i1,k,node0)
321     for (i1=0;i1<local_NE1;i1++) {
322     for (i0=0;i0<local_NE0;i0++) {
323    
324     k=i0+local_NE0*i1+faceNECount;
325     node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(NE2-1);
326    
327     out->FaceElements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+totalNECount;
328     out->FaceElements->Tag[k]=200;
329     out->FaceElements->Owner[k]=myRank;
330     if (useElementsOnFace) {
331     out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0+2*Nstride2 ;
332     out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+2*Nstride2+ 2*Nstride0;
333     out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
334     out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+2*Nstride2+2*Nstride1 ;
335    
336     out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0 ;
337     out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride0 ;
338     out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+ 2*Nstride1+2*Nstride0;
339     out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+ 2*Nstride1;
340    
341     out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+2*Nstride2+ 1*Nstride0;
342     out->FaceElements->Nodes[INDEX2(9,k,NN)] =node0+2*Nstride2+1*Nstride1+2*Nstride0;
343     out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
344     out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+2*Nstride2+1*Nstride1 ;
345    
346     out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+1*Nstride2;
347     out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+1*Nstride2 +2*Nstride0;
348     out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
349     out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+1*Nstride2+2*Nstride1 ;
350    
351     out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+ 1*Nstride0;
352     out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+ 1*Nstride1+2*Nstride0;
353     out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+ 2*Nstride1+1*Nstride0;
354     out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+ 1*Nstride1 ;
355    
356     } else {
357     out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0+2*Nstride2 ;
358     out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+2*Nstride2 +2*Nstride0;
359     out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
360     out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+2*Nstride2+2*Nstride1 ;
361     out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+2*Nstride2 +1*Nstride0;
362     out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride2+1*Nstride1+2*Nstride0;
363     out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+2*Nstride2+2*Nstride1+1*Nstride0;
364     out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+2*Nstride2+1*Nstride1 ;
365     if (useFullElementOrder){
366     out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+2*Nstride2+1*Nstride1+1*Nstride0;
367     }
368     }
369     }
370     }
371     faceNECount+=local_NE1*local_NE0;
372     }
373     totalNECount+=NE1*NE0;
374 jgs 82 }
375 gross 1733 if (!periodic[0] && (local_NE0>0)) {
376 ksteube 1312 /* ** elements on boundary 001 (x1=0): */
377 jgs 82
378 ksteube 1312 if (e_offset0 == 0) {
379     #pragma omp parallel for private(i1,i2,k,node0)
380     for (i2=0;i2<local_NE2;i2++) {
381     for (i1=0;i1<local_NE1;i1++) {
382    
383     k=i1+local_NE1*i2+faceNECount;
384     node0=Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
385     out->FaceElements->Id[k]=(i1+e_offset1)+NE1*(i2+e_offset2)+totalNECount;
386     out->FaceElements->Tag[k]=1;
387     out->FaceElements->Owner[k]=myRank;
388    
389     if (useElementsOnFace) {
390     out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0 ;
391     out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+2*Nstride2 ;
392     out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+2*Nstride2+2*Nstride1 ;
393     out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+2*Nstride1 ;
394    
395     out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+2*Nstride0 ;
396     out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride2+2*Nstride0 ;
397     out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
398     out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+2*Nstride1+2*Nstride0 ;
399    
400     out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+1*Nstride2 ;
401     out->FaceElements->Nodes[INDEX2(9,k,NN)] =node0+2*Nstride2+1*Nstride1 ;
402     out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+1*Nstride2+2*Nstride1 ;
403     out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+ 1*Nstride1 ;
404    
405     out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+ 1*Nstride0;
406     out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+2*Nstride2 +1*Nstride0;
407     out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
408     out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+2*Nstride1+ 1*Nstride0;
409    
410     out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+1*Nstride2+ 2*Nstride0;
411     out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
412     out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
413     out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+1*Nstride1+ 2*Nstride0;
414    
415     } else {
416     out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0 ;
417     out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+2*Nstride2 ;
418     out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+2*Nstride2+2*Nstride1 ;
419     out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+ 2*Nstride1 ;
420     out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+1*Nstride2 ;
421     out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride2+1*Nstride1 ;
422     out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+1*Nstride2+2*Nstride1 ;
423     out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+ 1*Nstride1 ;
424     if (useFullElementOrder){
425     out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+1*Nstride2+1*Nstride1 ;
426     }
427     }
428     }
429     }
430     faceNECount+=local_NE1*local_NE2;
431     }
432     totalNECount+=NE1*NE2;
433    
434     /* ** elements on boundary 002 (x1=1): */
435     if (local_NE0+e_offset0 == NE0) {
436     #pragma omp parallel for private(i1,i2,k,node0)
437     for (i2=0;i2<local_NE2;i2++) {
438     for (i1=0;i1<local_NE1;i1++) {
439     k=i1+local_NE1*i2+faceNECount;
440     node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
441     out->FaceElements->Id[k]=(i1+e_offset1)+NE1*(i2+e_offset2)+totalNECount;
442     out->FaceElements->Tag[k]=2;
443     out->FaceElements->Owner[k]=myRank;
444 jgs 82
445 ksteube 1312 if (useElementsOnFace) {
446     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+ 2*Nstride0;
447     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+ 2*Nstride1+2*Nstride0;
448     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
449     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride2+ 2*Nstride0;
450    
451     out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0 ;
452     out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+ 2*Nstride1 ;
453     out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride2+2*Nstride1 ;
454     out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+2*Nstride2 ;
455    
456     out->FaceElements->Nodes[INDEX2(8,k,NN)]=node0+ 1*Nstride1+2*Nstride0;
457     out->FaceElements->Nodes[INDEX2(9,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
458     out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
459     out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+1*Nstride2+ 2*Nstride0;
460    
461     out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+ 1*Nstride0;
462     out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+ 2*Nstride1+1*Nstride0;
463     out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
464     out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+2*Nstride2+ 1*Nstride0;
465    
466     out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+ 1*Nstride1 ;
467     out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+1*Nstride2+2*Nstride1 ;
468     out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+2*Nstride2+1*Nstride1 ;
469     out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+1*Nstride2 ;
470    
471     } else {
472     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0 +2*Nstride0;
473     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+ 2*Nstride1+2*Nstride0;
474     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
475     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride2+ 2*Nstride0;
476     out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+ 1*Nstride1+2*Nstride0;
477     out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
478     out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
479     out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+1*Nstride2 +2*Nstride0;
480     if (useFullElementOrder){
481     out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+1*Nstride2+1*Nstride1+2*Nstride0;
482     }
483     }
484    
485     }
486     }
487     faceNECount+=local_NE1*local_NE2;
488 jgs 82 }
489 ksteube 1312 totalNECount+=NE1*NE2;
490 jgs 82 }
491 gross 1733 if (!periodic[1] && (local_NE1>0)) {
492 ksteube 1312 /* ** elements on boundary 010 (x2=0): */
493     if (e_offset1 == 0) {
494     #pragma omp parallel for private(i0,i2,k,node0)
495     for (i2=0;i2<local_NE2;i2++) {
496     for (i0=0;i0<local_NE0;i0++) {
497     k=i0+local_NE0*i2+faceNECount;
498     node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride2*N_PER_E*(i2+e_offset2);
499    
500     out->FaceElements->Id[k]=(i2+e_offset2)+NE2*(e_offset0+i0)+totalNECount;
501     out->FaceElements->Tag[k]=10;
502     out->FaceElements->Owner[k]=myRank;
503     if (useElementsOnFace) {
504     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0 ;
505     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+ 2*Nstride0;
506     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2 +2*Nstride0;
507     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride2 ;
508    
509     out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+ 2*Nstride1 ;
510     out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+2*Nstride1+ 2*Nstride0;
511     out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
512     out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+2*Nstride2+2*Nstride1 ;
513    
514     out->FaceElements->Nodes[INDEX2(8,k,NN)]=node0+ 1*Nstride0;
515     out->FaceElements->Nodes[INDEX2(9,k,NN)]=node0+1*Nstride2+ 2*Nstride0;
516     out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+2*Nstride2+ 1*Nstride0;
517     out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+1*Nstride2 ;
518    
519     out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+ 1*Nstride1 ;
520     out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+ 1*Nstride1+2*Nstride0;
521     out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
522     out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+2*Nstride2+1*Nstride1 ;
523 bcumming 782
524 ksteube 1312 out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+ 2*Nstride1+1*Nstride0;
525     out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
526     out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
527     out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+1*Nstride2+2*Nstride1 ;
528    
529     } else {
530     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0 ;
531     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+ 2*Nstride0;
532     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2+ 2*Nstride0;
533     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride2 ;
534     out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+ 1*Nstride0;
535     out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+1*Nstride2+ 2*Nstride0;
536     out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride2+ 1*Nstride0;
537     out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+1*Nstride2 ;
538     if (useFullElementOrder){
539     out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+1*Nstride2+ 1*Nstride0;
540     }
541     }
542     }
543     }
544     faceNECount+=local_NE0*local_NE2;
545     }
546     totalNECount+=NE0*NE2;
547     /* ** elements on boundary 020 (x2=1): */
548     if (local_NE1+e_offset1 == NE1) {
549     #pragma omp parallel for private(i0,i2,k,node0)
550     for (i2=0;i2<local_NE2;i2++) {
551     for (i0=0;i0<local_NE0;i0++) {
552     k=i0+local_NE0*i2+faceNECount;
553     node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1)+Nstride2*N_PER_E*(i2+e_offset2);
554 bcumming 782
555 ksteube 1312 out->FaceElements->Id[k]=(i2+e_offset2)+NE2*(i0+e_offset0)+totalNECount;
556     out->FaceElements->Tag[k]=20;
557     out->FaceElements->Owner[k]=myRank;
558    
559     if (useElementsOnFace) {
560     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+ 2*Nstride1 ;
561     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride2+2*Nstride1 ;
562     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
563     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1+2*Nstride0 ;
564    
565     out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0 ;
566     out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+2*Nstride2 ;
567     out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride2+ 2*Nstride0;
568     out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+ 2*Nstride0;
569    
570     out->FaceElements->Nodes[INDEX2(8,k,NN)]=node0+1*Nstride2+2*Nstride1 ;
571     out->FaceElements->Nodes[INDEX2(9,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
572     out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
573     out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+ 2*Nstride1+1*Nstride0;
574    
575     out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+ 1*Nstride1 ;
576     out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+2*Nstride2+1*Nstride1 ;
577     out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
578     out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+ 1*Nstride1+2*Nstride0;
579    
580     out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+1*Nstride2 ;
581     out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+2*Nstride2 +1*Nstride0;
582     out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+1*Nstride2 +2*Nstride0;
583     out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+ 1*Nstride0;
584     } else {
585     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+ 2*Nstride1 ;
586     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride2+2*Nstride1 ;
587     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
588     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+ 2*Nstride1+2*Nstride0;
589     out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+1*Nstride2+2*Nstride1 ;
590     out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
591     out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
592     out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+ 2*Nstride1+1*Nstride0;
593     if (useFullElementOrder){
594     out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+1*Nstride2+2*Nstride1+1*Nstride0;
595     }
596     }
597     }
598     }
599     faceNECount+=local_NE0*local_NE2;
600     }
601     totalNECount+=NE0*NE2;
602 bcumming 782 }
603 ksteube 1312 /* add tag names */
604     Finley_Mesh_addTagMap(out,"top", 200);
605     Finley_Mesh_addTagMap(out,"bottom", 100);
606     Finley_Mesh_addTagMap(out,"left", 1);
607     Finley_Mesh_addTagMap(out,"right", 2);
608     Finley_Mesh_addTagMap(out,"front", 10);
609     Finley_Mesh_addTagMap(out,"back", 20);
610 bcumming 782
611 ksteube 1312 /* prepare mesh for further calculatuions:*/
612     if (Finley_noError()) {
613     Finley_Mesh_resolveNodeIds(out);
614 bcumming 782 }
615 ksteube 1312 if (Finley_noError()) {
616     Finley_Mesh_prepare(out, optimize);
617 bcumming 782 }
618     }
619    
620 ksteube 1312 if (!Finley_noError()) {
621     Finley_Mesh_free(out);
622     }
623 bcumming 782 /* free up memory */
624 ksteube 1312 Paso_MPIInfo_free( mpi_info );
625 bcumming 782 #ifdef Finley_TRACE
626     printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
627     #endif
628 ksteube 1312
629     return out;
630 bcumming 782 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26