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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 782 - (hide annotations)
Tue Jul 18 00:47:47 2006 UTC (13 years, 4 months ago) by bcumming
File MIME type: text/plain
File size: 47192 byte(s)
Large number of changes to Finley for meshing in MPI.

- optimisation and neatening up of rectcanglular mesh generation code
- first and second order 1D, 2D and 3D rectangular meshes are now
  available in finley and escript using MPI.
- reduced meshes now generated in MPI, and interpolation to and from 
  reduced data types now supported.  

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26