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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
  ViewVC Help
Powered by ViewVC 1.1.26