/[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 782 - (hide annotations)
Tue Jul 18 00:47:47 2006 UTC (13 years, 3 months ago) by bcumming
File MIME type: text/plain
File size: 27704 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] 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 782 #ifdef PASO_MPI
212     Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0*NDOF1, 0, 0 );
213 bcumming 730 #endif
214 jgs 82 Finley_ElementFile_allocTable(out->Elements,NE0*NE1);
215     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
216 jgs 150 if (! Finley_noError()) {
217 jgs 82 Finley_Mesh_dealloc(out);
218     return NULL;
219     }
220 bcumming 782
221 jgs 82 /* set nodes: */
222 jgs 153 #pragma omp parallel for private(i0,i1,k)
223 jgs 82 for (i1=0;i1<N1;i1++) {
224     for (i0=0;i0<N0;i0++) {
225 jgs 153 k=M0*i0+M1*i1;
226 jgs 82 out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(i0)/DBLE(N0-1)*Length[0];
227     out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
228 jgs 153 out->Nodes->Id[k]=i0+N0*i1;
229 jgs 82 out->Nodes->Tag[k]=0;
230 jgs 153 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1);
231 bcumming 782 #ifdef PASO_MPI
232     out->Nodes->Dom[k]=NODE_INTERNAL;
233     #endif
234 jgs 82 }
235     }
236     /* tags for the faces: */
237     if (!periodic[1]) {
238     for (i0=0;i0<N0;i0++) {
239 jgs 153 out->Nodes->Tag[M0*i0+M1*0]+=10;
240     out->Nodes->Tag[M0*i0+M1*(N1-1)]+=20;
241 jgs 82 }
242     }
243     if (!periodic[0]) {
244     for (i1=0;i1<N1;i1++) {
245 jgs 153 out->Nodes->Tag[M0*0+M1*i1]+=1;
246     out->Nodes->Tag[M0*(N0-1)+M1*i1]+=2;
247 jgs 82 }
248     }
249 jgs 153
250 jgs 82 /* set the elements: */
251    
252     #pragma omp parallel for private(i0,i1,k,node0)
253     for (i1=0;i1<NE1;i1++) {
254     for (i0=0;i0<NE0;i0++) {
255     k=i0+NE0*i1;
256     node0=2*i0+2*i1*N0;
257    
258     out->Elements->Id[k]=k;
259     out->Elements->Tag[k]=0;
260     out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);
261 bcumming 782 #ifdef PASO_MPI
262     out->Elements->Dom[k]=ELEMENT_INTERNAL;
263     #endif
264 jgs 82
265     out->Elements->Nodes[INDEX2(0,k,8)]=node0;
266     out->Elements->Nodes[INDEX2(1,k,8)]=node0+2;
267     out->Elements->Nodes[INDEX2(2,k,8)]=node0+2*N0+2;
268     out->Elements->Nodes[INDEX2(3,k,8)]=node0+2*N0;
269     out->Elements->Nodes[INDEX2(4,k,8)]=node0+1;
270     out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0+2;
271     out->Elements->Nodes[INDEX2(6,k,8)]=node0+2*N0+1;
272     out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0;
273    
274     }
275     }
276 jgs 123 out->Elements->minColor=0;
277     out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
278 jgs 82
279     /* face elements: */
280     if (useElementsOnFace) {
281     NUMNODES=8;
282     } else {
283     NUMNODES=3;
284     }
285    
286     totalNECount=NE0*NE1;
287     faceNECount=0;
288    
289     if (!periodic[1]) {
290     /* ** elements on boundary 010 (x2=0): */
291    
292     #pragma omp parallel for private(i0,k,node0)
293     for (i0=0;i0<NE0;i0++) {
294     k=i0+faceNECount;
295     node0=2*i0;
296    
297     out->FaceElements->Id[k]=i0+totalNECount;
298     out->FaceElements->Tag[k]=10;
299     out->FaceElements->Color[k]=i0%2;
300 bcumming 782 #ifdef PASO_MPI
301     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
302     #endif
303 jgs 82
304     if (useElementsOnFace) {
305     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
306     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;
307     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+2;
308     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0;
309     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
310     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+2;
311     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0+1;
312     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;
313     } else {
314     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
315     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;
316     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
317     }
318     }
319     totalNECount+=NE0;
320     faceNECount+=NE0;
321    
322     /* ** elements on boundary 020 (x2=1): */
323    
324     #pragma omp parallel for private(i0,k,node0)
325     for (i0=0;i0<NE0;i0++) {
326     k=i0+faceNECount;
327     node0=2*i0+2*(NE1-1)*N0;
328    
329     out->FaceElements->Id[k]=i0+totalNECount;
330     out->FaceElements->Tag[k]=20;
331     out->FaceElements->Color[k]=i0%2+2;
332 bcumming 782 #ifdef PASO_MPI
333     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
334     #endif
335 jgs 82 if (useElementsOnFace) {
336     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0+2;
337     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;
338     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;
339     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2;
340     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0+1;
341     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;
342     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+1;
343     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+2;
344     } else {
345     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0+2;
346     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;
347     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+1;
348     }
349     }
350     totalNECount+=NE0;
351     faceNECount+=NE0;
352     }
353     if (!periodic[0]) {
354     /* ** elements on boundary 001 (x1=0): */
355    
356     #pragma omp parallel for private(i1,k,node0)
357     for (i1=0;i1<NE1;i1++) {
358     k=i1+faceNECount;
359     node0=2*i1*N0;
360    
361     out->FaceElements->Id[k]=i1+totalNECount;
362     out->FaceElements->Tag[k]=1;
363     out->FaceElements->Color[k]=(i1%2)+4;
364 bcumming 782 #ifdef PASO_MPI
365     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
366     #endif
367 jgs 82
368     if (useElementsOnFace) {
369     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;
370     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
371     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2;
372     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0+2;
373     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;
374     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
375     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+2;
376     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0+1;
377     } else {
378     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;
379     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
380     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0;
381     }
382    
383     }
384     totalNECount+=NE1;
385     faceNECount+=NE1;
386    
387     /* ** elements on boundary 002 (x1=1): */
388    
389     #pragma omp parallel for private(i1,k,node0)
390     for (i1=0;i1<NE1;i1++) {
391     k=i1+faceNECount;
392     node0=2*(NE0-1)+2*i1*N0;
393    
394     out->FaceElements->Id[k]=i1+totalNECount;
395     out->FaceElements->Tag[k]=2;
396     out->FaceElements->Color[k]=(i1%2)+6;
397 bcumming 782 #ifdef PASO_MPI
398     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
399     #endif
400 jgs 82
401     if (useElementsOnFace) {
402     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;
403     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;
404     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0;
405     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
406     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0+2;
407     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0+1;
408     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0;
409     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
410     } else {
411     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;
412     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;
413     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+2;
414     }
415     }
416     totalNECount+=NE1;
417     faceNECount+=NE1;
418    
419     }
420 jgs 123 out->FaceElements->minColor=0;
421     out->FaceElements->maxColor=7;
422 jgs 82
423 bcumming 782 #ifdef PASO_MPI
424     Finley_ElementFile_setDomainFlags( out->Elements );
425     Finley_ElementFile_setDomainFlags( out->FaceElements );
426     Finley_ElementFile_setDomainFlags( out->ContactElements );
427     Finley_ElementFile_setDomainFlags( out->Points );
428    
429     /* reorder the degrees of freedom */
430     Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
431     #endif
432    
433 jgs 82 /* condense the nodes: */
434     Finley_Mesh_resolveNodeIds(out);
435 bcumming 782
436 jgs 82 /* prepare mesh for further calculatuions:*/
437     Finley_Mesh_prepare(out) ;
438    
439 bcumming 782
440 jgs 149 #ifdef Finley_TRACE
441 jgs 82 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
442 jgs 149 #endif
443 jgs 82
444 jgs 150 if (! Finley_noError()) {
445 jgs 82 Finley_Mesh_dealloc(out);
446     return NULL;
447     }
448     return out;
449     }
450 bcumming 782
451     #ifdef PASO_MPI
452     Finley_Mesh* Finley_RectangularMesh_Rec8(int* numElements,double* Length,int* periodic,int order,int useElementsOnFace) {
453     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];
454     dim_t N0t, NDOF0t;
455     index_t *numForward=NULL, *numBackward=NULL;
456     index_t kk, k,i,j,node0,firstNodeConstruct,firstNode=0, DOFcount=0, forwardDOF=NULL, backwardDOF=NULL, targetDomain=-1,faceNEBoundary;
457     Finley_Mesh* out=NULL;
458     char name[50];
459     index_t *indexForward=NULL, *indexBackward=NULL;
460     double time0=Finley_timer();
461     bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft=FALSE, boundaryRight=FALSE;
462     Paso_MPIInfo *mpi_info = NULL;
463    
464     /* these are used in constructing the face elements, and give the permutation
465     of the node order of a normal element for each face */
466     index_t *facePerm;
467     index_t face0[]={0, 1, 2, 3, 4, 5, 6, 7};
468     index_t face1[]={2, 3, 0, 1, 6, 7, 4, 5};
469     index_t face2[]={3, 0, 1, 2, 7, 4, 5, 6};
470     index_t face3[]={1, 2, 3, 0, 5, 6, 7, 4};
471     index_t face0nc[]={0, 1, 4};
472     index_t face1nc[]={2, 3, 6};
473     index_t face2nc[]={3, 0, 7};
474     index_t face3nc[]={1, 2, 5};
475    
476     NE0=MAX(1,numElements[0]);
477     NE1=MAX(1,numElements[1]);
478     N0=2*NE0+1;
479     N1=2*NE1+1;
480    
481     /* get MPI information */
482     mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
483     if (! Finley_noError()) {
484     Finley_Mesh_dealloc(out);
485     return NULL;
486     }
487    
488     // use the serial method for the single CPU case
489     if( mpi_info->size==1 ){
490     out = Finley_RectangularMesh_Rec8_singleCPU( numElements,Length,periodic,order,useElementsOnFace, mpi_info);
491     return out;
492     }
493    
494     if( mpi_info->rank==0 )
495     domLeft = TRUE;
496     if( mpi_info->rank==mpi_info->size-1 )
497     domRight = TRUE;
498     if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
499     domInternal = TRUE;
500    
501     /* dimensions of the local subdomain */
502     domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal, DOFBoundary );
503    
504     NFaceElements=0;
505     if (!periodic[0]) {
506     NDOF0=N0;
507     NFaceElements+=(domLeft + domRight)*NE1;
508     } else {
509     NDOF0=N0-1;
510     }
511     if (!periodic[1]) {
512     NDOF1=N1;
513     NFaceElements+=2*numElementsLocal;
514     } else {
515     NDOF1=N1-1;
516     }
517    
518     boundaryLeft = !domLeft || periodicLocal[0];
519     boundaryRight = !domRight || periodicLocal[1];
520     N0t = numNodesLocal + boundaryRight + boundaryLeft*2;
521     NDOF0t = numDOFLocal + boundaryRight + boundaryLeft*2;
522     firstNodeConstruct = firstNode - 2*boundaryLeft;
523     firstNodeConstruct = firstNodeConstruct<0 ? N0+firstNodeConstruct-1 : firstNodeConstruct;
524    
525     /* allocate mesh: */
526    
527     sprintf(name,"Rectangular %d x %d mesh",N0,N1);
528     out=Finley_Mesh_alloc(name,2,order,mpi_info);
529    
530     if (! Finley_noError()) return NULL;
531    
532     out->Elements=Finley_ElementFile_alloc(Rec8,out->order,mpi_info);
533     if (useElementsOnFace) {
534     out->FaceElements=Finley_ElementFile_alloc(Rec8Face,out->order,mpi_info);
535     out->ContactElements=Finley_ElementFile_alloc(Rec8Face_Contact,out->order,mpi_info);
536     } else {
537     out->FaceElements=Finley_ElementFile_alloc(Line3,out->order,mpi_info);
538     out->ContactElements=Finley_ElementFile_alloc(Line3_Contact,out->order,mpi_info);
539     }
540     out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
541    
542     if (! Finley_noError()) {
543     Finley_Mesh_dealloc(out);
544     return NULL;
545     }
546    
547     /* allocate tables: */
548     Finley_NodeFile_allocTable(out->Nodes,N0t*N1);
549     Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0t*NDOF1, NDOF1*3, 0 );
550     Finley_ElementFile_allocTable(out->Elements,numElementsLocal*NE1);
551     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
552     if (! Finley_noError()) {
553     Finley_Mesh_dealloc(out);
554     return NULL;
555     }
556    
557     /* set nodes: */
558     #pragma omp parallel for private(i0,i1,k)
559     for (k=0,i1=0;i1<N1;i1++) {
560     for (i0=0;i0<N0t;i0++, k++) {
561     out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE((firstNodeConstruct + i0) % N0)/DBLE(N0-1)*Length[0];
562     out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
563     out->Nodes->Id[k]=i0 + i1*N0t;
564     out->Nodes->Tag[k]=0;
565     out->Nodes->Dom[k]=NODE_INTERNAL;
566     out->Nodes->degreeOfFreedom[k]=i0 + (i1%NDOF1)*N0t;
567     }
568     /* take care to make the DOF reflect the internal/external position of the DOF */
569     if( boundaryLeft ){
570     out->Nodes->Dom[i1*N0t] = NODE_EXTERNAL;
571     out->Nodes->Dom[i1*N0t+1] = NODE_EXTERNAL;
572     out->Nodes->Dom[i1*N0t+2] = NODE_BOUNDARY;
573     if( periodicLocal[0] ){
574     out->Nodes->degreeOfFreedom[i1*N0t+3] = out->Nodes->degreeOfFreedom[i1*N0t+2];
575     out->Nodes->Dom[i1*N0t+3] = NODE_BOUNDARY;
576     }
577     }
578     if( boundaryRight ){
579     out->Nodes->Dom[(i1+1)*N0t-1] = NODE_EXTERNAL;
580     out->Nodes->Dom[(i1+1)*N0t-2] = NODE_BOUNDARY;
581     out->Nodes->Dom[(i1+1)*N0t-3] = NODE_BOUNDARY;
582     }
583     }
584     /* tags for the faces: */
585     if (!periodic[1]) {
586     for (i0=0;i0<N0t;i0++) {
587     out->Nodes->Tag[i0]+=10;
588     out->Nodes->Tag[i0+N0t*(N1-1)]+=20;
589     }
590     }
591     if (domLeft && !periodic[0]) {
592     for (i1=0;i1<N1;i1++) {
593     out->Nodes->Tag[i1*N0t]+=1;
594     }
595     }
596     if (domRight && !periodic[0]){
597     for (i1=0;i1<N1;i1++) {
598     out->Nodes->Tag[(i1+1)*N0t-1]+=2;
599     }
600     }
601    
602     /* setup the receive information */
603     indexBackward = TMPMEMALLOC( NDOF1*2, index_t );
604     indexForward = TMPMEMALLOC( NDOF1*2, index_t );
605     if( !(mpi_info->size==2 && periodicLocal[0])){
606     /* Backward (receive) DOF */
607     if( boundaryLeft ){
608     targetDomain = mpi_info->rank-1<0 ? mpi_info->size-1 : mpi_info->rank-1;
609    
610     for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){
611     indexBackward[i1] = out->Nodes->degreeOfFreedom[i*N0t];
612     indexBackward[i1+1] = out->Nodes->degreeOfFreedom[i*N0t+1];
613     indexForward[i] = out->Nodes->degreeOfFreedom[i*N0t+2];
614     }
615    
616     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );
617     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexBackward );
618     }
619    
620     if( boundaryRight ){
621     targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;
622    
623     for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){
624     indexBackward[i] = out->Nodes->degreeOfFreedom[(i+1)*N0t-1];
625     indexForward[i1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-3];
626     indexForward[i1+1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-2];
627     }
628    
629     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexForward );
630     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexBackward );
631     }
632     }
633     else{
634     targetDomain = 1;
635    
636     for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){
637     indexBackward[i] = out->Nodes->degreeOfFreedom[(i+1)*N0t-1];
638     indexForward[i1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-3];
639     indexForward[i1+1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-2];
640     }
641    
642     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexForward );
643     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexBackward );
644    
645     for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){
646     indexBackward[i1] = out->Nodes->degreeOfFreedom[i*N0t];
647     indexBackward[i1+1] = out->Nodes->degreeOfFreedom[i*N0t+1];
648     indexForward[i] = out->Nodes->degreeOfFreedom[i*N0t+2];
649     }
650    
651     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );
652     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexBackward );
653     }
654     TMPMEMFREE( indexBackward );
655     TMPMEMFREE( indexForward );
656    
657     /* set the elements: */
658     #pragma omp parallel for private(i0,i1,k,node0)
659     for (i1=0;i1<NE1;i1++) {
660     for (i0=0;i0<numElementsLocal;i0++) {
661     k=i0+numElementsLocal*i1;
662     node0 = (periodicLocal[0] && !i0) ? 2*i1*N0t : 2*i0+2*i1*N0t+periodicLocal[0] ;
663    
664     out->Elements->Id[k]=k;
665     out->Elements->Tag[k]=0;
666     out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);
667     out->Elements->Dom[k] = ELEMENT_INTERNAL;
668    
669     out->Elements->Nodes[INDEX2(0,k,8)]=node0;
670     out->Elements->Nodes[INDEX2(1,k,8)]=node0+2;
671     out->Elements->Nodes[INDEX2(2,k,8)]=node0+2*N0t+2;
672     out->Elements->Nodes[INDEX2(3,k,8)]=node0+2*N0t;
673     out->Elements->Nodes[INDEX2(4,k,8)]=node0+1;
674     out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0t+2;
675     out->Elements->Nodes[INDEX2(6,k,8)]=node0+2*N0t+1;
676     out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0t;
677     }
678     if( boundaryLeft )
679     out->Elements->Dom[i1*numElementsLocal] = ELEMENT_BOUNDARY;
680     if( boundaryRight )
681     out->Elements->Dom[(i1+1)*numElementsLocal-1] = ELEMENT_BOUNDARY;
682     }
683     out->Elements->minColor=0;
684     out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
685    
686     /* set the distribution information for the elements */
687     Finley_ElementFile_setDomainFlags( out->Elements );
688    
689     /* face elements: */
690     if (useElementsOnFace) {
691     NUMNODES=8;
692     } else {
693     NUMNODES=3;
694     }
695    
696     totalNECount=numElementsLocal*NE1;
697     faceNECount=0;
698     faceNEBoundary = totalNECount + numElementsInternal*(!periodic[1])*2 + NE1*(!periodic[0])*(domLeft+domRight);
699    
700     if (!periodic[1]) {
701     /* ** elements on boundary 010 (x2=0): */
702    
703     #pragma omp parallel for private(i0,k,node0)
704     for (i0=0;i0<numElementsLocal; i0++) {
705     k=i0+faceNECount;
706     kk=i0;
707     facePerm = useElementsOnFace ? face0 : face0nc;
708    
709     out->FaceElements->Id[k]=i0+totalNECount;
710     out->FaceElements->Tag[k]=10;
711     out->FaceElements->Color[k]=i0%2;
712     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
713    
714     for( j=0; j<NUMNODES; j++ )
715     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
716     }
717     if( boundaryLeft ){
718     out->FaceElements->Dom[faceNECount]=ELEMENT_BOUNDARY;
719     if( periodicLocal[0] )
720     out->FaceElements->Dom[faceNECount+1]=ELEMENT_BOUNDARY;
721     }
722     if( boundaryRight )
723     out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;
724     totalNECount+=numElementsLocal;
725     faceNECount+=numElementsLocal;
726    
727     /* ** elements on boundary 020 (x2=1): */
728    
729     #pragma omp parallel for private(i0,k,node0)
730     for (i0=0;i0<numElementsLocal;i0++) {
731     k=i0+faceNECount;
732     kk=i0 + numElementsLocal*(NE1-1);
733     facePerm = useElementsOnFace ? face1 : face1nc;
734    
735     out->FaceElements->Id[k]=i0+totalNECount;
736     out->FaceElements->Tag[k]=20;
737     out->FaceElements->Color[k]=i0%2+2;
738     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
739    
740     for( j=0; j<NUMNODES; j++ )
741     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
742     }
743     if( boundaryLeft ){
744     out->FaceElements->Dom[faceNECount]=ELEMENT_BOUNDARY;
745     if( periodicLocal[0] )
746     out->FaceElements->Dom[faceNECount+1]=ELEMENT_BOUNDARY;
747     }
748     if( boundaryRight )
749     out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;
750     totalNECount+=numElementsLocal;
751     faceNECount+=numElementsLocal;
752     }
753     if (domLeft && !periodicLocal[0]) {
754     /* ** elements on boundary 001 (x1=0): */
755    
756     #pragma omp parallel for private(i1,k,node0)
757     for (i1=0;i1<NE1;i1++) {
758     k=i1+faceNECount;
759     kk=i1*numElementsLocal;
760     facePerm = useElementsOnFace ? face2 : face2nc;
761    
762     out->FaceElements->Id[k]=i1+totalNECount;
763     out->FaceElements->Tag[k]=1;
764     out->FaceElements->Color[k]=(i1%2)+4;
765     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
766    
767     for( j=0; j<NUMNODES; j++ )
768     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
769     }
770     totalNECount+=NE1;
771     faceNECount+=NE1;
772    
773     /* ** elements on boundary 002 (x1=1): */
774     }
775     if (domRight && !periodicLocal[1]) {
776     #pragma omp parallel for private(i1,k,node0)
777     for (i1=0;i1<NE1;i1++) {
778     k=i1+faceNECount;
779     kk=(i1+1)*numElementsLocal-1;
780     facePerm = useElementsOnFace ? face3 : face3nc;
781    
782     out->FaceElements->Id[k]=i1+totalNECount;
783     out->FaceElements->Tag[k]=2;
784     out->FaceElements->Color[k]=(i1%2)+6;
785     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
786    
787     for( j=0; j<NUMNODES; j++ )
788     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
789     }
790     totalNECount+=NE1;
791     faceNECount+=NE1;
792     }
793     out->FaceElements->minColor=0;
794     out->FaceElements->maxColor=7;
795     out->FaceElements->numElements = faceNECount;
796    
797     Finley_ElementFile_setDomainFlags( out->FaceElements );
798    
799     /* setup distribution info for other elements */
800     Finley_ElementFile_setDomainFlags( out->ContactElements );
801     Finley_ElementFile_setDomainFlags( out->Points );
802    
803     /* reorder the degrees of freedom */
804     Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
805    
806     /* condense the nodes: */
807     Finley_Mesh_resolveNodeIds(out);
808     if ( !Finley_MPI_noError( mpi_info )) {
809     Paso_MPIInfo_dealloc( mpi_info );
810     Finley_Mesh_dealloc(out);
811     return NULL;
812     }
813    
814     /* prepare mesh for further calculatuions:*/
815     Finley_Mesh_prepare(out) ;
816    
817     #ifdef Finley_TRACE
818     printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
819     #endif
820    
821     if ( !Finley_MPI_noError( mpi_info )) {
822     Paso_MPIInfo_dealloc( mpi_info );
823     Finley_Mesh_dealloc(out);
824     return NULL;
825     }
826    
827     return out;
828     }
829     #endif

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26