/[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 751 - (hide annotations)
Mon Jun 26 01:46:34 2006 UTC (13 years, 4 months ago) by bcumming
File MIME type: text/plain
File size: 34191 byte(s)
Changes relating to the MPI version of escript
The standard OpenMP version of escript is unchanged

- updated data types (Finley_Mesh, Finley_NodeFile, etc) to store meshes
  over multiple MPI processes.
- added CommBuffer code in Paso for communication of Data associated
  with distributed meshes
- updates in Finley and Escript to support distributed data and operations
  on distributed data (such as interpolation).
- construction of RHS in MPI, so that simple explicit schemes (such as
  /docs/examples/wave.py without IO and the Locator) can run in MPI.
- updated mesh generation for first order line, rectangle and brick
  meshes and second order line meshes in MPI.        
- small changes to trunk/SConstruct and trunk/scons/ess_options.py to
  build the MPI version, these changes are turned off by default.

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     /* A bit messy, but it only has to be done once... */
47 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 )
48 bcumming 730 {
49     index_t i0;
50     dim_t numNodesGlobal = numElementsGlobal+1;
51    
52     (*numNodesLocal) = domain_MODdim( rank, size, numNodesGlobal );
53 bcumming 751
54 bcumming 730 numElementsLocal[0] = numNodesLocal[0]+1;
55     periodicLocal[0] = periodicLocal[1] = FALSE;
56     nodesExternal[0] = nodesExternal[1] = 1;
57     if( periodic )
58     {
59     if( size==1 )
60     {
61     numElementsLocal[0] = numElementsGlobal;
62     nodesExternal[0] = nodesExternal[1] = 0;
63     periodicLocal[0] = periodicLocal[1] = TRUE;
64     }
65     else
66     {
67     if( rank==0 )
68     {
69     periodicLocal[0] = TRUE;
70     numNodesLocal[0]++;
71     }
72     if( rank==(size-1) )
73     {
74     periodicLocal[1] = TRUE;
75     numNodesLocal[0]--;
76     numElementsLocal[0]--;
77     }
78     }
79     }
80     else if( !periodic )
81     {
82     if( rank==0 ){
83     nodesExternal[0]--;
84     numElementsLocal[0]--;
85     }
86     if( rank==(size-1) )
87     {
88     nodesExternal[1]--;
89     numElementsLocal[0]--;
90     }
91     }
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    
106     numDOFLocal[0] = numNodesLocal[0];
107     if( periodicLocal[0] )
108     numDOFLocal[0]--;
109 bcumming 751
110     numDOFInternal[0] = numDOFLocal[0];
111     if( size>1 )
112     {
113     if( rank>0 || periodic ){
114     numDOFInternal[0] -= 1;
115     }
116     if( rank<(size-1) || periodic ){
117     numDOFInternal[0] -= 1;
118     }
119     }
120 bcumming 730 DOFExternal[0] = nodesExternal[0];
121     DOFExternal[1] = nodesExternal[1];
122     }
123 bcumming 751 #endif
124 bcumming 730
125     Finley_Mesh* Finley_RectangularMesh_Rec4(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
126     #ifndef PASO_MPI
127     {
128 jgs 153 dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES,M0,M1;
129 jgs 123 index_t k,node0;
130 jgs 82 Finley_Mesh* out;
131     char name[50];
132     double time0=Finley_timer();
133     NE0=MAX(1,numElements[0]);
134     NE1=MAX(1,numElements[1]);
135     N0=NE0+1;
136     N1=NE1+1;
137    
138 jgs 153 if (N0<=N1) {
139     M0=1;
140     M1=N0;
141     } else {
142     M0=N1;
143     M1=1;
144     }
145    
146 jgs 82 NFaceElements=0;
147     if (!periodic[0]) {
148     NDOF0=N0;
149     NFaceElements+=2*NE1;
150     } else {
151     NDOF0=N0-1;
152     }
153     if (!periodic[1]) {
154     NDOF1=N1;
155     NFaceElements+=2*NE0;
156     } else {
157     NDOF1=N1-1;
158     }
159    
160     /* allocate mesh: */
161    
162     sprintf(name,"Rectangular %d x %d mesh",N0,N1);
163     out=Finley_Mesh_alloc(name,2,order);
164 bcumming 730
165 jgs 150 if (! Finley_noError()) return NULL;
166 jgs 82
167     out->Elements=Finley_ElementFile_alloc(Rec4,out->order);
168     if (useElementsOnFace) {
169     out->FaceElements=Finley_ElementFile_alloc(Rec4Face,out->order);
170     out->ContactElements=Finley_ElementFile_alloc(Rec4Face_Contact,out->order);
171     } else {
172     out->FaceElements=Finley_ElementFile_alloc(Line2,out->order);
173     out->ContactElements=Finley_ElementFile_alloc(Line2_Contact,out->order);
174     }
175     out->Points=Finley_ElementFile_alloc(Point1,out->order);
176 bcumming 730
177 jgs 150 if (! Finley_noError()) {
178 jgs 82 Finley_Mesh_dealloc(out);
179     return NULL;
180     }
181    
182     /* allocate tables: */
183     Finley_NodeFile_allocTable(out->Nodes,N0*N1);
184     Finley_ElementFile_allocTable(out->Elements,NE0*NE1);
185     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
186 bcumming 730
187 jgs 150 if (! Finley_noError()) {
188 jgs 82 Finley_Mesh_dealloc(out);
189     return NULL;
190     }
191     /* set nodes: */
192 jgs 153
193     #pragma omp parallel for private(i0,i1,k)
194 jgs 82 for (i1=0;i1<N1;i1++) {
195     for (i0=0;i0<N0;i0++) {
196 jgs 153 k=M0*i0+M1*i1;
197 jgs 82 out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(i0)/DBLE(N0-1)*Length[0];
198     out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
199 jgs 153 out->Nodes->Id[k]=i0+N0*i1;
200 jgs 82 out->Nodes->Tag[k]=0;
201 jgs 153 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1);
202 jgs 82 }
203     }
204     /* tags for the faces: */
205     if (!periodic[1]) {
206     for (i0=0;i0<N0;i0++) {
207 jgs 153 out->Nodes->Tag[M0*i0+M1*0]+=10;
208     out->Nodes->Tag[M0*i0+M1*(N1-1)]+=20;
209 jgs 82 }
210     }
211     if (!periodic[0]) {
212     for (i1=0;i1<N1;i1++) {
213 jgs 153 out->Nodes->Tag[M0*0+M1*i1]+=1;
214     out->Nodes->Tag[M0*(N0-1)+M1*i1]+=2;
215 jgs 82 }
216     }
217 jgs 153
218 jgs 82 /* set the elements: */
219    
220     #pragma omp parallel for private(i0,i1,k,node0)
221     for (i1=0;i1<NE1;i1++) {
222     for (i0=0;i0<NE0;i0++) {
223     k=i0+NE0*i1;
224     node0=i0+i1*N0;
225    
226     out->Elements->Id[k]=k;
227     out->Elements->Tag[k]=0;
228     out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);
229    
230     out->Elements->Nodes[INDEX2(0,k,4)]=node0;
231     out->Elements->Nodes[INDEX2(1,k,4)]=node0+1;
232     out->Elements->Nodes[INDEX2(2,k,4)]=node0+N0+1;
233     out->Elements->Nodes[INDEX2(3,k,4)]=node0+N0;
234    
235     }
236     }
237 jgs 123 out->Elements->minColor=0;
238     out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
239 jgs 82
240     /* face elements: */
241    
242     if (useElementsOnFace) {
243     NUMNODES=4;
244     } else {
245     NUMNODES=2;
246     }
247     totalNECount=NE0*NE1;
248     faceNECount=0;
249    
250     /* elements on boundary 010 (x2=0): */
251     if (!periodic[1]) {
252     #pragma omp parallel for private(i0,k,node0)
253     for (i0=0;i0<NE0;i0++) {
254     k=i0+faceNECount;
255     node0=i0;
256    
257     out->FaceElements->Id[k]=i0+totalNECount;
258     out->FaceElements->Tag[k]=10;
259     out->FaceElements->Color[k]=i0%2;
260    
261     if (useElementsOnFace) {
262     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
263     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
264     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
265     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
266     } else {
267     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
268     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
269     }
270     }
271     totalNECount+=NE0;
272     faceNECount+=NE0;
273    
274     /* ** elements on boundary 020 (x2=1): */
275    
276     #pragma omp parallel for private(i0,k,node0)
277     for (i0=0;i0<NE0;i0++) {
278     k=i0+faceNECount;
279     node0=i0+(NE1-1)*N0;
280    
281     out->FaceElements->Id[k]=i0+totalNECount;
282     out->FaceElements->Tag[k]=20;
283     out->FaceElements->Color[k]=i0%2+2;
284    
285     if (useElementsOnFace) {
286     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0+1;
287     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
288     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;
289     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
290     } else {
291     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0+1;
292     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
293     }
294     }
295     totalNECount+=NE0;
296     faceNECount+=NE0;
297     }
298     if (!periodic[0]) {
299     /* ** elements on boundary 001 (x1=0): */
300     #pragma omp parallel for private(i1,k,node0)
301     for (i1=0;i1<NE1;i1++) {
302     k=i1+faceNECount;
303     node0=i1*N0;
304    
305     out->FaceElements->Id[k]=i1+totalNECount;
306     out->FaceElements->Tag[k]=1;
307     out->FaceElements->Color[k]=i1%2+4;
308    
309     if (useElementsOnFace) {
310     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
311     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
312     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
313     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
314     } else {
315     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
316     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
317    
318     }
319     }
320     totalNECount+=NE1;
321     faceNECount+=NE1;
322    
323     /* ** elements on boundary 002 (x1=1): */
324    
325     #pragma omp parallel for private(i1,k,node0)
326     for (i1=0;i1<NE1;i1++) {
327     k=i1+faceNECount;
328     node0=(NE0-1)+i1*N0;
329    
330     out->FaceElements->Id[k]=i1+totalNECount;
331     out->FaceElements->Tag[k]=2;
332     out->FaceElements->Color[k]=i1%2+6;
333    
334     if (useElementsOnFace) {
335     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
336     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
337     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0;
338     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
339     } else {
340     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
341     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
342     }
343     }
344     totalNECount+=NE1;
345     faceNECount+=NE1;
346     }
347 jgs 123 out->FaceElements->minColor=0;
348     out->FaceElements->maxColor=7;
349 jgs 82
350 bcumming 730
351 jgs 82 /* face elements done: */
352    
353     /* condense the nodes: */
354    
355     Finley_Mesh_resolveNodeIds(out);
356    
357     /* prepare mesh for further calculatuions:*/
358    
359     Finley_Mesh_prepare(out) ;
360    
361 jgs 149 #ifdef Finley_TRACE
362 jgs 82 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
363 jgs 149 #endif
364 jgs 82
365 jgs 150 if (! Finley_noError()) {
366 jgs 82 Finley_Mesh_dealloc(out);
367     return NULL;
368     }
369     return out;
370     }
371 bcumming 730
372     /*----------------------------------------------------------------------------
373     MPI VERSION
374    
375    
376     ASSUMPTIONS
377     ===========
378    
379     - the domain dimension is large enough in the NE0 direction for each domain
380     to be 2 internal nodes wide in that direction. Face element calculation gets
381     buggered up otherwise.
382     - if dividing into two domains with periodic[0]=TRUE , then the global domain
383     must be at least 4 elements along the NE0 direction, otherwise redundant
384     nodes are generated.
385    
386     These problems can be avoided by ensuring that numElements[0]/mpiSize >= 2
387    
388     if domains are small enough to cause these problems, you probably don't
389     need to solve them in parallel. If you really do, rotate the domain so that the
390     long axis is aligned with NE0.
391     ------------------------------------------------------------------------------*/
392    
393    
394     #else
395     {
396 bcumming 751 dim_t N0,N1,NE0,NE1,i0,i1, NE0_local, NDOF0,NDOF1, NFaceElements, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal, faceNECount, totalNECount,M0,M1,numDOFInternal;
397 bcumming 730 index_t innerLoop=0;
398     index_t NUMNODES,k,firstNode=0, DOFcount=0, forwardDOF[2], backwardDOF[2], node0, node1;
399     index_t targetDomain=-1;
400     index_t *indexForward=NULL;
401     bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE;
402     Finley_Mesh* out;
403     char name[50];
404     double time0=Finley_timer();
405    
406     NE0=MAX(1,numElements[0]);
407     NE1=MAX(1,numElements[1]);
408     N0=NE0+1;
409     N1=NE1+1;
410     Paso_MPIInfo *mpi_info = NULL;
411    
412     /* get MPI information */
413     mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
414     if (! Finley_noError())
415     return NULL;
416    
417     if( mpi_info->rank==0 )
418     domLeft = TRUE;
419     if( mpi_info->rank==mpi_info->size-1 )
420     domRight = TRUE;
421     if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
422     domInternal = TRUE;
423    
424     /* dimensions of the local subdomain */
425 bcumming 751 domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal, &numDOFInternal );
426 bcumming 730
427     if (N0<=N1) {
428     M0=1;
429     M1=N0;
430     } else {
431     M0=N1;
432     M1=1;
433     }
434    
435     NFaceElements=0;
436     if (!periodic[0])
437     {
438     if( domLeft )
439     NFaceElements+=NE1;
440     if( domRight )
441     NFaceElements+=NE1;
442     }
443     if (!periodic[1]) {
444     NDOF1=N1;
445 bcumming 751 NFaceElements+=2*numElementsLocal;
446 bcumming 730 } else {
447     NDOF1=N1-1;
448     }
449    
450     /* allocate mesh: */
451    
452     sprintf(name,"Rectangular %d x %d mesh",N0,N1);
453     out=Finley_Mesh_alloc( name, 2, order, mpi_info );
454    
455     if (! Finley_noError()) return NULL;
456    
457     out->Elements=Finley_ElementFile_alloc( Rec4, out->order, mpi_info );
458     if (useElementsOnFace) {
459     out->FaceElements=Finley_ElementFile_alloc(Rec4Face,out->order, mpi_info );
460     out->ContactElements=Finley_ElementFile_alloc(Rec4Face_Contact,out->order, mpi_info );
461     } else {
462     out->FaceElements=Finley_ElementFile_alloc(Line2,out->order, mpi_info );
463     out->ContactElements=Finley_ElementFile_alloc(Line2_Contact,out->order, mpi_info );
464     }
465     out->Points=Finley_ElementFile_alloc(Point1,out->order, mpi_info );
466     if (! Finley_noError()) {
467     Finley_Mesh_dealloc(out);
468     return NULL;
469     }
470    
471     /* allocate tables: */
472     Finley_NodeFile_allocTable( out->Nodes, (numNodesLocal + numNodesExternal)*N1 );
473 bcumming 751 Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal*NDOF1, (DOFExternal[0]+DOFExternal[1])*NDOF1, 0 );
474 bcumming 730 Finley_ElementFile_allocTable(out->Elements,numElementsLocal*NE1);
475     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
476    
477     if (! Finley_noError()) {
478     Finley_Mesh_dealloc(out);
479     return NULL;
480     }
481     /* set nodes: */
482    
483     //#pragma omp parallel for private(i0,i1,k)
484     if( mpi_info->size==1 )
485     innerLoop = numNodesLocal;
486     else
487     innerLoop = numNodesLocal-periodicLocal[0];
488    
489     for (k=0,i1=0;i1<N1;i1++) {
490     for ( i0=0;i0<innerLoop;i0++,k++) {
491     out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(i0+firstNode)/DBLE(N0-1)*Length[0];
492     out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
493     out->Nodes->Id[k]=i0+innerLoop*i1;
494     out->Nodes->Tag[k]=0;
495     out->Nodes->degreeOfFreedom[k]=i0%numDOFLocal + numDOFLocal*(i1%NDOF1);
496     }
497     }
498    
499     /* tag 'em */
500     for( i0=0; i0<numNodesLocal; i0++ )
501     {
502     out->Nodes->Tag[i0] += 10; // bottom
503     out->Nodes->Tag[i0+numNodesLocal*(N1-1)] += 20; // top
504     }
505     if( domLeft && !periodicLocal[0] )
506     for( i1=0; i1<N1; i1++ )
507     out->Nodes->Tag[i1*numNodesLocal] += 1; // left
508    
509     if( domRight && !periodicLocal[0] )
510     for( i1=0; i1<N1; i1++ )
511     out->Nodes->Tag[(i1+1)*numNodesLocal-1] += 2; // right
512    
513     /* external nodes */
514     k = innerLoop*N1;
515     DOFcount = out->Nodes->degreeOfFreedom[k-1]+1;
516     if( mpi_info->size>1 )
517     {
518     /* add forward communication information for boundary nodes */
519     indexForward = TMPMEMALLOC( NDOF1, index_t );
520     if( domInternal || domRight || periodicLocal[0] )
521     {
522     for( int i=0; i<NDOF1; i++ )
523     indexForward[i] = i*numDOFLocal;
524     targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;
525     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );
526     }
527     if( domInternal || domLeft || periodicLocal[1] )
528     {
529     for( int i=0; i<NDOF1; i++ )
530     indexForward[i] = (i+1)*numDOFLocal - 1;
531     targetDomain = (mpi_info->rank+1) % mpi_info->size;
532     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );
533     }
534     TMPMEMFREE( indexForward );
535    
536     /* LHS */
537     if( periodicLocal[0] )
538     {
539     for (i1=0; i1<N1; i1++, k++)
540     {
541     out->Nodes->Coordinates[INDEX2(0,k,2)]=Length[0];
542     out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
543     out->Nodes->Id[k]=k;
544     out->Nodes->Tag[k]=0;
545     out->Nodes->degreeOfFreedom[k] = (i1 % NDOF1)*numDOFLocal;
546     }
547     out->Nodes->Tag[k-N1] = 10;
548     out->Nodes->Tag[k-1] = 20;
549     for (i1=0; i1<N1; i1++, k++)
550     {
551     out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(NE0-1)/DBLE(NE0)*Length[0];
552     out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
553     out->Nodes->Id[k]=k;
554     out->Nodes->Tag[k]=0;
555     out->Nodes->degreeOfFreedom[k] = DOFcount + i1%NDOF1;
556     }
557     out->Nodes->Tag[k-N1] = 10;
558     out->Nodes->Tag[k-1] = 20;
559     DOFcount = out->Nodes->degreeOfFreedom[k-1]+1;
560    
561     targetDomain = mpi_info->size-1;
562     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, out->Nodes->degreeOfFreedom+(numNodesLocal*N1) );
563     }
564     if( mpi_info->rank>0 )
565     {
566     for (i1=0; i1<N1; i1++, k++)
567     {
568     out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(firstNode-1)/DBLE(N0-1)*Length[0];
569     out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
570     out->Nodes->Id[k]=k;
571     out->Nodes->Tag[k]=0;
572     out->Nodes->degreeOfFreedom[k] = DOFcount + i1%NDOF1;
573     }
574     out->Nodes->Tag[k-N1] = 10;
575     out->Nodes->Tag[k-1] = 20;
576     DOFcount = out->Nodes->degreeOfFreedom[k-1]+1;
577    
578     targetDomain = mpi_info->rank-1;
579     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, out->Nodes->degreeOfFreedom+(k-N1) );
580     }
581    
582     /* RHS */
583     if( periodicLocal[1] || (mpi_info->rank < mpi_info->size-1) )
584     {
585     for (i1=0; i1<N1; i1++, k++)
586     {
587     out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(firstNode+numNodesLocal-periodicLocal[0])/DBLE(N0-1)*Length[0];
588     out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
589     out->Nodes->Id[k]=k;
590     out->Nodes->Tag[k]=0;
591     out->Nodes->degreeOfFreedom[k] = DOFcount + i1%NDOF1;
592     }
593     out->Nodes->Tag[k-N1] = 10;
594     out->Nodes->Tag[k-1] = 20;
595     DOFcount = out->Nodes->degreeOfFreedom[k-1]+1;
596    
597     targetDomain = (mpi_info->rank+1) % mpi_info->size;
598     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, out->Nodes->degreeOfFreedom+(k-N1) );
599     }
600     }
601 bcumming 751 out->Nodes->degreeOfFreedomDistribution->numInternal = NDOF1*numDOFInternal;
602     out->Nodes->degreeOfFreedomDistribution->numBoundary = out->Nodes->degreeOfFreedomDistribution->numLocal- out->Nodes->degreeOfFreedomDistribution->numInternal;
603 bcumming 730 /* set the elements: */
604    
605     /* internal */
606     //#pragma omp parallel for private(i0,i1,k,node0)
607     for ( i1=0, k=0; i1<NE1; i1++)
608     {
609     for ( i0=0; i0<numElementsInternal; i0++, k++)
610     {
611     node0=i0+i1*innerLoop;
612    
613     out->Elements->Id[k]=k;
614     out->Elements->Tag[k]=0;
615     out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
616    
617     out->Elements->Nodes[INDEX2(0,k,4)]=node0;
618     out->Elements->Nodes[INDEX2(1,k,4)]=node0+1;
619     out->Elements->Nodes[INDEX2(2,k,4)]=node0+innerLoop+1;
620     out->Elements->Nodes[INDEX2(3,k,4)]=node0+innerLoop;
621    
622     }
623     }
624     /* boundary */
625     if( mpi_info->size>1 )
626     {
627     if( domInternal )
628     {
629     /* left */
630     for( i1=0; i1<NE1; i1++, k++ )
631     {
632     node1=numNodesLocal*N1 + i1;
633     node0=i1*numNodesLocal;
634    
635     out->Elements->Id[k]=k;
636     out->Elements->Tag[k]=0;
637     out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
638    
639     out->Elements->Nodes[INDEX2(0,k,4)]=node1;
640     out->Elements->Nodes[INDEX2(1,k,4)]=node0;
641     out->Elements->Nodes[INDEX2(2,k,4)]=node0+numNodesLocal;
642     out->Elements->Nodes[INDEX2(3,k,4)]=node1+1;
643     }
644    
645     /* right */
646     for( i1=0; i1<NE1; i1++, k++ )
647     {
648     node0 = (i1+1)*numNodesLocal-1;
649     node1 = (numNodesLocal+1)*N1 + i1;
650    
651     out->Elements->Id[k]=k;
652     out->Elements->Tag[k]=0;
653     out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
654    
655     out->Elements->Nodes[INDEX2(0,k,4)]=node0;
656     out->Elements->Nodes[INDEX2(1,k,4)]=node1;
657     out->Elements->Nodes[INDEX2(2,k,4)]=node1+1;
658     out->Elements->Nodes[INDEX2(3,k,4)]=node0+numNodesLocal;
659     }
660     }
661     else if( domLeft )
662     {
663     /* left */
664     if( periodicLocal[0] )
665     {
666     node1 = numNodesLocal*N1;
667     node0 = innerLoop*N1;
668     for( i1=0; i1<NE1; i1++, k++, node1++, node0++ )
669     {
670     out->Elements->Id[k]=k;
671     out->Elements->Tag[k]=0;
672     out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
673    
674     out->Elements->Nodes[INDEX2(0,k,4)]=node1;
675     out->Elements->Nodes[INDEX2(1,k,4)]=node0;
676     out->Elements->Nodes[INDEX2(2,k,4)]=node0+1;
677     out->Elements->Nodes[INDEX2(3,k,4)]=node1+1;
678     }
679     }
680     /* right */
681     for( i1=0; i1<NE1; i1++, k++ )
682     {
683     node0 = (i1+1)*innerLoop-1;
684     node1 = (numNodesLocal+periodicLocal[0])*N1 + i1;
685    
686     out->Elements->Id[k]=k;
687     out->Elements->Tag[k]=0;
688     out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
689    
690     out->Elements->Nodes[INDEX2(0,k,4)]=node0;
691     out->Elements->Nodes[INDEX2(1,k,4)]=node1;
692     out->Elements->Nodes[INDEX2(2,k,4)]=node1+1;
693     out->Elements->Nodes[INDEX2(3,k,4)]=node0+numNodesLocal-periodicLocal[0];
694     }
695     }
696     else if( domRight )
697     {
698     /* left */
699     for( i1=0; i1<NE1; i1++, k++ )
700     {
701     node1=numNodesLocal*N1 + i1;
702     node0=i1*numNodesLocal;
703    
704     out->Elements->Id[k]=k;
705     out->Elements->Tag[k]=0;
706     out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
707    
708     out->Elements->Nodes[INDEX2(0,k,4)]=node1;
709     out->Elements->Nodes[INDEX2(1,k,4)]=node0;
710     out->Elements->Nodes[INDEX2(2,k,4)]=node0+numNodesLocal;
711     out->Elements->Nodes[INDEX2(3,k,4)]=node1+1;
712     }
713     /* right */
714     if( periodicLocal[1] )
715     {
716     for( i1=0; i1<NE1; i1++, k++ )
717     {
718     node0=(i1+1)*numNodesLocal-1;
719     node1 = (numNodesLocal+1)*N1 + i1;
720    
721     out->Elements->Id[k]=k;
722     out->Elements->Tag[k]=0;
723     out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
724    
725     out->Elements->Nodes[INDEX2(0,k,4)]=node0;
726     out->Elements->Nodes[INDEX2(1,k,4)]=node1;
727     out->Elements->Nodes[INDEX2(2,k,4)]=node1+1;
728     out->Elements->Nodes[INDEX2(3,k,4)]=node0+numNodesLocal;
729     }
730     }
731     }
732     }
733     out->Elements->minColor=0;
734     out->Elements->maxColor=0; //COLOR_MOD(0)+3*COLOR_MOD(0);
735    
736     if( domInternal )
737     {
738     out->Elements->elementDistribution->numInternal = NE1*(numElementsLocal-2);
739     out->Elements->elementDistribution->numBoundary = NE1*2;
740     }
741     else if( mpi_info->size==1 )
742     {
743     out->Elements->elementDistribution->numInternal = NE1*numElementsLocal;
744     }
745     else
746     {
747     out->Elements->elementDistribution->numInternal += NE1*(numElementsLocal-1-periodic[0]);
748     out->Elements->elementDistribution->numBoundary += NE1*(1 + periodic[0]);
749     }
750    
751     out->FaceElements->elementDistribution->numLocal = out->FaceElements->elementDistribution->numInternal + out->FaceElements->elementDistribution->numBoundary;
752    
753     /* face elements: */
754    
755     if (useElementsOnFace) {
756     NUMNODES=4;
757     } else {
758     NUMNODES=2;
759     }
760     totalNECount=numElementsLocal*NE1;
761     faceNECount=0;
762    
763     /* do the internal elements first */
764     if (!periodic[1])
765     {
766     /* elements on boundary 010 (x1=0): */
767     #pragma omp parallel for private(i0,k,node0)
768     for (i0=0;i0<numElementsLocal-numNodesExternal;i0++) {
769     k=i0+faceNECount;
770     node0=i0;
771    
772     out->FaceElements->Id[k]=i0+totalNECount;
773     out->FaceElements->Tag[k] = 10;
774     out->FaceElements->Color[k] = 0; // i0%2;
775 bcumming 751 /* WARNING */
776     /* should be using numDOFLocal instead of numNodesLocal for the case of left hand domain with periodic[0]=true */
777 bcumming 730 if (useElementsOnFace) {
778     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
779     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
780     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal+1;
781     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal;
782     } else {
783     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
784     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
785     }
786     }
787     totalNECount += numElementsLocal-numNodesExternal;
788     faceNECount += numElementsLocal-numNodesExternal;
789    
790     /* ** elements on boundary 020 (x1=1): */
791    
792     #pragma omp parallel for private(i0,k,node0)
793     for (i0=0;i0<numElementsLocal-numNodesExternal;i0++) {
794     k=i0+faceNECount;
795     node0=i0+(NE1-1)*numNodesLocal;
796    
797     out->FaceElements->Id[k]=i0+totalNECount;
798     out->FaceElements->Tag[k] = 20;
799     out->FaceElements->Color[k] = 0; // i0%2+2;
800    
801     if (useElementsOnFace) {
802     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numNodesLocal+1;
803     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal;
804     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;
805     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
806     } else {
807     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numNodesLocal+1;
808     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal;
809     }
810     }
811     totalNECount += numElementsLocal-numNodesExternal;
812     faceNECount += numElementsLocal-numNodesExternal;
813     }
814     if ( !periodic[0])
815     {
816     /* ** elements on boundary 001 (x0=0): */
817     if( domLeft )
818     {
819     #pragma omp parallel for private(i1,k,node0)
820     for (i1=0;i1<NE1;i1++) {
821     k=i1+faceNECount;
822     node0=i1*numNodesLocal;
823    
824     out->FaceElements->Id[k]=i1+totalNECount;
825     out->FaceElements->Tag[k] = 1;
826     out->FaceElements->Color[k] = 0; //i1%2+4;
827    
828     if (useElementsOnFace) {
829     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numNodesLocal;
830     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
831     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
832     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal+1;
833     } else {
834     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numNodesLocal;
835     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
836    
837     }
838     }
839     totalNECount+=NE1;
840     faceNECount+=NE1;
841     }
842     /* ** elements on boundary 002 (x0=1): */
843     if( domRight )
844     {
845     #pragma omp parallel for private(i1,k,node0)
846     for (i1=0;i1<NE1;i1++) {
847     k=i1+faceNECount;
848     node0=(numNodesLocal-2)+i1*numNodesLocal;
849    
850     out->FaceElements->Id[k]=i1+totalNECount;
851     out->FaceElements->Tag[k] = 2;
852     out->FaceElements->Color[k] = 0; //i1%2+6;
853    
854     if (useElementsOnFace) {
855     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
856     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal+1;
857     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal;
858     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
859     } else {
860     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
861     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal+1;
862     }
863     }
864     totalNECount+=NE1;
865     faceNECount+=NE1;
866     }
867     }
868    
869     if( mpi_info->size>1 )
870 bcumming 751 {
871     if( !periodic[1] && (domInternal || domRight) )
872     {
873     // bottom left
874     k = faceNECount;
875     node0 = numNodesLocal*N1;
876 bcumming 730
877 bcumming 751 out->FaceElements->Id[k]=totalNECount;
878     out->FaceElements->Tag[k] = 10;
879     out->FaceElements->Color[k] = 0; //i1%2+6;
880 bcumming 730
881 bcumming 751 if (useElementsOnFace) {
882     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
883     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=0;
884     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=numNodesLocal;
885     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
886     } else {
887     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
888     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=0;
889     }
890     totalNECount++;
891     faceNECount++;
892 bcumming 730
893 bcumming 751 // top left
894     k = faceNECount;
895     node0 = (numNodesLocal+1)*N1 - 2;
896 bcumming 730
897 bcumming 751 out->FaceElements->Id[k]=totalNECount;
898     out->FaceElements->Tag[k] = 20;
899     out->FaceElements->Color[k] = 0; //i1%2+6;
900 bcumming 730
901 bcumming 751 if (useElementsOnFace) {
902     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=numNodesLocal*(N1-1);
903     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
904     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;
905     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=numNodesLocal*(N1-2);
906     } else {
907     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=numNodesLocal*(N1-1);
908     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
909     }
910     totalNECount++;
911     faceNECount++;
912     }
913     if( !periodic[1] && (domInternal || domLeft) )
914     {
915     // bottom right
916     k = faceNECount;
917     node0 = (numNodesLocal+nodesExternal[0])*N1;
918 bcumming 730
919 bcumming 751 out->FaceElements->Id[k]=totalNECount;
920     out->FaceElements->Tag[k] = 10;
921     out->FaceElements->Color[k] = 0; //i1%2+6;
922 bcumming 730
923 bcumming 751 if (useElementsOnFace) {
924     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=numNodesLocal-1;
925     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
926     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
927     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=2*numNodesLocal-1;
928     } else {
929     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=numNodesLocal-1;
930     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
931     }
932     totalNECount++;
933     faceNECount++;
934 bcumming 730
935 bcumming 751 // top right
936     k = faceNECount;
937     node0 = (numNodesLocal+1+nodesExternal[0])*N1-2;
938 bcumming 730
939 bcumming 751 out->FaceElements->Id[k]=totalNECount;
940     out->FaceElements->Tag[k] = 20;
941     out->FaceElements->Color[k] = 0; //i1%2+6;
942 bcumming 730
943 bcumming 751 if (useElementsOnFace) {
944     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
945     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=numNodesLocal*N1-1;
946     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=numNodesLocal*(N1-1)-1;
947     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
948     } else {
949     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
950     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=numNodesLocal*N1-1;
951     }
952     totalNECount++;
953     faceNECount++;
954     }
955     }
956 bcumming 730 out->FaceElements->minColor = 0;
957     out->FaceElements->maxColor = 0; //7;
958    
959 bcumming 751 out->FaceElements->elementDistribution->numInternal = 0;
960     out->FaceElements->elementDistribution->numBoundary = 0;
961 bcumming 730 if( domInternal )
962     {
963     if( !periodic[1] )
964     {
965     out->FaceElements->elementDistribution->numInternal = 2*(numElementsLocal-2);
966     out->FaceElements->elementDistribution->numBoundary = 4;
967     }
968     }
969     else if( mpi_info->size==1 )
970     {
971     if( !periodic[0] )
972     out->FaceElements->elementDistribution->numInternal += NE1*2;
973     if( !periodic[1] )
974     out->FaceElements->elementDistribution->numInternal += numElementsLocal*2;
975     }
976     else
977     {
978     if( !periodic[0] )
979     out->FaceElements->elementDistribution->numInternal += NE1;
980     if( !periodic[1] )
981     {
982     out->FaceElements->elementDistribution->numInternal += 2*(numElementsLocal-1-periodic[0]);
983     out->FaceElements->elementDistribution->numBoundary += 2*(1 + periodic[0]);
984     }
985     }
986    
987 bcumming 751 out->FaceElements->elementDistribution->numLocal = faceNECount;
988 bcumming 730
989     /* face elements done: */
990 bcumming 751
991     /* setup distribution info for other elements */
992     out->ContactElements->elementDistribution->numLocal = out->ContactElements->elementDistribution->numInternal = out->ContactElements->elementDistribution->numInternal = 0;
993     out->Points->elementDistribution->numLocal = out->Points->elementDistribution->numInternal = out->Points->elementDistribution->numInternal = 0;
994 bcumming 730
995     /* condense the nodes: */
996     Finley_Mesh_resolveNodeIds( out );
997    
998 bcumming 751 /* setup the CommBuffer */
999     Finley_NodeDistribution_formCommBuffer( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
1000     if ( !Finley_MPI_noError( mpi_info )) {
1001     if( Finley_noError() )
1002     Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );
1003     Paso_MPIInfo_dealloc( mpi_info );
1004     Finley_Mesh_dealloc(out);
1005     return NULL;
1006     }
1007    
1008     Finley_NodeDistribution_calculateIndexExternal( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
1009    
1010 bcumming 730 /* prepare mesh for further calculatuions:*/
1011 bcumming 751 Finley_Mesh_prepare(out) ;
1012 bcumming 730
1013     #ifdef Finley_TRACE
1014     printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
1015     #endif
1016    
1017 bcumming 751 if( !Finley_MPI_noError(mpi_info) )
1018     {
1019     if( Finley_noError() )
1020     Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );
1021     Paso_MPIInfo_dealloc( mpi_info );
1022     Finley_Mesh_dealloc(out);
1023     return NULL;
1024     }
1025    
1026     /* free up memory */
1027 bcumming 730 Paso_MPIInfo_dealloc( mpi_info );
1028    
1029     return out;
1030     }
1031     #endif
1032    
1033    
1034    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26