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

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26