/[escript]/trunk/finley/src/Mesh_hex20.c
ViewVC logotype

Annotation of /trunk/finley/src/Mesh_hex20.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1981 - (hide annotations)
Thu Nov 6 05:27:33 2008 UTC (11 years, 1 month ago) by jfenwick
File MIME type: text/plain
File size: 34889 byte(s)
More warning removal.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26