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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 787 - (hide annotations)
Wed Jul 26 01:46:45 2006 UTC (13 years, 4 months ago) by bcumming
File MIME type: text/plain
File size: 28903 byte(s)
MPI update
Each element (normal elements, faceElements, ContactElements and points)
are now assigned a unique global id to streamline per-element
calculations and file IO of element data.



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] mesh with second order elements (Rec8) in the rectangle */
18     /* [0,Length[0]] x [0,Length[1]]. order is the desired accuracy of the integration scheme. */
19    
20     /**************************************************************/
21    
22     /* Author: gross@access.edu.au */
23     /* Version: $Id$ */
24    
25     /**************************************************************/
26    
27     #include "RectangularMesh.h"
28    
29 bcumming 782 #ifdef PASO_MPI
30     /* get the number of nodes/elements for domain with rank=rank, of size processors
31     where n is the total number of nodes/elements in the global domain */
32     static index_t domain_MODdim( index_t rank, index_t size, index_t n )
33     {
34     rank = size-rank-1;
35    
36     if( rank < n%size )
37     return (index_t)floor(n/size)+1;
38     return (index_t)floor(n/size);
39     }
40    
41    
42     /* Determines the number of nodes/elements etc along an axis which is numElementsGlobal long for domain rank */
43     /* A bit messy, but it only has to be done once... */
44     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 )
45     {
46     index_t i0;
47     dim_t numNodesGlobal = numElementsGlobal+1;
48    
49     (*numNodesLocal) = 2*domain_MODdim( rank, size, numNodesGlobal ) - 1;
50     if( rank<size-1 ) // add on node for right hand boundary
51     (*numNodesLocal) += 1;
52    
53     numElementsLocal[0] = domain_MODdim( rank, size, numNodesGlobal )+1;
54     periodicLocal[0] = periodicLocal[1] = FALSE;
55     nodesExternal[0] = nodesExternal[1] = 1;
56     if( periodic )
57     {
58     if( size==1 )
59     {
60     numElementsLocal[0] = numElementsGlobal;
61     nodesExternal[0] = nodesExternal[1] = 0;
62     periodicLocal[0] = periodicLocal[1] = TRUE;
63     }
64     else
65     {
66     if( rank==0 )
67     {
68     periodicLocal[0] = TRUE;
69     numNodesLocal[0]++;
70     }
71     if( rank==(size-1) )
72     {
73     periodicLocal[1] = TRUE;
74     numNodesLocal[0]--;
75     numElementsLocal[0]--;
76     }
77     }
78     }
79     else if( !periodic )
80     {
81     if( rank==0 ){
82     nodesExternal[0]--;
83     numElementsLocal[0]--;
84     }
85     if( rank==(size-1) )
86     {
87     nodesExternal[1]--;
88     numElementsLocal[0]--;
89     }
90     }
91     nodesExternal[0]*=2;
92     numNodesExternal[0] = nodesExternal[0]+nodesExternal[1];
93    
94     numElementsInternal[0] = numElementsLocal[0];
95     if( (rank==0) && (rank==size-1) );
96    
97     else if( !periodic && ( (rank==0) ^ (rank==size-1) ) )
98     numElementsInternal[0] -= 1;
99     else
100     numElementsInternal[0] -= 2;
101    
102     firstNode[0] = 0;
103     for( i0=0; i0<rank; i0++ )
104     firstNode[0] += domain_MODdim( i0, size, numNodesGlobal );
105     firstNode[0] *= 2;
106    
107     numDOFLocal[0] = numNodesLocal[0];
108     if( periodicLocal[0] )
109     numDOFLocal[0]--;
110    
111     DOFExternal[0] = nodesExternal[0];
112     DOFExternal[1] = nodesExternal[1];
113    
114     if( size>1 )
115     {
116     DOFBoundary[0] = periodicLocal[0]*1 + (rank>0)*1;
117     DOFBoundary[1] = periodicLocal[1]*2 + (rank<(size-1))*2;
118     }
119     else
120     {
121     DOFBoundary[0] = DOFBoundary[1] = 0;
122     }
123     }
124    
125    
126     #endif
127 jgs 82 /**************************************************************/
128    
129 bcumming 782 #ifdef PASO_MPI
130     Finley_Mesh* Finley_RectangularMesh_Rec8_singleCPU(int* numElements,double* Length,int* periodic,int order,int useElementsOnFace, Paso_MPIInfo *mpi_info)
131     #else
132     Finley_Mesh* Finley_RectangularMesh_Rec8(int* numElements,double* Length,int* periodic,int order,int useElementsOnFace)
133     #endif
134     {
135 jgs 153 dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES,M0,M1;
136 jgs 123 index_t k,node0;
137 bcumming 730 Finley_Mesh* out=NULL;
138 jgs 82 char name[50];
139     double time0=Finley_timer();
140     NE0=MAX(1,numElements[0]);
141     NE1=MAX(1,numElements[1]);
142     N0=2*NE0+1;
143     N1=2*NE1+1;
144    
145 bcumming 782 #ifndef PASO_MPI
146 jgs 153 if (N0<=N1) {
147     M0=1;
148     M1=N0;
149     } else {
150     M0=N1;
151     M1=1;
152     }
153 bcumming 782 #else
154     M0=1;
155     M1=N0;
156     #endif
157 jgs 153
158 jgs 82 NFaceElements=0;
159     if (!periodic[0]) {
160     NDOF0=N0;
161     NFaceElements+=2*NE1;
162     } else {
163     NDOF0=N0-1;
164     }
165     if (!periodic[1]) {
166     NDOF1=N1;
167     NFaceElements+=2*NE0;
168     } else {
169     NDOF1=N1-1;
170     }
171    
172     /* allocate mesh: */
173    
174     sprintf(name,"Rectangular %d x %d mesh",N0,N1);
175 bcumming 782 #ifdef PASO_MPI
176     out=Finley_Mesh_alloc(name,2,order,mpi_info);
177     #else
178 jgs 82 out=Finley_Mesh_alloc(name,2,order);
179 bcumming 782 #endif
180 jgs 150 if (! Finley_noError()) return NULL;
181 jgs 82
182 bcumming 782 #ifdef PASO_MPI
183     out->Elements=Finley_ElementFile_alloc(Rec8,out->order,mpi_info);
184     if (useElementsOnFace) {
185     out->FaceElements=Finley_ElementFile_alloc(Rec8Face,out->order,mpi_info);
186     out->ContactElements=Finley_ElementFile_alloc(Rec8Face_Contact,out->order,mpi_info);
187     } else {
188     out->FaceElements=Finley_ElementFile_alloc(Line3,out->order,mpi_info);
189     out->ContactElements=Finley_ElementFile_alloc(Line3_Contact,out->order,mpi_info);
190     }
191     out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
192     #else
193 jgs 82 out->Elements=Finley_ElementFile_alloc(Rec8,out->order);
194     if (useElementsOnFace) {
195     out->FaceElements=Finley_ElementFile_alloc(Rec8Face,out->order);
196     out->ContactElements=Finley_ElementFile_alloc(Rec8Face_Contact,out->order);
197     } else {
198     out->FaceElements=Finley_ElementFile_alloc(Line3,out->order);
199     out->ContactElements=Finley_ElementFile_alloc(Line3_Contact,out->order);
200     }
201     out->Points=Finley_ElementFile_alloc(Point1,out->order);
202 bcumming 730 #endif
203 bcumming 782
204 jgs 150 if (! Finley_noError()) {
205 jgs 82 Finley_Mesh_dealloc(out);
206     return NULL;
207     }
208    
209     /* allocate tables: */
210     Finley_NodeFile_allocTable(out->Nodes,N0*N1);
211 bcumming 787 Finley_ElementFile_allocTable(out->Elements,NE0*NE1);
212     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
213 bcumming 782 #ifdef PASO_MPI
214     Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0*NDOF1, 0, 0 );
215 bcumming 787 Finley_ElementDistribution_allocTable( out->Elements->elementDistribution, NE0*NE1, NE0*NE1);
216     Finley_ElementDistribution_allocTable( out->FaceElements->elementDistribution, NFaceElements, NFaceElements );
217 bcumming 730 #endif
218 jgs 150 if (! Finley_noError()) {
219 jgs 82 Finley_Mesh_dealloc(out);
220     return NULL;
221     }
222 bcumming 782
223 jgs 82 /* set nodes: */
224 jgs 153 #pragma omp parallel for private(i0,i1,k)
225 jgs 82 for (i1=0;i1<N1;i1++) {
226     for (i0=0;i0<N0;i0++) {
227 jgs 153 k=M0*i0+M1*i1;
228 jgs 82 out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(i0)/DBLE(N0-1)*Length[0];
229     out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
230 jgs 153 out->Nodes->Id[k]=i0+N0*i1;
231 jgs 82 out->Nodes->Tag[k]=0;
232 jgs 153 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1);
233 bcumming 782 #ifdef PASO_MPI
234     out->Nodes->Dom[k]=NODE_INTERNAL;
235     #endif
236 jgs 82 }
237     }
238     /* tags for the faces: */
239     if (!periodic[1]) {
240     for (i0=0;i0<N0;i0++) {
241 jgs 153 out->Nodes->Tag[M0*i0+M1*0]+=10;
242     out->Nodes->Tag[M0*i0+M1*(N1-1)]+=20;
243 jgs 82 }
244     }
245     if (!periodic[0]) {
246     for (i1=0;i1<N1;i1++) {
247 jgs 153 out->Nodes->Tag[M0*0+M1*i1]+=1;
248     out->Nodes->Tag[M0*(N0-1)+M1*i1]+=2;
249 jgs 82 }
250     }
251 jgs 153
252 jgs 82 /* set the elements: */
253    
254     #pragma omp parallel for private(i0,i1,k,node0)
255     for (i1=0;i1<NE1;i1++) {
256     for (i0=0;i0<NE0;i0++) {
257     k=i0+NE0*i1;
258     node0=2*i0+2*i1*N0;
259    
260     out->Elements->Id[k]=k;
261     out->Elements->Tag[k]=0;
262     out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);
263 bcumming 782 #ifdef PASO_MPI
264     out->Elements->Dom[k]=ELEMENT_INTERNAL;
265     #endif
266 jgs 82
267     out->Elements->Nodes[INDEX2(0,k,8)]=node0;
268     out->Elements->Nodes[INDEX2(1,k,8)]=node0+2;
269     out->Elements->Nodes[INDEX2(2,k,8)]=node0+2*N0+2;
270     out->Elements->Nodes[INDEX2(3,k,8)]=node0+2*N0;
271     out->Elements->Nodes[INDEX2(4,k,8)]=node0+1;
272     out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0+2;
273     out->Elements->Nodes[INDEX2(6,k,8)]=node0+2*N0+1;
274     out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0;
275    
276     }
277     }
278 jgs 123 out->Elements->minColor=0;
279     out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
280 jgs 82
281     /* face elements: */
282     if (useElementsOnFace) {
283     NUMNODES=8;
284     } else {
285     NUMNODES=3;
286     }
287    
288     totalNECount=NE0*NE1;
289     faceNECount=0;
290    
291     if (!periodic[1]) {
292     /* ** elements on boundary 010 (x2=0): */
293    
294     #pragma omp parallel for private(i0,k,node0)
295     for (i0=0;i0<NE0;i0++) {
296     k=i0+faceNECount;
297     node0=2*i0;
298    
299     out->FaceElements->Id[k]=i0+totalNECount;
300     out->FaceElements->Tag[k]=10;
301     out->FaceElements->Color[k]=i0%2;
302 bcumming 782 #ifdef PASO_MPI
303     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
304     #endif
305 jgs 82
306     if (useElementsOnFace) {
307     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
308     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;
309     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+2;
310     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0;
311     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
312     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+2;
313     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0+1;
314     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;
315     } else {
316     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
317     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;
318     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
319     }
320     }
321     totalNECount+=NE0;
322     faceNECount+=NE0;
323    
324     /* ** elements on boundary 020 (x2=1): */
325    
326     #pragma omp parallel for private(i0,k,node0)
327     for (i0=0;i0<NE0;i0++) {
328     k=i0+faceNECount;
329     node0=2*i0+2*(NE1-1)*N0;
330    
331     out->FaceElements->Id[k]=i0+totalNECount;
332     out->FaceElements->Tag[k]=20;
333     out->FaceElements->Color[k]=i0%2+2;
334 bcumming 782 #ifdef PASO_MPI
335     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
336     #endif
337 jgs 82 if (useElementsOnFace) {
338     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0+2;
339     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;
340     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;
341     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2;
342     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0+1;
343     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;
344     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+1;
345     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+2;
346     } else {
347     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0+2;
348     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;
349     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+1;
350     }
351     }
352     totalNECount+=NE0;
353     faceNECount+=NE0;
354     }
355     if (!periodic[0]) {
356     /* ** elements on boundary 001 (x1=0): */
357    
358     #pragma omp parallel for private(i1,k,node0)
359     for (i1=0;i1<NE1;i1++) {
360     k=i1+faceNECount;
361     node0=2*i1*N0;
362    
363     out->FaceElements->Id[k]=i1+totalNECount;
364     out->FaceElements->Tag[k]=1;
365     out->FaceElements->Color[k]=(i1%2)+4;
366 bcumming 782 #ifdef PASO_MPI
367     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
368     #endif
369 jgs 82
370     if (useElementsOnFace) {
371     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;
372     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
373     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2;
374     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0+2;
375     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;
376     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
377     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+2;
378     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0+1;
379     } else {
380     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;
381     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
382     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0;
383     }
384    
385     }
386     totalNECount+=NE1;
387     faceNECount+=NE1;
388    
389     /* ** elements on boundary 002 (x1=1): */
390    
391     #pragma omp parallel for private(i1,k,node0)
392     for (i1=0;i1<NE1;i1++) {
393     k=i1+faceNECount;
394     node0=2*(NE0-1)+2*i1*N0;
395    
396     out->FaceElements->Id[k]=i1+totalNECount;
397     out->FaceElements->Tag[k]=2;
398     out->FaceElements->Color[k]=(i1%2)+6;
399 bcumming 782 #ifdef PASO_MPI
400     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
401     #endif
402 jgs 82
403     if (useElementsOnFace) {
404     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;
405     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;
406     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0;
407     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
408     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0+2;
409     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0+1;
410     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0;
411     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
412     } else {
413     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;
414     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;
415     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+2;
416     }
417     }
418     totalNECount+=NE1;
419     faceNECount+=NE1;
420    
421     }
422 jgs 123 out->FaceElements->minColor=0;
423     out->FaceElements->maxColor=7;
424 jgs 82
425 bcumming 782 #ifdef PASO_MPI
426     Finley_ElementFile_setDomainFlags( out->Elements );
427     Finley_ElementFile_setDomainFlags( out->FaceElements );
428     Finley_ElementFile_setDomainFlags( out->ContactElements );
429     Finley_ElementFile_setDomainFlags( out->Points );
430    
431     /* reorder the degrees of freedom */
432     Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
433     #endif
434    
435 jgs 82 /* condense the nodes: */
436     Finley_Mesh_resolveNodeIds(out);
437 bcumming 782
438 bcumming 787 if (! Finley_noError()) {
439     Finley_Mesh_dealloc(out);
440     return NULL;
441     }
442 jgs 82 /* prepare mesh for further calculatuions:*/
443     Finley_Mesh_prepare(out) ;
444    
445 bcumming 782
446 jgs 149 #ifdef Finley_TRACE
447 jgs 82 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
448 jgs 149 #endif
449 jgs 82
450 jgs 150 if (! Finley_noError()) {
451 jgs 82 Finley_Mesh_dealloc(out);
452     return NULL;
453     }
454     return out;
455     }
456 bcumming 782
457     #ifdef PASO_MPI
458     Finley_Mesh* Finley_RectangularMesh_Rec8(int* numElements,double* Length,int* periodic,int order,int useElementsOnFace) {
459     dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal, DOFBoundary[2];
460     dim_t N0t, NDOF0t;
461     index_t *numForward=NULL, *numBackward=NULL;
462     index_t kk, k,i,j,node0,firstNodeConstruct,firstNode=0, DOFcount=0, forwardDOF=NULL, backwardDOF=NULL, targetDomain=-1,faceNEBoundary;
463     Finley_Mesh* out=NULL;
464     char name[50];
465     index_t *indexForward=NULL, *indexBackward=NULL;
466     double time0=Finley_timer();
467     bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft=FALSE, boundaryRight=FALSE;
468     Paso_MPIInfo *mpi_info = NULL;
469 bcumming 787 index_t faceStart, numFaceLeft;
470 bcumming 782
471     /* these are used in constructing the face elements, and give the permutation
472     of the node order of a normal element for each face */
473     index_t *facePerm;
474     index_t face0[]={0, 1, 2, 3, 4, 5, 6, 7};
475     index_t face1[]={2, 3, 0, 1, 6, 7, 4, 5};
476     index_t face2[]={3, 0, 1, 2, 7, 4, 5, 6};
477     index_t face3[]={1, 2, 3, 0, 5, 6, 7, 4};
478     index_t face0nc[]={0, 1, 4};
479     index_t face1nc[]={2, 3, 6};
480     index_t face2nc[]={3, 0, 7};
481     index_t face3nc[]={1, 2, 5};
482    
483     NE0=MAX(1,numElements[0]);
484     NE1=MAX(1,numElements[1]);
485     N0=2*NE0+1;
486     N1=2*NE1+1;
487    
488     /* get MPI information */
489     mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
490     if (! Finley_noError()) {
491     Finley_Mesh_dealloc(out);
492     return NULL;
493     }
494    
495     // use the serial method for the single CPU case
496     if( mpi_info->size==1 ){
497     out = Finley_RectangularMesh_Rec8_singleCPU( numElements,Length,periodic,order,useElementsOnFace, mpi_info);
498     return out;
499     }
500    
501     if( mpi_info->rank==0 )
502     domLeft = TRUE;
503     if( mpi_info->rank==mpi_info->size-1 )
504     domRight = TRUE;
505     if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
506     domInternal = TRUE;
507    
508     /* dimensions of the local subdomain */
509     domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal, DOFBoundary );
510    
511     NFaceElements=0;
512     if (!periodic[0]) {
513     NDOF0=N0;
514     NFaceElements+=(domLeft + domRight)*NE1;
515     } else {
516     NDOF0=N0-1;
517     }
518     if (!periodic[1]) {
519     NDOF1=N1;
520     NFaceElements+=2*numElementsLocal;
521     } else {
522     NDOF1=N1-1;
523     }
524    
525     boundaryLeft = !domLeft || periodicLocal[0];
526     boundaryRight = !domRight || periodicLocal[1];
527     N0t = numNodesLocal + boundaryRight + boundaryLeft*2;
528     NDOF0t = numDOFLocal + boundaryRight + boundaryLeft*2;
529     firstNodeConstruct = firstNode - 2*boundaryLeft;
530     firstNodeConstruct = firstNodeConstruct<0 ? N0+firstNodeConstruct-1 : firstNodeConstruct;
531    
532     /* allocate mesh: */
533    
534     sprintf(name,"Rectangular %d x %d mesh",N0,N1);
535     out=Finley_Mesh_alloc(name,2,order,mpi_info);
536    
537     if (! Finley_noError()) return NULL;
538    
539     out->Elements=Finley_ElementFile_alloc(Rec8,out->order,mpi_info);
540     if (useElementsOnFace) {
541     out->FaceElements=Finley_ElementFile_alloc(Rec8Face,out->order,mpi_info);
542     out->ContactElements=Finley_ElementFile_alloc(Rec8Face_Contact,out->order,mpi_info);
543     } else {
544     out->FaceElements=Finley_ElementFile_alloc(Line3,out->order,mpi_info);
545     out->ContactElements=Finley_ElementFile_alloc(Line3_Contact,out->order,mpi_info);
546     }
547     out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
548    
549     if (! Finley_noError()) {
550     Finley_Mesh_dealloc(out);
551     return NULL;
552     }
553    
554     /* allocate tables: */
555     Finley_NodeFile_allocTable(out->Nodes,N0t*N1);
556     Finley_ElementFile_allocTable(out->Elements,numElementsLocal*NE1);
557     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
558 bcumming 787
559     Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0t*NDOF1, NDOF1*3, 0 );
560     Finley_ElementDistribution_allocTable( out->Elements->elementDistribution, numElementsLocal*NE1, NE1*(numElementsLocal-boundaryRight*(!periodic[1])) );
561     Finley_ElementDistribution_allocTable( out->FaceElements->elementDistribution, NFaceElements, NFaceElements-2*boundaryRight*(!periodic[1]) );
562 bcumming 782 if (! Finley_noError()) {
563     Finley_Mesh_dealloc(out);
564     return NULL;
565     }
566    
567     /* set nodes: */
568     #pragma omp parallel for private(i0,i1,k)
569 bcumming 787 for (i1=0;i1<N1;i1++) {
570 bcumming 782 for (i0=0;i0<N0t;i0++, k++) {
571 bcumming 787 k=i0+i1*N0t;
572 bcumming 782 out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE((firstNodeConstruct + i0) % N0)/DBLE(N0-1)*Length[0];
573     out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
574     out->Nodes->Id[k]=i0 + i1*N0t;
575     out->Nodes->Tag[k]=0;
576     out->Nodes->Dom[k]=NODE_INTERNAL;
577     out->Nodes->degreeOfFreedom[k]=i0 + (i1%NDOF1)*N0t;
578     }
579     /* take care to make the DOF reflect the internal/external position of the DOF */
580     if( boundaryLeft ){
581     out->Nodes->Dom[i1*N0t] = NODE_EXTERNAL;
582     out->Nodes->Dom[i1*N0t+1] = NODE_EXTERNAL;
583     out->Nodes->Dom[i1*N0t+2] = NODE_BOUNDARY;
584     if( periodicLocal[0] ){
585     out->Nodes->degreeOfFreedom[i1*N0t+3] = out->Nodes->degreeOfFreedom[i1*N0t+2];
586     out->Nodes->Dom[i1*N0t+3] = NODE_BOUNDARY;
587     }
588     }
589     if( boundaryRight ){
590     out->Nodes->Dom[(i1+1)*N0t-1] = NODE_EXTERNAL;
591     out->Nodes->Dom[(i1+1)*N0t-2] = NODE_BOUNDARY;
592     out->Nodes->Dom[(i1+1)*N0t-3] = NODE_BOUNDARY;
593     }
594     }
595     /* tags for the faces: */
596     if (!periodic[1]) {
597     for (i0=0;i0<N0t;i0++) {
598     out->Nodes->Tag[i0]+=10;
599     out->Nodes->Tag[i0+N0t*(N1-1)]+=20;
600     }
601     }
602     if (domLeft && !periodic[0]) {
603     for (i1=0;i1<N1;i1++) {
604     out->Nodes->Tag[i1*N0t]+=1;
605     }
606     }
607     if (domRight && !periodic[0]){
608     for (i1=0;i1<N1;i1++) {
609     out->Nodes->Tag[(i1+1)*N0t-1]+=2;
610     }
611     }
612    
613     /* setup the receive information */
614     indexBackward = TMPMEMALLOC( NDOF1*2, index_t );
615     indexForward = TMPMEMALLOC( NDOF1*2, index_t );
616     if( !(mpi_info->size==2 && periodicLocal[0])){
617     /* Backward (receive) DOF */
618     if( boundaryLeft ){
619     targetDomain = mpi_info->rank-1<0 ? mpi_info->size-1 : mpi_info->rank-1;
620    
621     for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){
622     indexBackward[i1] = out->Nodes->degreeOfFreedom[i*N0t];
623     indexBackward[i1+1] = out->Nodes->degreeOfFreedom[i*N0t+1];
624     indexForward[i] = out->Nodes->degreeOfFreedom[i*N0t+2];
625     }
626    
627     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );
628     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexBackward );
629     }
630    
631     if( boundaryRight ){
632     targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;
633    
634     for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){
635     indexBackward[i] = out->Nodes->degreeOfFreedom[(i+1)*N0t-1];
636     indexForward[i1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-3];
637     indexForward[i1+1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-2];
638     }
639    
640     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexForward );
641     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexBackward );
642     }
643     }
644     else{
645     targetDomain = 1;
646    
647     for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){
648     indexBackward[i] = out->Nodes->degreeOfFreedom[(i+1)*N0t-1];
649     indexForward[i1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-3];
650     indexForward[i1+1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-2];
651     }
652    
653     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexForward );
654     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexBackward );
655    
656     for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){
657     indexBackward[i1] = out->Nodes->degreeOfFreedom[i*N0t];
658     indexBackward[i1+1] = out->Nodes->degreeOfFreedom[i*N0t+1];
659     indexForward[i] = out->Nodes->degreeOfFreedom[i*N0t+2];
660     }
661    
662     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );
663     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexBackward );
664     }
665     TMPMEMFREE( indexBackward );
666     TMPMEMFREE( indexForward );
667    
668     /* set the elements: */
669     #pragma omp parallel for private(i0,i1,k,node0)
670     for (i1=0;i1<NE1;i1++) {
671     for (i0=0;i0<numElementsLocal;i0++) {
672     k=i0+numElementsLocal*i1;
673     node0 = (periodicLocal[0] && !i0) ? 2*i1*N0t : 2*i0+2*i1*N0t+periodicLocal[0] ;
674    
675 bcumming 787 out->Elements->Id[k]=((firstNodeConstruct/2+i0)%NE0) * NE1 + i1;
676 bcumming 782 out->Elements->Tag[k]=0;
677     out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);
678     out->Elements->Dom[k] = ELEMENT_INTERNAL;
679    
680     out->Elements->Nodes[INDEX2(0,k,8)]=node0;
681     out->Elements->Nodes[INDEX2(1,k,8)]=node0+2;
682     out->Elements->Nodes[INDEX2(2,k,8)]=node0+2*N0t+2;
683     out->Elements->Nodes[INDEX2(3,k,8)]=node0+2*N0t;
684     out->Elements->Nodes[INDEX2(4,k,8)]=node0+1;
685     out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0t+2;
686     out->Elements->Nodes[INDEX2(6,k,8)]=node0+2*N0t+1;
687     out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0t;
688     }
689     if( boundaryLeft )
690     out->Elements->Dom[i1*numElementsLocal] = ELEMENT_BOUNDARY;
691     if( boundaryRight )
692     out->Elements->Dom[(i1+1)*numElementsLocal-1] = ELEMENT_BOUNDARY;
693     }
694     out->Elements->minColor=0;
695     out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
696    
697     /* set the distribution information for the elements */
698     Finley_ElementFile_setDomainFlags( out->Elements );
699    
700     /* face elements: */
701     if (useElementsOnFace) {
702     NUMNODES=8;
703     } else {
704     NUMNODES=3;
705     }
706    
707 bcumming 787 Finley_ElementFile_setDomainFlags( out->FaceElements );
708 bcumming 782 totalNECount=numElementsLocal*NE1;
709     faceNECount=0;
710     faceNEBoundary = totalNECount + numElementsInternal*(!periodic[1])*2 + NE1*(!periodic[0])*(domLeft+domRight);
711    
712 bcumming 787 numFaceLeft = domLeft*(!periodic[0])*NE1;
713     faceStart = out->FaceElements->elementDistribution->vtxdist[mpi_info->rank];
714 bcumming 782 if (!periodic[1]) {
715 bcumming 787 /* ** elements on boundary 010 (x1=0): */
716 bcumming 782
717 bcumming 787 #pragma omp parallel for private(i0,k,kk,facePerm)
718 bcumming 782 for (i0=0;i0<numElementsLocal; i0++) {
719     k=i0+faceNECount;
720     kk=i0;
721     facePerm = useElementsOnFace ? face0 : face0nc;
722    
723 bcumming 787 out->FaceElements->Id[k]=faceStart + numFaceLeft + i0*2;
724 bcumming 782 out->FaceElements->Tag[k]=10;
725     out->FaceElements->Color[k]=i0%2;
726     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
727    
728     for( j=0; j<NUMNODES; j++ )
729     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
730     }
731     if( boundaryLeft ){
732     out->FaceElements->Dom[faceNECount]=ELEMENT_BOUNDARY;
733     if( periodicLocal[0] )
734     out->FaceElements->Dom[faceNECount+1]=ELEMENT_BOUNDARY;
735     }
736 bcumming 787 if( boundaryRight ){
737 bcumming 782 out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;
738 bcumming 787 out->FaceElements->Id[faceNECount+numElementsLocal-1]=out->FaceElements->elementDistribution->vtxdist[mpi_info->rank+1];
739     }
740 bcumming 782 totalNECount+=numElementsLocal;
741     faceNECount+=numElementsLocal;
742    
743 bcumming 787 /* ** elements on boundary 020 (x1=1): */
744 bcumming 782
745 bcumming 787 #pragma omp parallel for private(i0,k,kk,facePerm)
746 bcumming 782 for (i0=0;i0<numElementsLocal;i0++) {
747     k=i0+faceNECount;
748     kk=i0 + numElementsLocal*(NE1-1);
749     facePerm = useElementsOnFace ? face1 : face1nc;
750    
751 bcumming 787 out->FaceElements->Id[k]=faceStart + numFaceLeft + i0*2+1;
752 bcumming 782 out->FaceElements->Tag[k]=20;
753     out->FaceElements->Color[k]=i0%2+2;
754     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
755    
756     for( j=0; j<NUMNODES; j++ )
757     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
758     }
759     if( boundaryLeft ){
760     out->FaceElements->Dom[faceNECount]=ELEMENT_BOUNDARY;
761     if( periodicLocal[0] )
762     out->FaceElements->Dom[faceNECount+1]=ELEMENT_BOUNDARY;
763     }
764 bcumming 787 if( boundaryRight ){
765 bcumming 782 out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;
766 bcumming 787 out->FaceElements->Id[faceNECount+numElementsLocal-1]=out->FaceElements->elementDistribution->vtxdist[mpi_info->rank+1]+1;
767     }
768 bcumming 782 totalNECount+=numElementsLocal;
769     faceNECount+=numElementsLocal;
770     }
771     if (domLeft && !periodicLocal[0]) {
772 bcumming 787 /* ** elements on boundary 001 (x0=0): */
773 bcumming 782
774 bcumming 787 #pragma omp parallel for private(i1,k,kk,facePerm)
775 bcumming 782 for (i1=0;i1<NE1;i1++) {
776     k=i1+faceNECount;
777     kk=i1*numElementsLocal;
778     facePerm = useElementsOnFace ? face2 : face2nc;
779    
780 bcumming 787 out->FaceElements->Id[k]=faceStart + i1;
781 bcumming 782 out->FaceElements->Tag[k]=1;
782     out->FaceElements->Color[k]=(i1%2)+4;
783     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
784    
785     for( j=0; j<NUMNODES; j++ )
786     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
787     }
788     totalNECount+=NE1;
789     faceNECount+=NE1;
790    
791 bcumming 787 /* ** elements on boundary 002 (x0=1): */
792 bcumming 782 }
793     if (domRight && !periodicLocal[1]) {
794 bcumming 787 #pragma omp parallel for private(i1,k,kk,facePerm)
795 bcumming 782 for (i1=0;i1<NE1;i1++) {
796     k=i1+faceNECount;
797     kk=(i1+1)*numElementsLocal-1;
798     facePerm = useElementsOnFace ? face3 : face3nc;
799    
800 bcumming 787 out->FaceElements->Id[k]=faceStart + numElementsLocal*2*(!periodic[1])+i1;
801 bcumming 782 out->FaceElements->Tag[k]=2;
802     out->FaceElements->Color[k]=(i1%2)+6;
803     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
804    
805     for( j=0; j<NUMNODES; j++ )
806     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
807     }
808     totalNECount+=NE1;
809     faceNECount+=NE1;
810     }
811     out->FaceElements->minColor=0;
812     out->FaceElements->maxColor=7;
813     out->FaceElements->numElements = faceNECount;
814    
815     /* setup distribution info for other elements */
816     Finley_ElementFile_setDomainFlags( out->ContactElements );
817     Finley_ElementFile_setDomainFlags( out->Points );
818    
819 bcumming 787 Finley_Mesh_prepareElementDistribution( out );
820    
821 bcumming 782 /* reorder the degrees of freedom */
822     Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
823    
824     /* condense the nodes: */
825     Finley_Mesh_resolveNodeIds(out);
826     if ( !Finley_MPI_noError( mpi_info )) {
827     Paso_MPIInfo_dealloc( mpi_info );
828     Finley_Mesh_dealloc(out);
829     return NULL;
830     }
831    
832     /* prepare mesh for further calculatuions:*/
833     Finley_Mesh_prepare(out) ;
834    
835     #ifdef Finley_TRACE
836     printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
837     #endif
838    
839     if ( !Finley_MPI_noError( mpi_info )) {
840     Paso_MPIInfo_dealloc( mpi_info );
841     Finley_Mesh_dealloc(out);
842     return NULL;
843     }
844    
845     return out;
846     }
847     #endif

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26