/[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 2722 - (hide annotations)
Fri Oct 16 06:45:01 2009 UTC (10 years, 5 months ago) by gross
File MIME type: text/plain
File size: 35702 byte(s)
First steps towards macro elements. at the monent a  macro element is identical to its counterpart of 2nd order.




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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26