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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26