/[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 964 - (hide annotations)
Tue Feb 13 05:10:26 2007 UTC (12 years, 9 months ago) by gross
File MIME type: text/plain
File size: 29242 byte(s)
The set/getRefVal functions of Data objects have been removed (mainly to avoid later problems with MPI).
Moreover, a faster access to the reference id of samples has been introduced. I don't think that anybody will
profit form this at this stage but it will allow a faster dump of data objects.


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     index_t kk, k,i,j,node0,firstNodeConstruct,firstNode=0, DOFcount=0, forwardDOF=NULL, backwardDOF=NULL, targetDomain=-1,faceNEBoundary;
466     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