/[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 971 - (hide annotations)
Wed Feb 14 04:40:49 2007 UTC (12 years, 5 months ago) by ksteube
File MIME type: text/plain
File size: 29242 byte(s)
Had to undo commit to new MPI branch. The changes went into the original and
not the branch. The files committed here are exactly the same as revision 969.


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 gross 964 if (Finley_noError()) {
451     if ( ! Finley_Mesh_isPrepared(out)) {
452     Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");
453     }
454     } else {
455 jgs 82 Finley_Mesh_dealloc(out);
456     }
457     return out;
458     }
459 bcumming 782
460     #ifdef PASO_MPI
461     Finley_Mesh* Finley_RectangularMesh_Rec8(int* numElements,double* Length,int* periodic,int order,int useElementsOnFace) {
462     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];
463     dim_t N0t, NDOF0t;
464     index_t *numForward=NULL, *numBackward=NULL;
465 ksteube 971 index_t kk, k,i,j,node0,firstNodeConstruct,firstNode=0, DOFcount=0, forwardDOF=NULL, backwardDOF=NULL, targetDomain=-1,faceNEBoundary;
466 bcumming 782 Finley_Mesh* out=NULL;
467     char name[50];
468     index_t *indexForward=NULL, *indexBackward=NULL;
469     double time0=Finley_timer();
470     bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft=FALSE, boundaryRight=FALSE;
471     Paso_MPIInfo *mpi_info = NULL;
472 bcumming 787 index_t faceStart, numFaceLeft;
473 bcumming 782
474     /* these are used in constructing the face elements, and give the permutation
475     of the node order of a normal element for each face */
476     index_t *facePerm;
477     index_t face0[]={0, 1, 2, 3, 4, 5, 6, 7};
478     index_t face1[]={2, 3, 0, 1, 6, 7, 4, 5};
479     index_t face2[]={3, 0, 1, 2, 7, 4, 5, 6};
480     index_t face3[]={1, 2, 3, 0, 5, 6, 7, 4};
481     index_t face0nc[]={0, 1, 4};
482     index_t face1nc[]={2, 3, 6};
483     index_t face2nc[]={3, 0, 7};
484     index_t face3nc[]={1, 2, 5};
485    
486     NE0=MAX(1,numElements[0]);
487     NE1=MAX(1,numElements[1]);
488     N0=2*NE0+1;
489     N1=2*NE1+1;
490    
491     /* get MPI information */
492     mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
493     if (! Finley_noError()) {
494     Finley_Mesh_dealloc(out);
495     return NULL;
496     }
497    
498     // use the serial method for the single CPU case
499     if( mpi_info->size==1 ){
500     out = Finley_RectangularMesh_Rec8_singleCPU( numElements,Length,periodic,order,useElementsOnFace, mpi_info);
501     return out;
502     }
503    
504     if( mpi_info->rank==0 )
505     domLeft = TRUE;
506     if( mpi_info->rank==mpi_info->size-1 )
507     domRight = TRUE;
508     if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
509     domInternal = TRUE;
510    
511     /* dimensions of the local subdomain */
512     domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal, DOFBoundary );
513    
514     NFaceElements=0;
515     if (!periodic[0]) {
516     NDOF0=N0;
517     NFaceElements+=(domLeft + domRight)*NE1;
518     } else {
519     NDOF0=N0-1;
520     }
521     if (!periodic[1]) {
522     NDOF1=N1;
523     NFaceElements+=2*numElementsLocal;
524     } else {
525     NDOF1=N1-1;
526     }
527    
528     boundaryLeft = !domLeft || periodicLocal[0];
529     boundaryRight = !domRight || periodicLocal[1];
530     N0t = numNodesLocal + boundaryRight + boundaryLeft*2;
531     NDOF0t = numDOFLocal + boundaryRight + boundaryLeft*2;
532     firstNodeConstruct = firstNode - 2*boundaryLeft;
533     firstNodeConstruct = firstNodeConstruct<0 ? N0+firstNodeConstruct-1 : firstNodeConstruct;
534    
535     /* allocate mesh: */
536    
537     sprintf(name,"Rectangular %d x %d mesh",N0,N1);
538     out=Finley_Mesh_alloc(name,2,order,mpi_info);
539    
540     if (! Finley_noError()) return NULL;
541    
542     out->Elements=Finley_ElementFile_alloc(Rec8,out->order,mpi_info);
543     if (useElementsOnFace) {
544     out->FaceElements=Finley_ElementFile_alloc(Rec8Face,out->order,mpi_info);
545     out->ContactElements=Finley_ElementFile_alloc(Rec8Face_Contact,out->order,mpi_info);
546     } else {
547     out->FaceElements=Finley_ElementFile_alloc(Line3,out->order,mpi_info);
548     out->ContactElements=Finley_ElementFile_alloc(Line3_Contact,out->order,mpi_info);
549     }
550     out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
551    
552     if (! Finley_noError()) {
553     Finley_Mesh_dealloc(out);
554     return NULL;
555     }
556    
557     /* allocate tables: */
558     Finley_NodeFile_allocTable(out->Nodes,N0t*N1);
559     Finley_ElementFile_allocTable(out->Elements,numElementsLocal*NE1);
560     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
561 bcumming 787
562     Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0t*NDOF1, NDOF1*3, 0 );
563     Finley_ElementDistribution_allocTable( out->Elements->elementDistribution, numElementsLocal*NE1, NE1*(numElementsLocal-boundaryRight*(!periodic[1])) );
564     Finley_ElementDistribution_allocTable( out->FaceElements->elementDistribution, NFaceElements, NFaceElements-2*boundaryRight*(!periodic[1]) );
565 bcumming 782 if (! Finley_noError()) {
566     Finley_Mesh_dealloc(out);
567     return NULL;
568     }
569    
570     /* set nodes: */
571     #pragma omp parallel for private(i0,i1,k)
572 bcumming 787 for (i1=0;i1<N1;i1++) {
573 bcumming 782 for (i0=0;i0<N0t;i0++, k++) {
574 bcumming 787 k=i0+i1*N0t;
575 bcumming 782 out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE((firstNodeConstruct + i0) % N0)/DBLE(N0-1)*Length[0];
576     out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
577     out->Nodes->Id[k]=i0 + i1*N0t;
578     out->Nodes->Tag[k]=0;
579     out->Nodes->Dom[k]=NODE_INTERNAL;
580     out->Nodes->degreeOfFreedom[k]=i0 + (i1%NDOF1)*N0t;
581     }
582     /* take care to make the DOF reflect the internal/external position of the DOF */
583     if( boundaryLeft ){
584     out->Nodes->Dom[i1*N0t] = NODE_EXTERNAL;
585     out->Nodes->Dom[i1*N0t+1] = NODE_EXTERNAL;
586     out->Nodes->Dom[i1*N0t+2] = NODE_BOUNDARY;
587     if( periodicLocal[0] ){
588     out->Nodes->degreeOfFreedom[i1*N0t+3] = out->Nodes->degreeOfFreedom[i1*N0t+2];
589     out->Nodes->Dom[i1*N0t+3] = NODE_BOUNDARY;
590     }
591     }
592     if( boundaryRight ){
593     out->Nodes->Dom[(i1+1)*N0t-1] = NODE_EXTERNAL;
594     out->Nodes->Dom[(i1+1)*N0t-2] = NODE_BOUNDARY;
595     out->Nodes->Dom[(i1+1)*N0t-3] = NODE_BOUNDARY;
596     }
597     }
598     /* tags for the faces: */
599     if (!periodic[1]) {
600     for (i0=0;i0<N0t;i0++) {
601     out->Nodes->Tag[i0]+=10;
602     out->Nodes->Tag[i0+N0t*(N1-1)]+=20;
603     }
604     }
605     if (domLeft && !periodic[0]) {
606     for (i1=0;i1<N1;i1++) {
607     out->Nodes->Tag[i1*N0t]+=1;
608     }
609     }
610     if (domRight && !periodic[0]){
611     for (i1=0;i1<N1;i1++) {
612     out->Nodes->Tag[(i1+1)*N0t-1]+=2;
613     }
614     }
615    
616     /* setup the receive information */
617     indexBackward = TMPMEMALLOC( NDOF1*2, index_t );
618     indexForward = TMPMEMALLOC( NDOF1*2, index_t );
619     if( !(mpi_info->size==2 && periodicLocal[0])){
620     /* Backward (receive) DOF */
621     if( boundaryLeft ){
622     targetDomain = mpi_info->rank-1<0 ? mpi_info->size-1 : mpi_info->rank-1;
623    
624     for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){
625     indexBackward[i1] = out->Nodes->degreeOfFreedom[i*N0t];
626     indexBackward[i1+1] = out->Nodes->degreeOfFreedom[i*N0t+1];
627     indexForward[i] = out->Nodes->degreeOfFreedom[i*N0t+2];
628     }
629    
630     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );
631     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexBackward );
632     }
633    
634     if( boundaryRight ){
635     targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;
636    
637     for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){
638     indexBackward[i] = out->Nodes->degreeOfFreedom[(i+1)*N0t-1];
639     indexForward[i1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-3];
640     indexForward[i1+1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-2];
641     }
642    
643     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexForward );
644     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexBackward );
645     }
646     }
647     else{
648     targetDomain = 1;
649    
650     for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){
651     indexBackward[i] = out->Nodes->degreeOfFreedom[(i+1)*N0t-1];
652     indexForward[i1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-3];
653     indexForward[i1+1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-2];
654     }
655    
656     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexForward );
657     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexBackward );
658    
659     for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){
660     indexBackward[i1] = out->Nodes->degreeOfFreedom[i*N0t];
661     indexBackward[i1+1] = out->Nodes->degreeOfFreedom[i*N0t+1];
662     indexForward[i] = out->Nodes->degreeOfFreedom[i*N0t+2];
663     }
664    
665     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );
666     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexBackward );
667     }
668     TMPMEMFREE( indexBackward );
669     TMPMEMFREE( indexForward );
670    
671     /* set the elements: */
672     #pragma omp parallel for private(i0,i1,k,node0)
673     for (i1=0;i1<NE1;i1++) {
674     for (i0=0;i0<numElementsLocal;i0++) {
675     k=i0+numElementsLocal*i1;
676     node0 = (periodicLocal[0] && !i0) ? 2*i1*N0t : 2*i0+2*i1*N0t+periodicLocal[0] ;
677    
678 bcumming 787 out->Elements->Id[k]=((firstNodeConstruct/2+i0)%NE0) * NE1 + i1;
679 bcumming 782 out->Elements->Tag[k]=0;
680     out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);
681     out->Elements->Dom[k] = ELEMENT_INTERNAL;
682    
683     out->Elements->Nodes[INDEX2(0,k,8)]=node0;
684     out->Elements->Nodes[INDEX2(1,k,8)]=node0+2;
685     out->Elements->Nodes[INDEX2(2,k,8)]=node0+2*N0t+2;
686     out->Elements->Nodes[INDEX2(3,k,8)]=node0+2*N0t;
687     out->Elements->Nodes[INDEX2(4,k,8)]=node0+1;
688     out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0t+2;
689     out->Elements->Nodes[INDEX2(6,k,8)]=node0+2*N0t+1;
690     out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0t;
691     }
692     if( boundaryLeft )
693     out->Elements->Dom[i1*numElementsLocal] = ELEMENT_BOUNDARY;
694     if( boundaryRight )
695     out->Elements->Dom[(i1+1)*numElementsLocal-1] = ELEMENT_BOUNDARY;
696     }
697     out->Elements->minColor=0;
698     out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
699    
700     /* set the distribution information for the elements */
701     Finley_ElementFile_setDomainFlags( out->Elements );
702    
703     /* face elements: */
704     if (useElementsOnFace) {
705     NUMNODES=8;
706     } else {
707     NUMNODES=3;
708     }
709    
710 bcumming 787 Finley_ElementFile_setDomainFlags( out->FaceElements );
711 bcumming 782 totalNECount=numElementsLocal*NE1;
712     faceNECount=0;
713     faceNEBoundary = totalNECount + numElementsInternal*(!periodic[1])*2 + NE1*(!periodic[0])*(domLeft+domRight);
714    
715 bcumming 787 numFaceLeft = domLeft*(!periodic[0])*NE1;
716     faceStart = out->FaceElements->elementDistribution->vtxdist[mpi_info->rank];
717 bcumming 782 if (!periodic[1]) {
718 bcumming 787 /* ** elements on boundary 010 (x1=0): */
719 bcumming 782
720 bcumming 787 #pragma omp parallel for private(i0,k,kk,facePerm)
721 bcumming 782 for (i0=0;i0<numElementsLocal; i0++) {
722     k=i0+faceNECount;
723     kk=i0;
724     facePerm = useElementsOnFace ? face0 : face0nc;
725    
726 bcumming 787 out->FaceElements->Id[k]=faceStart + numFaceLeft + i0*2;
727 bcumming 782 out->FaceElements->Tag[k]=10;
728     out->FaceElements->Color[k]=i0%2;
729     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
730    
731     for( j=0; j<NUMNODES; j++ )
732     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
733     }
734     if( boundaryLeft ){
735     out->FaceElements->Dom[faceNECount]=ELEMENT_BOUNDARY;
736     if( periodicLocal[0] )
737     out->FaceElements->Dom[faceNECount+1]=ELEMENT_BOUNDARY;
738     }
739 bcumming 787 if( boundaryRight ){
740 bcumming 782 out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;
741 bcumming 787 out->FaceElements->Id[faceNECount+numElementsLocal-1]=out->FaceElements->elementDistribution->vtxdist[mpi_info->rank+1];
742     }
743 bcumming 782 totalNECount+=numElementsLocal;
744     faceNECount+=numElementsLocal;
745    
746 bcumming 787 /* ** elements on boundary 020 (x1=1): */
747 bcumming 782
748 bcumming 787 #pragma omp parallel for private(i0,k,kk,facePerm)
749 bcumming 782 for (i0=0;i0<numElementsLocal;i0++) {
750     k=i0+faceNECount;
751     kk=i0 + numElementsLocal*(NE1-1);
752     facePerm = useElementsOnFace ? face1 : face1nc;
753    
754 bcumming 787 out->FaceElements->Id[k]=faceStart + numFaceLeft + i0*2+1;
755 bcumming 782 out->FaceElements->Tag[k]=20;
756     out->FaceElements->Color[k]=i0%2+2;
757     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
758    
759     for( j=0; j<NUMNODES; j++ )
760     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
761     }
762     if( boundaryLeft ){
763     out->FaceElements->Dom[faceNECount]=ELEMENT_BOUNDARY;
764     if( periodicLocal[0] )
765     out->FaceElements->Dom[faceNECount+1]=ELEMENT_BOUNDARY;
766     }
767 bcumming 787 if( boundaryRight ){
768 bcumming 782 out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;
769 bcumming 787 out->FaceElements->Id[faceNECount+numElementsLocal-1]=out->FaceElements->elementDistribution->vtxdist[mpi_info->rank+1]+1;
770     }
771 bcumming 782 totalNECount+=numElementsLocal;
772     faceNECount+=numElementsLocal;
773     }
774     if (domLeft && !periodicLocal[0]) {
775 bcumming 787 /* ** elements on boundary 001 (x0=0): */
776 bcumming 782
777 bcumming 787 #pragma omp parallel for private(i1,k,kk,facePerm)
778 bcumming 782 for (i1=0;i1<NE1;i1++) {
779     k=i1+faceNECount;
780     kk=i1*numElementsLocal;
781     facePerm = useElementsOnFace ? face2 : face2nc;
782    
783 bcumming 787 out->FaceElements->Id[k]=faceStart + i1;
784 bcumming 782 out->FaceElements->Tag[k]=1;
785     out->FaceElements->Color[k]=(i1%2)+4;
786     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
787    
788     for( j=0; j<NUMNODES; j++ )
789     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
790     }
791     totalNECount+=NE1;
792     faceNECount+=NE1;
793    
794 bcumming 787 /* ** elements on boundary 002 (x0=1): */
795 bcumming 782 }
796     if (domRight && !periodicLocal[1]) {
797 bcumming 787 #pragma omp parallel for private(i1,k,kk,facePerm)
798 bcumming 782 for (i1=0;i1<NE1;i1++) {
799     k=i1+faceNECount;
800     kk=(i1+1)*numElementsLocal-1;
801     facePerm = useElementsOnFace ? face3 : face3nc;
802    
803 bcumming 787 out->FaceElements->Id[k]=faceStart + numElementsLocal*2*(!periodic[1])+i1;
804 bcumming 782 out->FaceElements->Tag[k]=2;
805     out->FaceElements->Color[k]=(i1%2)+6;
806     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
807    
808     for( j=0; j<NUMNODES; j++ )
809     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
810     }
811     totalNECount+=NE1;
812     faceNECount+=NE1;
813     }
814     out->FaceElements->minColor=0;
815     out->FaceElements->maxColor=7;
816     out->FaceElements->numElements = faceNECount;
817    
818     /* setup distribution info for other elements */
819     Finley_ElementFile_setDomainFlags( out->ContactElements );
820     Finley_ElementFile_setDomainFlags( out->Points );
821    
822 bcumming 787 Finley_Mesh_prepareElementDistribution( out );
823    
824 bcumming 782 /* reorder the degrees of freedom */
825     Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
826    
827     /* condense the nodes: */
828     Finley_Mesh_resolveNodeIds(out);
829     if ( !Finley_MPI_noError( mpi_info )) {
830     Paso_MPIInfo_dealloc( mpi_info );
831     Finley_Mesh_dealloc(out);
832     return NULL;
833     }
834    
835     /* prepare mesh for further calculatuions:*/
836     Finley_Mesh_prepare(out) ;
837    
838     #ifdef Finley_TRACE
839     printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
840     #endif
841    
842     if ( !Finley_MPI_noError( mpi_info )) {
843     Paso_MPIInfo_dealloc( mpi_info );
844     Finley_Mesh_dealloc(out);
845     return NULL;
846     }
847 gross 964 if (Finley_noError()) {
848     if (! Finley_Mesh_isPrepared(out)) {
849     Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");
850     }
851     }
852     return out;
853 bcumming 782 }
854     #endif

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26