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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1062 - (hide annotations)
Mon Mar 26 06:17:53 2007 UTC (12 years, 6 months ago) by gross
File MIME type: text/plain
File size: 38037 byte(s)
reduced integration schemes are implemented now for grad, integrate, etc. Tests still to be added.
1 jgs 150 /*
2 elspeth 616 ************************************************************
3     * Copyright 2006 by ACcESS MNRF *
4     * *
5     * http://www.access.edu.au *
6     * Primary Business: Queensland, Australia *
7     * Licensed under the Open Software License version 3.0 *
8     * http://www.opensource.org/licenses/osl-3.0.php *
9     * *
10     ************************************************************
11 jgs 150 */
12    
13 jgs 82 /**************************************************************/
14    
15     /* Finley: generates rectangular meshes */
16    
17     /* Generates a numElements[0] x numElements[1] x numElements[2] mesh with first order elements (Hex8) in the brick */
18     /* [0,Length[0]] x [0,Length[1]] x [0,Length[2]]. order is the desired accuracy of the */
19     /* integration scheme. */
20    
21     /**************************************************************/
22    
23 jgs 150 /* Author: gross@access.edu.au */
24     /* Version: $Id$ */
25 jgs 82
26     /**************************************************************/
27    
28     #include "RectangularMesh.h"
29    
30     /**************************************************************/
31    
32 bcumming 751 #ifdef PASO_MPI
33     /* get the number of nodes/elements for domain with rank=rank, of size processors
34     where n is the total number of nodes/elements in the global domain */
35     static index_t domain_MODdim( index_t rank, index_t size, index_t n )
36     {
37     rank = size-rank-1;
38    
39     if( rank < n%size )
40     return (index_t)floor(n/size)+1;
41     return (index_t)floor(n/size);
42     }
43    
44    
45     /* Determines the number of nodes/elements etc along an axis which is numElementsGlobal long for domain rank */
46     /* A bit messy, but it only has to be done once... */
47     static void domain_calculateDimension( index_t rank, dim_t size, dim_t numElementsGlobal, bool_t periodic, dim_t *numNodesLocal, dim_t *numDOFLocal, dim_t *numElementsLocal, dim_t *numElementsInternal, dim_t *firstNode, dim_t *nodesExternal, dim_t *DOFExternal, dim_t *numNodesExternal, bool_t *periodicLocal )
48     {
49     index_t i0;
50     dim_t numNodesGlobal = numElementsGlobal+1;
51    
52     (*numNodesLocal) = domain_MODdim( rank, size, numNodesGlobal );
53    
54     numElementsLocal[0] = numNodesLocal[0]+1;
55     periodicLocal[0] = periodicLocal[1] = FALSE;
56     nodesExternal[0] = nodesExternal[1] = 1;
57     if( periodic )
58     {
59     if( size==1 )
60     {
61     numElementsLocal[0] = numElementsGlobal;
62     nodesExternal[0] = nodesExternal[1] = 0;
63     periodicLocal[0] = periodicLocal[1] = TRUE;
64     }
65     else
66     {
67     if( rank==0 )
68     {
69     periodicLocal[0] = TRUE;
70     numNodesLocal[0]++;
71     }
72     if( rank==(size-1) )
73     {
74     periodicLocal[1] = TRUE;
75     numNodesLocal[0]--;
76     numElementsLocal[0]--;
77     }
78     }
79     }
80     else if( !periodic )
81     {
82     if( rank==0 ){
83     nodesExternal[0]--;
84     numElementsLocal[0]--;
85     }
86     if( rank==(size-1) )
87     {
88     nodesExternal[1]--;
89     numElementsLocal[0]--;
90     }
91     }
92     numNodesExternal[0] = nodesExternal[0]+nodesExternal[1];
93    
94     numElementsInternal[0] = numElementsLocal[0];
95     if( (rank==0) && (rank==size-1) );
96    
97     else if( !periodic && ( (rank==0) ^ (rank==size-1) ) )
98     numElementsInternal[0] -= 1;
99     else
100     numElementsInternal[0] -= 2;
101    
102     firstNode[0] = 0;
103     for( i0=0; i0<rank; i0++ )
104     firstNode[0] += domain_MODdim( i0, size, numNodesGlobal );
105    
106     numDOFLocal[0] = numNodesLocal[0];
107     if( periodicLocal[0] )
108     {
109     numDOFLocal[0]--;
110     }
111     DOFExternal[0] = nodesExternal[0];
112     DOFExternal[1] = nodesExternal[1];
113     }
114    
115     #endif
116    
117 bcumming 782 #ifdef PASO_MPI
118 gross 1062 Finley_Mesh* Finley_RectangularMesh_Hex8_singleCPU(dim_t* numElements,double* Length,bool_t* periodic, index_t order, index_t reduced_order, bool_t useElementsOnFace,Paso_MPIInfo *mpi_info)
119 bcumming 782 #else
120 gross 1062 Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order, index_t reduced_order, bool_t useElementsOnFace)
121 bcumming 751 #endif
122     {
123 jgs 153 dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES,M0,M1,M2;
124 jgs 123 index_t node0;
125 jgs 82 Finley_Mesh* out;
126     char name[50];
127     double time0=Finley_timer();
128     NE0=MAX(1,numElements[0]);
129     NE1=MAX(1,numElements[1]);
130     NE2=MAX(1,numElements[2]);
131     N0=NE0+1;
132     N1=NE1+1;
133     N2=NE2+1;
134    
135 jgs 153 if (N0<=MIN(N1,N2)) {
136     if (N1 <= N2) {
137     M0=1;
138     M1=N0;
139     M2=N0*N1;
140     } else {
141     M0=1;
142     M2=N0;
143     M1=N0*N2;
144     }
145     } else if (N1<=MIN(N2,N0)) {
146     if (N2 <= N0) {
147     M1=1;
148     M2=N1;
149     M0=N2*N1;
150     } else {
151     M1=1;
152     M0=N1;
153     M2=N1*N0;
154     }
155     } else {
156     if (N0 <= N1) {
157     M2=1;
158     M0=N2;
159     M1=N2*N0;
160     } else {
161     M2=1;
162     M1=N2;
163     M0=N1*N2;
164     }
165     }
166    
167    
168 jgs 82 NFaceElements=0;
169     if (!periodic[0]) {
170     NDOF0=N0;
171     NFaceElements+=2*NE1*NE2;
172     } else {
173     NDOF0=N0-1;
174     }
175     if (!periodic[1]) {
176     NDOF1=N1;
177     NFaceElements+=2*NE0*NE2;
178     } else {
179     NDOF1=N1-1;
180     }
181     if (!periodic[2]) {
182     NDOF2=N2;
183     NFaceElements+=2*NE0*NE1;
184     } else {
185     NDOF2=N2-1;
186     }
187    
188     /* allocate mesh: */
189    
190     sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
191 bcumming 751
192 bcumming 730 #ifndef PASO_MPI
193 gross 1062 out=Finley_Mesh_alloc(name,3,order,reduced_order);
194 bcumming 730 #else
195 gross 1062 out=Finley_Mesh_alloc(name,3,order,reduced_order,mpi_info);
196 bcumming 730 #endif
197 jgs 150 if (! Finley_noError()) return NULL;
198 jgs 82
199 bcumming 751 #ifdef PASO_MPI
200 gross 1062 out->Elements=Finley_ElementFile_alloc(Hex8,out->order,out->reduced_order,mpi_info);
201 bcumming 751 if (useElementsOnFace) {
202 gross 1062 out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,out->reduced_order,mpi_info);
203     out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,out->reduced_order,mpi_info);
204 bcumming 751 } else {
205 gross 1062 out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,out->reduced_order,mpi_info);
206     out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,out->reduced_order,mpi_info);
207 bcumming 751 }
208 gross 1062 out->Points=Finley_ElementFile_alloc(Point1,out->order,out->reduced_order,mpi_info);
209 bcumming 751 #else
210 gross 1062 out->Elements=Finley_ElementFile_alloc(Hex8,out->order, out->reduced_order);
211 jgs 82 if (useElementsOnFace) {
212 gross 1062 out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,out->reduced_order);
213     out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order, out->reduced_order);
214 jgs 82 } else {
215 gross 1062 out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order, out->reduced_order);
216     out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order, out->reduced_order);
217 jgs 82 }
218 gross 1062 out->Points=Finley_ElementFile_alloc(Point1,out->order, out->reduced_order);
219 bcumming 730 #endif
220 jgs 150 if (! Finley_noError()) {
221 jgs 82 Finley_Mesh_dealloc(out);
222     return NULL;
223     }
224    
225    
226 bcumming 751 /* allocate tables: */
227 jgs 82 Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
228 bcumming 788 Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
229     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
230 bcumming 751 #ifdef PASO_MPI
231     Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0*NDOF1*NDOF2, 0, 0 );
232 bcumming 788 Finley_ElementDistribution_allocTable( out->Elements->elementDistribution, NE0*NE1*NE2, NE0*NE1*NE2);
233     Finley_ElementDistribution_allocTable( out->FaceElements->elementDistribution, NFaceElements, NFaceElements );
234 bcumming 730 #endif
235 jgs 150 if (! Finley_noError()) {
236 jgs 82 Finley_Mesh_dealloc(out);
237     return NULL;
238     }
239    
240     /* set nodes: */
241 jgs 153
242     #pragma omp parallel for private(i0,i1,i2,k)
243 jgs 82 for (i2=0;i2<N2;i2++) {
244     for (i1=0;i1<N1;i1++) {
245     for (i0=0;i0<N0;i0++) {
246 jgs 153 k=M0*i0+M1*i1+M2*i2;
247 jgs 82 out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0)/DBLE(N0-1)*Length[0];
248     out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
249     out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
250 jgs 153 out->Nodes->Id[k]=i0+N0*i1+N0*N1*i2;
251 jgs 82 out->Nodes->Tag[k]=0;
252 jgs 153 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1) +M2*(i2%NDOF2);
253 bcumming 782 #ifdef PASO_MPI
254     out->Nodes->Dom[k]=NODE_INTERNAL;
255     #endif
256 jgs 82 }
257     }
258     }
259     /* tags for the faces: */
260 jgs 153 /* tags for the faces: */
261 jgs 82 if (!periodic[2]) {
262 jgs 153 for (i1=0;i1<N1;i1++) {
263     for (i0=0;i0<N0;i0++) {
264     out->Nodes->Tag[M0*i0+M1*i1+M2*0]+=100;
265     out->Nodes->Tag[M0*i0+M1*i1+M2*(N2-1)]+=200;
266     }
267     }
268 jgs 82 }
269     if (!periodic[1]) {
270     for (i2=0;i2<N2;i2++) {
271     for (i0=0;i0<N0;i0++) {
272 jgs 153 out->Nodes->Tag[M0*i0+M1*0+M2*i2]+=10;
273     out->Nodes->Tag[M0*i0+M1*(N1-1)+M2*i2]+=20;
274 jgs 82 }
275     }
276     }
277     if (!periodic[0]) {
278     for (i2=0;i2<N2;i2++) {
279     for (i1=0;i1<N1;i1++) {
280 jgs 153 out->Nodes->Tag[M0*0+M1*i1+M2*i2]+=1;
281     out->Nodes->Tag[M0*(N0-1)+M1*i1+M2*i2]+=2;
282 jgs 82 }
283     }
284     }
285 jgs 153
286 jgs 82 /* set the elements: */
287    
288     #pragma omp parallel for private(i0,i1,i2,k,node0)
289     for (i2=0;i2<NE2;i2++) {
290     for (i1=0;i1<NE1;i1++) {
291     for (i0=0;i0<NE0;i0++) {
292     k=i0+NE0*i1+NE0*NE1*i2;
293     node0=i0+i1*N0+N0*N1*i2;
294    
295     out->Elements->Id[k]=k;
296     out->Elements->Tag[k]=0;
297     out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
298 bcumming 782 #ifdef PASO_MPI
299     out->Elements->Dom[k]=ELEMENT_INTERNAL;
300     #endif
301 jgs 82
302     out->Elements->Nodes[INDEX2(0,k,8)]=node0;
303     out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
304     out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0+1;
305     out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0;
306     out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0*N1;
307     out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0*N1+1;
308     out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0*N1+N0+1;
309     out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0*N1+N0;
310    
311     }
312     }
313     }
314 jgs 123 out->Elements->minColor=0;
315     out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
316 jgs 82
317     /* face elements: */
318    
319     if (useElementsOnFace) {
320     NUMNODES=8;
321     } else {
322     NUMNODES=4;
323     }
324     totalNECount=NE0*NE1*NE2;
325     faceNECount=0;
326    
327     /* these are the quadrilateral elements on boundary 1 (x3=0): */
328     if (!periodic[2]) {
329     /* ** elements on boundary 100 (x3=0): */
330    
331     #pragma omp parallel for private(i0,i1,k,node0)
332     for (i1=0;i1<NE1;i1++) {
333     for (i0=0;i0<NE0;i0++) {
334     k=i0+NE0*i1+faceNECount;
335     node0=i0+i1*N0;
336    
337     out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
338     out->FaceElements->Tag[k]=100;
339     out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
340 bcumming 782 #ifdef PASO_MPI
341     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
342     #endif
343 jgs 82
344     if (useElementsOnFace) {
345     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
346     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
347     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
348     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
349     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;
350     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+N0;
351     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
352     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+1;
353     } else {
354     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
355     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
356     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
357     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
358     }
359     }
360     }
361     totalNECount+=NE1*NE0;
362     faceNECount+=NE1*NE0;
363    
364     /* ** elements on boundary 200 (x3=1) */
365    
366     #pragma omp parallel for private(i0,i1,k,node0)
367     for (i1=0;i1<NE1;i1++) {
368     for (i0=0;i0<NE0;i0++) {
369     k=i0+NE0*i1+faceNECount;
370     node0=i0+i1*N0+N0*N1*(NE2-1);
371    
372     out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
373     out->FaceElements->Tag[k]=200;
374     out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
375 bcumming 782 #ifdef PASO_MPI
376     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
377     #endif
378 jgs 82
379     if (useElementsOnFace) {
380     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
381     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
382     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
383     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
384     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
385     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
386     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+1;
387     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;
388     } else {
389     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
390     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
391     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
392     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
393     }
394    
395    
396     }
397     }
398     totalNECount+=NE1*NE0;
399     faceNECount+=NE1*NE0;
400     }
401     if (!periodic[0]) {
402     /* ** elements on boundary 001 (x1=0): */
403    
404     #pragma omp parallel for private(i1,i2,k,node0)
405     for (i2=0;i2<NE2;i2++) {
406     for (i1=0;i1<NE1;i1++) {
407     k=i1+NE1*i2+faceNECount;
408     node0=i1*N0+N0*N1*i2;
409    
410     out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
411     out->FaceElements->Tag[k]=1;
412     out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
413 bcumming 782 #ifdef PASO_MPI
414     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
415     #endif
416 jgs 82
417     if (useElementsOnFace) {
418     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
419     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
420     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
421     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
422     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
423     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+1;
424     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
425     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+1;
426     } else {
427     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
428     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
429     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
430     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
431     }
432     }
433     }
434     totalNECount+=NE1*NE2;
435     faceNECount+=NE1*NE2;
436    
437     /* ** elements on boundary 002 (x1=1): */
438    
439     #pragma omp parallel for private(i1,i2,k,node0)
440     for (i2=0;i2<NE2;i2++) {
441     for (i1=0;i1<NE1;i1++) {
442     k=i1+NE1*i2+faceNECount;
443     node0=(NE0-1)+i1*N0+N0*N1*i2 ;
444    
445     out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
446     out->FaceElements->Tag[k]=2;
447     out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
448 bcumming 782 #ifdef PASO_MPI
449     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
450     #endif
451 jgs 82
452     if (useElementsOnFace) {
453     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
454     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
455     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
456     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
457     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
458     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;
459     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0;
460     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1;
461     } else {
462     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
463     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
464     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
465     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
466     }
467     }
468     }
469     totalNECount+=NE1*NE2;
470     faceNECount+=NE1*NE2;
471     }
472     if (!periodic[1]) {
473     /* ** elements on boundary 010 (x2=0): */
474    
475     #pragma omp parallel for private(i0,i2,k,node0)
476     for (i2=0;i2<NE2;i2++) {
477     for (i0=0;i0<NE0;i0++) {
478     k=i0+NE0*i2+faceNECount;
479     node0=i0+N0*N1*i2;
480    
481     out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
482     out->FaceElements->Tag[k]=10;
483     out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
484 bcumming 782 #ifdef PASO_MPI
485     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
486     #endif
487 jgs 82
488     if (useElementsOnFace) {
489     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
490     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
491     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
492     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
493     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;
494     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+1;
495     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*N0+N0+1;
496     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0+N0;
497     } else {
498     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
499     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
500     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
501     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
502     }
503     }
504     }
505     totalNECount+=NE0*NE2;
506     faceNECount+=NE0*NE2;
507    
508     /* ** elements on boundary 020 (x2=1): */
509    
510     #pragma omp parallel for private(i0,i2,k,node0)
511     for (i2=0;i2<NE2;i2++) {
512     for (i0=0;i0<NE0;i0++) {
513     k=i0+NE0*i2+faceNECount;
514     node0=i0+(NE1-1)*N0+N0*N1*i2;
515    
516     out->FaceElements->Tag[k]=20;
517     out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
518     out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
519 bcumming 782 #ifdef PASO_MPI
520     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
521     #endif
522 jgs 82
523     if (useElementsOnFace) {
524     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
525     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
526     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
527     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
528     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
529     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1;
530     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+1;
531     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
532     } else {
533     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
534     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
535     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
536     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
537     }
538    
539     }
540     }
541     totalNECount+=NE0*NE2;
542     faceNECount+=NE0*NE2;
543     }
544 jgs 123 out->FaceElements->minColor=0;
545     out->FaceElements->maxColor=23;
546 jgs 82
547 bcumming 751 #ifdef PASO_MPI
548 bcumming 782 Finley_ElementFile_setDomainFlags( out->Elements );
549     Finley_ElementFile_setDomainFlags( out->FaceElements );
550     Finley_ElementFile_setDomainFlags( out->ContactElements );
551     Finley_ElementFile_setDomainFlags( out->Points );
552    
553     /* reorder the degrees of freedom */
554     Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
555     #endif
556 bcumming 751
557 jgs 82 /* condense the nodes: */
558     Finley_Mesh_resolveNodeIds(out);
559    
560     /* prepare mesh for further calculatuions:*/
561     Finley_Mesh_prepare(out) ;
562    
563 jgs 149 #ifdef Finley_TRACE
564 jgs 82 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
565 jgs 149 #endif
566 gross 964 if (Finley_noError()) {
567     if ( ! Finley_Mesh_isPrepared(out)) {
568     Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");
569     }
570     } else {
571 jgs 82 Finley_Mesh_dealloc(out);
572     }
573     return out;
574     }
575 bcumming 751
576     #ifdef PASO_MPI
577 gross 1062 Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order, index_t reduced_order, bool_t useElementsOnFace)
578 bcumming 751 {
579 bcumming 782 dim_t N0,N1,N2,N0t,NDOF0t,NE0,NE1,NE2,i0,i1,i2,kk,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES;//,M0,M1,M2;
580 bcumming 751 dim_t idCount, NE0_local, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal;
581     bool_t dom_left, dom_right, dom_internal;
582 bcumming 782 index_t firstNode=0, DOFcount=0, node0, node1, node2;
583     index_t targetDomain=-1, firstNodeConstruct, j;
584     bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft=FALSE, boundaryRight=FALSE;
585     index_t *indexBackward=NULL, *indexForward=NULL,*facePerm=NULL, *forwardDOF=NULL, *backwardDOF=NULL;
586 bcumming 751 Finley_Mesh* out;
587    
588     char name[50];
589     Paso_MPIInfo *mpi_info = NULL;
590     double time0=Finley_timer();
591    
592 bcumming 782 index_t face0[] = {3, 0, 4, 7, 2, 1, 5, 6};
593     index_t face1[] = {1, 2, 6, 5, 0, 3, 7, 4};
594     index_t face2[] = {0, 3, 2, 1, 4, 7, 6, 5};
595     index_t face3[] = {4, 5, 6, 7, 0, 1, 2, 3};
596     index_t face4[] = {0, 1, 5, 4, 3, 2, 6, 7};
597     index_t face5[] = {3, 7, 6, 2, 0, 4, 5, 1};
598 bcumming 751 NE0=MAX(1,numElements[0]);
599     NE1=MAX(1,numElements[1]);
600     NE2=MAX(1,numElements[2]);
601     N0=NE0+1;
602     N1=NE1+1;
603     N2=NE2+1;
604    
605    
606     /* get MPI information */
607     mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
608     if (! Finley_noError())
609     return NULL;
610    
611     /* use the serial version to generate the mesh for the 1-CPU case */
612     if( mpi_info->size==1 )
613     {
614 gross 1062 out = Finley_RectangularMesh_Hex8_singleCPU( numElements, Length, periodic, order, reduced_order, useElementsOnFace, mpi_info );
615 bcumming 751 return out;
616     }
617    
618     if( mpi_info->rank==0 )
619     domLeft = TRUE;
620     if( mpi_info->rank==mpi_info->size-1 )
621     domRight = TRUE;
622     if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
623     domInternal = TRUE;
624    
625     /* dimensions of the local subdomain */
626     domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal );
627    
628     /* count Degrees of Freedom along each axis */
629     NDOF0 = N0 - periodic[0];
630     NDOF1 = N1 - periodic[1];
631     NDOF2 = N2 - periodic[2];
632    
633     /* count face elements */
634     /* internal face elements */
635     NFaceElements = 0;
636     if( !periodic[0] )
637     NFaceElements += (domLeft+domRight)*NE1*NE2;
638     if( !periodic[1] )
639 bcumming 782 NFaceElements += 2*numElementsLocal*NE2;
640 bcumming 751 if( !periodic[2] )
641 bcumming 782 NFaceElements += 2*numElementsLocal*NE1;
642    
643     boundaryLeft = !domLeft || periodicLocal[0];
644     boundaryRight = !domRight || periodicLocal[1];
645     N0t = numNodesLocal + boundaryRight + boundaryLeft;
646     NDOF0t = numDOFLocal + boundaryRight + boundaryLeft;
647     firstNodeConstruct = firstNode - boundaryLeft;
648     firstNodeConstruct = firstNodeConstruct<0 ? N0-2 : firstNodeConstruct;
649 bcumming 751
650     /* allocate mesh: */
651     sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
652    
653 gross 1062 out=Finley_Mesh_alloc(name,3,order,reduced_order,mpi_info);
654 bcumming 751 if (! Finley_noError()) return NULL;
655    
656 gross 1062 out->Elements=Finley_ElementFile_alloc(Hex8,out->order,out->reduced_order, mpi_info);
657 bcumming 751 if (useElementsOnFace) {
658 gross 1062 out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,out->reduced_order, mpi_info);
659     out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,out->reduced_order, mpi_info);
660 bcumming 751 } else {
661 gross 1062 out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,out->reduced_order, mpi_info);
662     out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,out->reduced_order, mpi_info);
663 bcumming 751 }
664 gross 1062 out->Points=Finley_ElementFile_alloc(Point1,out->order,out->reduced_order, mpi_info);
665 bcumming 751
666     if (! Finley_noError()) {
667     Finley_Mesh_dealloc(out);
668     return NULL;
669     }
670    
671     /* allocate tables: */
672 bcumming 782 Finley_NodeFile_allocTable(out->Nodes,N0t*N1*N2);
673 bcumming 751 Finley_ElementFile_allocTable(out->Elements,(numElementsLocal)*NE1*NE2);
674     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
675 bcumming 788
676     Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal*NDOF1*NDOF2, NDOF1*NDOF2*2, 0 );
677     Finley_ElementDistribution_allocTable( out->Elements->elementDistribution, numElementsLocal*NE1*NE2, NE1*NE2*(numElementsLocal-boundaryRight*(!periodic[1])) );
678     Finley_ElementDistribution_allocTable( out->FaceElements->elementDistribution, NFaceElements, NFaceElements-2*boundaryRight*(NE2*(!periodic[1])+NE1*(!periodic[2])) );
679 bcumming 751 if (! Finley_noError()) {
680     Finley_Mesh_dealloc(out);
681     return NULL;
682     }
683    
684     /* set nodes: */
685     /* INTERNAL/BOUNDARY NODES */
686 bcumming 759 k=0;
687 bcumming 751 #pragma omp parallel for private(i0,i1,i2,k)
688 bcumming 759 for (i2=0;i2<N2;i2++) {
689 bcumming 751 for (i1=0;i1<N1;i1++) {
690 bcumming 782 for (i0=0;i0<N0t;i0++,k++) {
691     out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE((i0+firstNodeConstruct) % N0)/DBLE(N0-1)*Length[0];
692 bcumming 751 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
693     out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
694     out->Nodes->Id[k]=k;
695     out->Nodes->Tag[k]=0;
696 bcumming 782 out->Nodes->degreeOfFreedom[k]=i0 + (i1%NDOF1)*N0t + (i2%NDOF2)*N0t*N1;
697     out->Nodes->Dom[k]=NODE_INTERNAL;
698 bcumming 751 }
699     }
700     }
701 bcumming 782
702     /* mark the nodes that reference external and boundary DOF as such */
703     if( boundaryLeft ){
704     for( i1=0; i1<N1; i1++ )
705     for( i2=0; i2<N2; i2++ ) {
706     out->Nodes->Dom[N1*N0t*i2+N0t*i1] = NODE_EXTERNAL;
707     out->Nodes->Dom[N1*N0t*i2+N0t*i1+1] = NODE_BOUNDARY;
708     }
709     }
710     if( boundaryRight ){
711     for( i1=0; i1<N1; i1++ )
712     for( i2=0; i2<N2; i2++ ) {
713     out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-1] = NODE_EXTERNAL;
714     out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-2] = NODE_BOUNDARY;
715     }
716     }
717     if( periodicLocal[0] ){
718     for( i1=0; i1<N1; i1++ )
719     for( i2=0; i2<N2; i2++ ) {
720     out->Nodes->degreeOfFreedom[N1*N0t*i2+i1*N0t+3] = out->Nodes->degreeOfFreedom[N1*N0t*i2+i1*N0t+2];
721     out->Nodes->Dom[N1*N0t*i2+N0t*i1+3] = NODE_BOUNDARY;
722     }
723     }
724    
725     /* tag Nodes that are referenced by face elements */
726 bcumming 751 if (!periodic[2]) {
727     for (i1=0;i1<N1;i1++) {
728 bcumming 782 for (i0=0;i0<N0t;i0++) {
729     out->Nodes->Tag[i0 + N0t*i1]+=100;
730     out->Nodes->Tag[i0 + N0t*i1 + N0t*N1*(N2-1)]+=200;
731 bcumming 751 }
732     }
733     }
734     if (!periodic[1]) {
735     for (i2=0;i2<N2;i2++) {
736 bcumming 782 for (i0=0;i0<N0t;i0++) {
737     out->Nodes->Tag[i0 + i2*N1*N0t]+=10;
738     out->Nodes->Tag[i0 + (i2+1)*N1*N0t-N0t]+=20;
739 bcumming 751 }
740     }
741     }
742     if (!periodic[0] && !domInternal ) {
743     for (i2=0;i2<N2;i2++) {
744     for (i1=0;i1<N1;i1++) {
745     if( domLeft )
746 bcumming 782 out->Nodes->Tag[i1*N0t + i2*N0t*N1]+=1;
747 bcumming 751 if( domRight )
748 bcumming 782 out->Nodes->Tag[(i1+1)*N0t-1 + i2*N0t*N1]+=2;
749 bcumming 751 }
750     }
751     }
752 bcumming 782
753     /* form the boudary communication information */
754     forwardDOF = MEMALLOC(NDOF1*NDOF2,index_t);
755     backwardDOF = MEMALLOC(NDOF1*NDOF2,index_t);
756     if( !(mpi_info->size==2 && periodicLocal[0])){
757     if( boundaryLeft ) {
758     targetDomain = mpi_info->rank-1 < 0 ? mpi_info->size-1 : mpi_info->rank-1;
759     for( i2=0; i2<NDOF2; i2++ ){
760     for( i1=0; i1<NDOF1; i1++ ){
761     forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];
762     backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];
763     }
764     }
765     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
766     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
767 bcumming 751 }
768 bcumming 782 if( boundaryRight ) {
769     targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;
770     for( i2=0; i2<NDOF2; i2++ ){
771     for( i1=0; i1<NDOF1; i1++ ){
772     forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];
773     backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];
774     }
775     }
776     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
777     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
778 bcumming 751 }
779 bcumming 782 } else{
780     /* periodic boundary conditions with 2 domains, need to change the order in which domain 0 passes boundary data */
781     targetDomain = 1;
782    
783     for( i2=0; i2<NDOF2; i2++ ){
784     for( i1=0; i1<NDOF1; i1++ ){
785     forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];
786     backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];
787     }
788     }
789     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
790     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
791 bcumming 751
792 bcumming 782 for( i2=0; i2<NDOF2; i2++ ){
793     for( i1=0; i1<NDOF1; i1++ ){
794     forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];
795     backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];
796     }
797     }
798     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
799     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
800 bcumming 751 }
801 bcumming 782 MEMFREE( forwardDOF );
802     MEMFREE( backwardDOF );
803 bcumming 751 /* set the elements: */
804    
805     /* INTERNAL elements */
806     k = 0;
807     #pragma omp parallel for private(i0,i1,i2,k,node0)
808     for (i2=0;i2<NE2;i2++) {
809     for (i1=0;i1<NE1;i1++) {
810 bcumming 782 for (i0=0;i0<numElementsLocal;i0++,k++) {
811     node0 = (periodicLocal[0] && !i0) ? i1*N0t + i2*N1*N0t : i1*N0t + i2*N1*N0t + i0 + periodicLocal[0];
812 bcumming 751
813 bcumming 787 out->Elements->Id[k]=((firstNodeConstruct+i0)%NE0)*NE1*NE2 + NE1*i2 + i1;
814 bcumming 751 out->Elements->Tag[k]=0;
815 bcumming 782 out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
816     out->Elements->Dom[k]=ELEMENT_INTERNAL;
817 bcumming 751
818     out->Elements->Nodes[INDEX2(0,k,8)]=node0;
819     out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
820 bcumming 782 out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0t+1;
821     out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0t;
822     out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0t*N1;
823     out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0t*N1+1;
824     out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0t*N1+N0t+1;
825     out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0t*N1+N0t;
826 bcumming 751
827     }
828     }
829     }
830     out->Elements->minColor=0;
831 bcumming 782 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
832     if( boundaryLeft )
833     for( i2=0; i2<NE2; i2++ )
834     for( i1=0; i1<NE1; i1++ )
835     out->Elements->Dom[i2*NE1*numElementsLocal+i1*numElementsLocal]=ELEMENT_BOUNDARY;
836     if( boundaryRight )
837     for( i2=0; i2<NE2; i2++ )
838     for( i1=0; i1<NE1; i1++ )
839     out->Elements->Dom[i2*NE1*numElementsLocal+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
840    
841     Finley_ElementFile_setDomainFlags( out->Elements );
842 bcumming 751
843     /* face elements: */
844     if (useElementsOnFace) {
845     NUMNODES=8;
846     } else {
847     NUMNODES=4;
848     }
849 bcumming 782 totalNECount=out->Elements->numElements;
850 bcumming 751 faceNECount=0;
851     idCount = totalNECount;
852    
853     /* these are the quadrilateral elements on boundary 1 (x3=0): */
854     numElementsInternal = numElementsLocal-nodesExternal[0]-nodesExternal[1];
855     if (!periodic[2]) {
856     /* elements on boundary 100 (x3=0): */
857    
858 bcumming 782 #pragma omp parallel for private(i0,i1,k)
859 bcumming 751 for (i1=0;i1<NE1;i1++) {
860 bcumming 782 for (i0=0; i0<numElementsLocal; i0++) {
861     k=i0+numElementsLocal*i1+faceNECount;
862     kk=i0 + i1*numElementsLocal;
863     facePerm = face2;
864 bcumming 751
865     out->FaceElements->Id[k]=idCount++;
866     out->FaceElements->Tag[k]=100;
867 bcumming 782 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
868     out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
869 bcumming 751
870 bcumming 782 for( j=0; j<NUMNODES; j++ )
871     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
872 bcumming 751 }
873     }
874 bcumming 782 if( boundaryLeft ){
875     for( i1=0; i1<NE1; i1++ )
876     out->FaceElements->Dom[faceNECount+i1*numElementsLocal]=ELEMENT_BOUNDARY;
877     if( periodicLocal[0] )
878     for( i1=0; i1<NE1; i1++ )
879     out->FaceElements->Dom[faceNECount+i1*numElementsLocal+1]=ELEMENT_BOUNDARY;
880     }
881     if( boundaryRight )
882     for( i1=0; i1<NE1; i1++ )
883     out->FaceElements->Dom[faceNECount+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
884     totalNECount+=NE1*numElementsLocal;
885     faceNECount+=NE1*numElementsLocal;
886 bcumming 751
887     /* ** elements on boundary 200 (x3=1) */
888    
889 bcumming 782 #pragma omp parallel for private(i0,i1,k)
890 bcumming 751 for (i1=0;i1<NE1;i1++) {
891 bcumming 782 for (i0=0;i0<numElementsLocal;i0++) {
892     k=i0+numElementsLocal*i1+faceNECount;
893     kk=i0+i1*numElementsLocal+numElementsLocal*NE1*(NE2-1);
894     facePerm = face3;
895 bcumming 751
896     out->FaceElements->Id[k]=idCount++;
897     out->FaceElements->Tag[k]=200;
898 bcumming 782 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
899     out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
900 bcumming 751
901 bcumming 782 for( j=0; j<NUMNODES; j++ )
902     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
903 bcumming 751 }
904     }
905 bcumming 782 if( boundaryLeft ){
906     for( i1=0; i1<NE1; i1++ )
907     out->FaceElements->Dom[faceNECount+i1*numElementsLocal]=ELEMENT_BOUNDARY;
908     if( periodicLocal[0] )
909     for( i1=0; i1<NE1; i1++ )
910     out->FaceElements->Dom[faceNECount+i1*numElementsLocal+1]=ELEMENT_BOUNDARY;
911     }
912     if( boundaryRight )
913     for( i1=0; i1<NE1; i1++ )
914     out->FaceElements->Dom[faceNECount+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
915     totalNECount+=NE1*numElementsLocal;
916     faceNECount+=NE1*numElementsLocal;
917 bcumming 751 }
918     if (!periodic[0] && !domInternal) {
919     /* ** elements on boundary 001 (x1=0): */
920     if( domLeft ){
921 bcumming 782 #pragma omp parallel for private(i1,i2,k)
922 bcumming 751 for (i2=0;i2<NE2;i2++) {
923     for (i1=0;i1<NE1;i1++) {
924     k=i1+NE1*i2+faceNECount;
925 bcumming 782 kk=i1*numElementsLocal + i2*numElementsLocal*NE1;
926     facePerm = face0;
927 bcumming 751
928     out->FaceElements->Id[k]=idCount++;
929     out->FaceElements->Tag[k]=1;
930 bcumming 782 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
931     out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
932 bcumming 751
933 bcumming 782 for( j=0; j<NUMNODES; j++ )
934     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
935 bcumming 751 }
936     }
937     totalNECount+=NE1*NE2;
938     faceNECount+=NE1*NE2;
939     }
940     /* ** elements on boundary 002 (x1=1): */
941     if( domRight ) {
942 bcumming 782 #pragma omp parallel for private(i1,i2,k)
943 bcumming 751 for (i2=0;i2<NE2;i2++) {
944     for (i1=0;i1<NE1;i1++) {
945     k=i1+NE1*i2+faceNECount;
946 bcumming 782 kk=(i1+1)*numElementsLocal + i2*numElementsLocal*NE1 - 1;
947     facePerm = face1;
948 bcumming 751
949     out->FaceElements->Id[k]=idCount++;
950     out->FaceElements->Tag[k]=2;
951 bcumming 782 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
952     out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
953 bcumming 751
954 bcumming 782 for( j=0; j<NUMNODES; j++ )
955     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
956 bcumming 751 }
957     }
958     totalNECount+=NE1*NE2;
959     faceNECount+=NE1*NE2;
960     }
961     }
962     if (!periodic[1]) {
963     /* ** elements on boundary 010 (x2=0): */
964    
965 bcumming 782 #pragma omp parallel for private(i0,i2,k)
966 bcumming 751 for (i2=0;i2<NE2;i2++) {
967 bcumming 782 for (i0=0;i0<numElementsLocal;i0++) {
968     k=i0+numElementsLocal*i2+faceNECount;
969     kk=i0+numElementsLocal*NE1*i2;
970     facePerm = face4;
971 bcumming 751
972     out->FaceElements->Id[k]=idCount++;
973     out->FaceElements->Tag[k]=10;
974 bcumming 782 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
975     out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
976 bcumming 751
977 bcumming 782 for( j=0; j<NUMNODES; j++ )
978     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
979 bcumming 751 }
980     }
981 bcumming 782 if( boundaryLeft ){
982     for( i2=0; i2<NE2; i2++ )
983     out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;
984     if( periodicLocal[0] )
985     for( i2=0; i2<NE2; i2++ )
986     out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;
987     }
988     if( boundaryRight )
989     for( i2=0; i2<NE2; i2++ )
990     out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
991     totalNECount+=numElementsLocal*NE2;
992     faceNECount+=numElementsLocal*NE2;
993 bcumming 751
994     /* ** elements on boundary 020 (x2=1): */
995    
996 bcumming 782 #pragma omp parallel for private(i0,i2,k)
997 bcumming 751 for (i2=0;i2<NE2;i2++) {
998 bcumming 782 for (i0=0;i0<numElementsLocal;i0++) {
999     k=i0+numElementsLocal*i2+faceNECount;
1000     kk=i0+numElementsLocal*NE1*(i2+1)-numElementsLocal;
1001     facePerm = face5;
1002 bcumming 751
1003     out->FaceElements->Tag[k]=20;
1004     out->FaceElements->Id[k]=idCount++;
1005 bcumming 782 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
1006 bcumming 751 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
1007    
1008 bcumming 782 for( j=0; j<NUMNODES; j++ )
1009     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
1010 bcumming 751 }
1011     }
1012 bcumming 782 if( boundaryLeft ){
1013     for( i2=0; i2<NE2; i2++ )
1014     out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;
1015     if( periodicLocal[0] )
1016     for( i2=0; i2<NE2; i2++ )
1017     out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;
1018     }
1019     if( boundaryRight )
1020     for( i2=0; i2<NE2; i2++ )
1021     out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
1022     totalNECount+=numElementsLocal*NE2;
1023     faceNECount+=numElementsLocal*NE2;
1024 bcumming 751 }
1025     out->FaceElements->elementDistribution->numInternal = faceNECount;
1026    
1027     out->FaceElements->minColor=0;
1028 bcumming 782 out->FaceElements->maxColor=23;
1029     out->FaceElements->numElements=faceNECount;
1030    
1031     Finley_ElementFile_setDomainFlags( out->FaceElements );
1032 bcumming 751
1033 bcumming 782 /* setup distribution info for other elements */
1034     Finley_ElementFile_setDomainFlags( out->ContactElements );
1035     Finley_ElementFile_setDomainFlags( out->Points );
1036 bcumming 751
1037 bcumming 782 /* reorder the degrees of freedom */
1038     Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
1039 bcumming 751
1040     /* condense the nodes: */
1041 bcumming 782 Finley_Mesh_resolveNodeIds(out);
1042     if( !Finley_MPI_noError(mpi_info) )
1043     {
1044 bcumming 751 Paso_MPIInfo_dealloc( mpi_info );
1045     Finley_Mesh_dealloc(out);
1046     return NULL;
1047 bcumming 782 }
1048 bcumming 751
1049     /* prepare mesh for further calculatuions:*/
1050 bcumming 782 Finley_Mesh_prepare(out);
1051     if( !Finley_MPI_noError(mpi_info) )
1052     {
1053     Paso_MPIInfo_dealloc( mpi_info );
1054     Finley_Mesh_dealloc(out);
1055     return NULL;
1056     }
1057 bcumming 751
1058 bcumming 782 /* free up memory */
1059     Paso_MPIInfo_dealloc( mpi_info );
1060 bcumming 751
1061     #ifdef Finley_TRACE
1062     printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
1063     #endif
1064 gross 964 if (Finley_noError()) {
1065     if ( ! Finley_Mesh_isPrepared(out) ) {
1066     Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");
1067     }
1068     }
1069 bcumming 782 return out;
1070 bcumming 751 }
1071     #endif
1072    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26