/[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 788 - (hide annotations)
Wed Jul 26 05:12:15 2006 UTC (13 years ago) by bcumming
File MIME type: text/plain
File size: 37216 byte(s)
small corrections to code in previous commit, and some header files that
should have been included in the last commit.


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     Finley_Mesh* Finley_RectangularMesh_Hex8_singleCPU(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace,Paso_MPIInfo *mpi_info)
119     #else
120 bcumming 751 Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
121     #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 jgs 82 out=Finley_Mesh_alloc(name,3,order);
194 bcumming 730 #else
195 bcumming 751 out=Finley_Mesh_alloc(name,3,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     out->Elements=Finley_ElementFile_alloc(Hex8,out->order,mpi_info);
201     if (useElementsOnFace) {
202     out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,mpi_info);
203     out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,mpi_info);
204     } else {
205     out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,mpi_info);
206     out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,mpi_info);
207     }
208     out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
209     #else
210 jgs 82 out->Elements=Finley_ElementFile_alloc(Hex8,out->order);
211     if (useElementsOnFace) {
212     out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order);
213     out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order);
214     } else {
215     out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order);
216     out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order);
217     }
218     out->Points=Finley_ElementFile_alloc(Point1,out->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 jgs 82
567 jgs 150 if (! Finley_noError()) {
568 jgs 82 Finley_Mesh_dealloc(out);
569     return NULL;
570     }
571     return out;
572     }
573 bcumming 751
574     #ifdef PASO_MPI
575     Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
576     {
577 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;
578 bcumming 751 dim_t idCount, NE0_local, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal;
579     bool_t dom_left, dom_right, dom_internal;
580 bcumming 782 index_t firstNode=0, DOFcount=0, node0, node1, node2;
581     index_t targetDomain=-1, firstNodeConstruct, j;
582     bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft=FALSE, boundaryRight=FALSE;
583     index_t *indexBackward=NULL, *indexForward=NULL,*facePerm=NULL, *forwardDOF=NULL, *backwardDOF=NULL;
584 bcumming 751 Finley_Mesh* out;
585    
586     char name[50];
587     Paso_MPIInfo *mpi_info = NULL;
588     double time0=Finley_timer();
589    
590 bcumming 782 index_t face0[] = {3, 0, 4, 7, 2, 1, 5, 6};
591     index_t face1[] = {1, 2, 6, 5, 0, 3, 7, 4};
592     index_t face2[] = {0, 3, 2, 1, 4, 7, 6, 5};
593     index_t face3[] = {4, 5, 6, 7, 0, 1, 2, 3};
594     index_t face4[] = {0, 1, 5, 4, 3, 2, 6, 7};
595     index_t face5[] = {3, 7, 6, 2, 0, 4, 5, 1};
596 bcumming 751 NE0=MAX(1,numElements[0]);
597     NE1=MAX(1,numElements[1]);
598     NE2=MAX(1,numElements[2]);
599     N0=NE0+1;
600     N1=NE1+1;
601     N2=NE2+1;
602    
603    
604     /* get MPI information */
605     mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
606     if (! Finley_noError())
607     return NULL;
608    
609     /* use the serial version to generate the mesh for the 1-CPU case */
610     if( mpi_info->size==1 )
611     {
612 bcumming 782 out = Finley_RectangularMesh_Hex8_singleCPU( numElements, Length, periodic, order, useElementsOnFace, mpi_info );
613 bcumming 751 return out;
614     }
615    
616     if( mpi_info->rank==0 )
617     domLeft = TRUE;
618     if( mpi_info->rank==mpi_info->size-1 )
619     domRight = TRUE;
620     if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
621     domInternal = TRUE;
622    
623     /* dimensions of the local subdomain */
624     domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal );
625    
626     /* count Degrees of Freedom along each axis */
627     NDOF0 = N0 - periodic[0];
628     NDOF1 = N1 - periodic[1];
629     NDOF2 = N2 - periodic[2];
630    
631     /* count face elements */
632     /* internal face elements */
633     NFaceElements = 0;
634     if( !periodic[0] )
635     NFaceElements += (domLeft+domRight)*NE1*NE2;
636     if( !periodic[1] )
637 bcumming 782 NFaceElements += 2*numElementsLocal*NE2;
638 bcumming 751 if( !periodic[2] )
639 bcumming 782 NFaceElements += 2*numElementsLocal*NE1;
640    
641     boundaryLeft = !domLeft || periodicLocal[0];
642     boundaryRight = !domRight || periodicLocal[1];
643     N0t = numNodesLocal + boundaryRight + boundaryLeft;
644     NDOF0t = numDOFLocal + boundaryRight + boundaryLeft;
645     firstNodeConstruct = firstNode - boundaryLeft;
646     firstNodeConstruct = firstNodeConstruct<0 ? N0-2 : firstNodeConstruct;
647 bcumming 751
648     /* allocate mesh: */
649     sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
650    
651     out=Finley_Mesh_alloc(name,3,order,mpi_info);
652     if (! Finley_noError()) return NULL;
653    
654     out->Elements=Finley_ElementFile_alloc(Hex8,out->order,mpi_info);
655     if (useElementsOnFace) {
656     out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,mpi_info);
657     out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,mpi_info);
658     } else {
659     out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,mpi_info);
660     out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,mpi_info);
661     }
662     out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
663    
664     if (! Finley_noError()) {
665     Finley_Mesh_dealloc(out);
666     return NULL;
667     }
668    
669     /* allocate tables: */
670 bcumming 782 Finley_NodeFile_allocTable(out->Nodes,N0t*N1*N2);
671 bcumming 751 Finley_ElementFile_allocTable(out->Elements,(numElementsLocal)*NE1*NE2);
672     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
673 bcumming 788
674     Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal*NDOF1*NDOF2, NDOF1*NDOF2*2, 0 );
675     Finley_ElementDistribution_allocTable( out->Elements->elementDistribution, numElementsLocal*NE1*NE2, NE1*NE2*(numElementsLocal-boundaryRight*(!periodic[1])) );
676     Finley_ElementDistribution_allocTable( out->FaceElements->elementDistribution, NFaceElements, NFaceElements-2*boundaryRight*(NE2*(!periodic[1])+NE1*(!periodic[2])) );
677 bcumming 751 if (! Finley_noError()) {
678     Finley_Mesh_dealloc(out);
679     return NULL;
680     }
681    
682     /* set nodes: */
683     /* INTERNAL/BOUNDARY NODES */
684 bcumming 759 k=0;
685 bcumming 751 #pragma omp parallel for private(i0,i1,i2,k)
686 bcumming 759 for (i2=0;i2<N2;i2++) {
687 bcumming 751 for (i1=0;i1<N1;i1++) {
688 bcumming 782 for (i0=0;i0<N0t;i0++,k++) {
689     out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE((i0+firstNodeConstruct) % N0)/DBLE(N0-1)*Length[0];
690 bcumming 751 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
691     out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
692     out->Nodes->Id[k]=k;
693     out->Nodes->Tag[k]=0;
694 bcumming 782 out->Nodes->degreeOfFreedom[k]=i0 + (i1%NDOF1)*N0t + (i2%NDOF2)*N0t*N1;
695     out->Nodes->Dom[k]=NODE_INTERNAL;
696 bcumming 751 }
697     }
698     }
699 bcumming 782
700     /* mark the nodes that reference external and boundary DOF as such */
701     if( boundaryLeft ){
702     for( i1=0; i1<N1; i1++ )
703     for( i2=0; i2<N2; i2++ ) {
704     out->Nodes->Dom[N1*N0t*i2+N0t*i1] = NODE_EXTERNAL;
705     out->Nodes->Dom[N1*N0t*i2+N0t*i1+1] = NODE_BOUNDARY;
706     }
707     }
708     if( boundaryRight ){
709     for( i1=0; i1<N1; i1++ )
710     for( i2=0; i2<N2; i2++ ) {
711     out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-1] = NODE_EXTERNAL;
712     out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-2] = NODE_BOUNDARY;
713     }
714     }
715     if( periodicLocal[0] ){
716     for( i1=0; i1<N1; i1++ )
717     for( i2=0; i2<N2; i2++ ) {
718     out->Nodes->degreeOfFreedom[N1*N0t*i2+i1*N0t+3] = out->Nodes->degreeOfFreedom[N1*N0t*i2+i1*N0t+2];
719     out->Nodes->Dom[N1*N0t*i2+N0t*i1+3] = NODE_BOUNDARY;
720     }
721     }
722    
723     /* tag Nodes that are referenced by face elements */
724 bcumming 751 if (!periodic[2]) {
725     for (i1=0;i1<N1;i1++) {
726 bcumming 782 for (i0=0;i0<N0t;i0++) {
727     out->Nodes->Tag[i0 + N0t*i1]+=100;
728     out->Nodes->Tag[i0 + N0t*i1 + N0t*N1*(N2-1)]+=200;
729 bcumming 751 }
730     }
731     }
732     if (!periodic[1]) {
733     for (i2=0;i2<N2;i2++) {
734 bcumming 782 for (i0=0;i0<N0t;i0++) {
735     out->Nodes->Tag[i0 + i2*N1*N0t]+=10;
736     out->Nodes->Tag[i0 + (i2+1)*N1*N0t-N0t]+=20;
737 bcumming 751 }
738     }
739     }
740     if (!periodic[0] && !domInternal ) {
741     for (i2=0;i2<N2;i2++) {
742     for (i1=0;i1<N1;i1++) {
743     if( domLeft )
744 bcumming 782 out->Nodes->Tag[i1*N0t + i2*N0t*N1]+=1;
745 bcumming 751 if( domRight )
746 bcumming 782 out->Nodes->Tag[(i1+1)*N0t-1 + i2*N0t*N1]+=2;
747 bcumming 751 }
748     }
749     }
750 bcumming 782
751     /* form the boudary communication information */
752     forwardDOF = MEMALLOC(NDOF1*NDOF2,index_t);
753     backwardDOF = MEMALLOC(NDOF1*NDOF2,index_t);
754     if( !(mpi_info->size==2 && periodicLocal[0])){
755     if( boundaryLeft ) {
756     targetDomain = mpi_info->rank-1 < 0 ? mpi_info->size-1 : mpi_info->rank-1;
757     for( i2=0; i2<NDOF2; i2++ ){
758     for( i1=0; i1<NDOF1; i1++ ){
759     forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];
760     backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];
761     }
762     }
763     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
764     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
765 bcumming 751 }
766 bcumming 782 if( boundaryRight ) {
767     targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;
768     for( i2=0; i2<NDOF2; i2++ ){
769     for( i1=0; i1<NDOF1; i1++ ){
770     forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];
771     backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];
772     }
773     }
774     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
775     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
776 bcumming 751 }
777 bcumming 782 } else{
778     /* periodic boundary conditions with 2 domains, need to change the order in which domain 0 passes boundary data */
779     targetDomain = 1;
780    
781     for( i2=0; i2<NDOF2; i2++ ){
782     for( i1=0; i1<NDOF1; i1++ ){
783     forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];
784     backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];
785     }
786     }
787     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
788     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
789 bcumming 751
790 bcumming 782 for( i2=0; i2<NDOF2; i2++ ){
791     for( i1=0; i1<NDOF1; i1++ ){
792     forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];
793     backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];
794     }
795     }
796     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
797     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
798 bcumming 751 }
799 bcumming 782 MEMFREE( forwardDOF );
800     MEMFREE( backwardDOF );
801 bcumming 751 /* set the elements: */
802    
803     /* INTERNAL elements */
804     k = 0;
805     #pragma omp parallel for private(i0,i1,i2,k,node0)
806     for (i2=0;i2<NE2;i2++) {
807     for (i1=0;i1<NE1;i1++) {
808 bcumming 782 for (i0=0;i0<numElementsLocal;i0++,k++) {
809     node0 = (periodicLocal[0] && !i0) ? i1*N0t + i2*N1*N0t : i1*N0t + i2*N1*N0t + i0 + periodicLocal[0];
810 bcumming 751
811 bcumming 787 out->Elements->Id[k]=((firstNodeConstruct+i0)%NE0)*NE1*NE2 + NE1*i2 + i1;
812 bcumming 751 out->Elements->Tag[k]=0;
813 bcumming 782 out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
814     out->Elements->Dom[k]=ELEMENT_INTERNAL;
815 bcumming 751
816     out->Elements->Nodes[INDEX2(0,k,8)]=node0;
817     out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
818 bcumming 782 out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0t+1;
819     out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0t;
820     out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0t*N1;
821     out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0t*N1+1;
822     out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0t*N1+N0t+1;
823     out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0t*N1+N0t;
824 bcumming 751
825     }
826     }
827     }
828     out->Elements->minColor=0;
829 bcumming 782 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
830     if( boundaryLeft )
831     for( i2=0; i2<NE2; i2++ )
832     for( i1=0; i1<NE1; i1++ )
833     out->Elements->Dom[i2*NE1*numElementsLocal+i1*numElementsLocal]=ELEMENT_BOUNDARY;
834     if( boundaryRight )
835     for( i2=0; i2<NE2; i2++ )
836     for( i1=0; i1<NE1; i1++ )
837     out->Elements->Dom[i2*NE1*numElementsLocal+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
838    
839     Finley_ElementFile_setDomainFlags( out->Elements );
840 bcumming 751
841     /* face elements: */
842     if (useElementsOnFace) {
843     NUMNODES=8;
844     } else {
845     NUMNODES=4;
846     }
847 bcumming 782 totalNECount=out->Elements->numElements;
848 bcumming 751 faceNECount=0;
849     idCount = totalNECount;
850    
851     /* these are the quadrilateral elements on boundary 1 (x3=0): */
852     numElementsInternal = numElementsLocal-nodesExternal[0]-nodesExternal[1];
853     if (!periodic[2]) {
854     /* elements on boundary 100 (x3=0): */
855    
856 bcumming 782 #pragma omp parallel for private(i0,i1,k)
857 bcumming 751 for (i1=0;i1<NE1;i1++) {
858 bcumming 782 for (i0=0; i0<numElementsLocal; i0++) {
859     k=i0+numElementsLocal*i1+faceNECount;
860     kk=i0 + i1*numElementsLocal;
861     facePerm = face2;
862 bcumming 751
863     out->FaceElements->Id[k]=idCount++;
864     out->FaceElements->Tag[k]=100;
865 bcumming 782 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
866     out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
867 bcumming 751
868 bcumming 782 for( j=0; j<NUMNODES; j++ )
869     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
870 bcumming 751 }
871     }
872 bcumming 782 if( boundaryLeft ){
873     for( i1=0; i1<NE1; i1++ )
874     out->FaceElements->Dom[faceNECount+i1*numElementsLocal]=ELEMENT_BOUNDARY;
875     if( periodicLocal[0] )
876     for( i1=0; i1<NE1; i1++ )
877     out->FaceElements->Dom[faceNECount+i1*numElementsLocal+1]=ELEMENT_BOUNDARY;
878     }
879     if( boundaryRight )
880     for( i1=0; i1<NE1; i1++ )
881     out->FaceElements->Dom[faceNECount+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
882     totalNECount+=NE1*numElementsLocal;
883     faceNECount+=NE1*numElementsLocal;
884 bcumming 751
885     /* ** elements on boundary 200 (x3=1) */
886    
887 bcumming 782 #pragma omp parallel for private(i0,i1,k)
888 bcumming 751 for (i1=0;i1<NE1;i1++) {
889 bcumming 782 for (i0=0;i0<numElementsLocal;i0++) {
890     k=i0+numElementsLocal*i1+faceNECount;
891     kk=i0+i1*numElementsLocal+numElementsLocal*NE1*(NE2-1);
892     facePerm = face3;
893 bcumming 751
894     out->FaceElements->Id[k]=idCount++;
895     out->FaceElements->Tag[k]=200;
896 bcumming 782 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
897     out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
898 bcumming 751
899 bcumming 782 for( j=0; j<NUMNODES; j++ )
900     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
901 bcumming 751 }
902     }
903 bcumming 782 if( boundaryLeft ){
904     for( i1=0; i1<NE1; i1++ )
905     out->FaceElements->Dom[faceNECount+i1*numElementsLocal]=ELEMENT_BOUNDARY;
906     if( periodicLocal[0] )
907     for( i1=0; i1<NE1; i1++ )
908     out->FaceElements->Dom[faceNECount+i1*numElementsLocal+1]=ELEMENT_BOUNDARY;
909     }
910     if( boundaryRight )
911     for( i1=0; i1<NE1; i1++ )
912     out->FaceElements->Dom[faceNECount+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
913     totalNECount+=NE1*numElementsLocal;
914     faceNECount+=NE1*numElementsLocal;
915 bcumming 751 }
916     if (!periodic[0] && !domInternal) {
917     /* ** elements on boundary 001 (x1=0): */
918     if( domLeft ){
919 bcumming 782 #pragma omp parallel for private(i1,i2,k)
920 bcumming 751 for (i2=0;i2<NE2;i2++) {
921     for (i1=0;i1<NE1;i1++) {
922     k=i1+NE1*i2+faceNECount;
923 bcumming 782 kk=i1*numElementsLocal + i2*numElementsLocal*NE1;
924     facePerm = face0;
925 bcumming 751
926     out->FaceElements->Id[k]=idCount++;
927     out->FaceElements->Tag[k]=1;
928 bcumming 782 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
929     out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
930 bcumming 751
931 bcumming 782 for( j=0; j<NUMNODES; j++ )
932     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
933 bcumming 751 }
934     }
935     totalNECount+=NE1*NE2;
936     faceNECount+=NE1*NE2;
937     }
938     /* ** elements on boundary 002 (x1=1): */
939     if( domRight ) {
940 bcumming 782 #pragma omp parallel for private(i1,i2,k)
941 bcumming 751 for (i2=0;i2<NE2;i2++) {
942     for (i1=0;i1<NE1;i1++) {
943     k=i1+NE1*i2+faceNECount;
944 bcumming 782 kk=(i1+1)*numElementsLocal + i2*numElementsLocal*NE1 - 1;
945     facePerm = face1;
946 bcumming 751
947     out->FaceElements->Id[k]=idCount++;
948     out->FaceElements->Tag[k]=2;
949 bcumming 782 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
950     out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
951 bcumming 751
952 bcumming 782 for( j=0; j<NUMNODES; j++ )
953     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
954 bcumming 751 }
955     }
956     totalNECount+=NE1*NE2;
957     faceNECount+=NE1*NE2;
958     }
959     }
960     if (!periodic[1]) {
961     /* ** elements on boundary 010 (x2=0): */
962    
963 bcumming 782 #pragma omp parallel for private(i0,i2,k)
964 bcumming 751 for (i2=0;i2<NE2;i2++) {
965 bcumming 782 for (i0=0;i0<numElementsLocal;i0++) {
966     k=i0+numElementsLocal*i2+faceNECount;
967     kk=i0+numElementsLocal*NE1*i2;
968     facePerm = face4;
969 bcumming 751
970     out->FaceElements->Id[k]=idCount++;
971     out->FaceElements->Tag[k]=10;
972 bcumming 782 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
973     out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
974 bcumming 751
975 bcumming 782 for( j=0; j<NUMNODES; j++ )
976     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
977 bcumming 751 }
978     }
979 bcumming 782 if( boundaryLeft ){
980     for( i2=0; i2<NE2; i2++ )
981     out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;
982     if( periodicLocal[0] )
983     for( i2=0; i2<NE2; i2++ )
984     out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;
985     }
986     if( boundaryRight )
987     for( i2=0; i2<NE2; i2++ )
988     out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
989     totalNECount+=numElementsLocal*NE2;
990     faceNECount+=numElementsLocal*NE2;
991 bcumming 751
992     /* ** elements on boundary 020 (x2=1): */
993    
994 bcumming 782 #pragma omp parallel for private(i0,i2,k)
995 bcumming 751 for (i2=0;i2<NE2;i2++) {
996 bcumming 782 for (i0=0;i0<numElementsLocal;i0++) {
997     k=i0+numElementsLocal*i2+faceNECount;
998     kk=i0+numElementsLocal*NE1*(i2+1)-numElementsLocal;
999     facePerm = face5;
1000 bcumming 751
1001     out->FaceElements->Tag[k]=20;
1002     out->FaceElements->Id[k]=idCount++;
1003 bcumming 782 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
1004 bcumming 751 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
1005    
1006 bcumming 782 for( j=0; j<NUMNODES; j++ )
1007     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
1008 bcumming 751 }
1009     }
1010 bcumming 782 if( boundaryLeft ){
1011     for( i2=0; i2<NE2; i2++ )
1012     out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;
1013     if( periodicLocal[0] )
1014     for( i2=0; i2<NE2; i2++ )
1015     out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;
1016     }
1017     if( boundaryRight )
1018     for( i2=0; i2<NE2; i2++ )
1019     out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
1020     totalNECount+=numElementsLocal*NE2;
1021     faceNECount+=numElementsLocal*NE2;
1022 bcumming 751 }
1023     out->FaceElements->elementDistribution->numInternal = faceNECount;
1024    
1025     out->FaceElements->minColor=0;
1026 bcumming 782 out->FaceElements->maxColor=23;
1027     out->FaceElements->numElements=faceNECount;
1028    
1029     Finley_ElementFile_setDomainFlags( out->FaceElements );
1030 bcumming 751
1031 bcumming 782 /* setup distribution info for other elements */
1032     Finley_ElementFile_setDomainFlags( out->ContactElements );
1033     Finley_ElementFile_setDomainFlags( out->Points );
1034 bcumming 751
1035 bcumming 782 /* reorder the degrees of freedom */
1036     Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
1037 bcumming 751
1038     /* condense the nodes: */
1039 bcumming 782 Finley_Mesh_resolveNodeIds(out);
1040     if( !Finley_MPI_noError(mpi_info) )
1041     {
1042 bcumming 751 Paso_MPIInfo_dealloc( mpi_info );
1043     Finley_Mesh_dealloc(out);
1044     return NULL;
1045 bcumming 782 }
1046 bcumming 751
1047     /* prepare mesh for further calculatuions:*/
1048 bcumming 782 Finley_Mesh_prepare(out);
1049     if( !Finley_MPI_noError(mpi_info) )
1050     {
1051     Paso_MPIInfo_dealloc( mpi_info );
1052     Finley_Mesh_dealloc(out);
1053     return NULL;
1054     }
1055 bcumming 751
1056 bcumming 782 /* free up memory */
1057     Paso_MPIInfo_dealloc( mpi_info );
1058 bcumming 751
1059     #ifdef Finley_TRACE
1060     printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
1061     #endif
1062    
1063 bcumming 782 return out;
1064 bcumming 751 }
1065     #endif
1066    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26