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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
  ViewVC Help
Powered by ViewVC 1.1.26