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

Annotation of /trunk/finley/src/Mesh_line3.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: 27631 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 numElements[0] mesh with second order elements (Line3) in the interval */
18     /* [0,Length[0]]. order is the desired accuracy of the integration scheme. */
19    
20     /**************************************************************/
21    
22 jgs 150 /* Author: gross@access.edu.au */
23     /* Version: $Id$ */
24 jgs 82
25     /**************************************************************/
26    
27     #include "RectangularMesh.h"
28    
29    
30 bcumming 751
31     /**************************************************************
32     The code for Mesh_line[2/3].c is messy, it was used as a test
33     bed for getting the MPI distributed mesh figured out before
34     doing the rectangle/brick domains. The code is efficient, it just
35     uses too many if/else if/else statements. Unfortunately, calculating
36     a distributed mesh is an easy thing to visualise, but a devil to implement!
37    
38     Sorry if I don't get around
39     to neatening this code up, but rest assured, the more important
40     (and complex) rectangle and brick codes are much better organised.
41     **************************************************************/
42    
43     #ifdef PASO_MPI
44     /* get the number of nodes/elements for domain with rank=rank, of size processors
45     where n is the total number of nodes/elements in the global domain */
46     static index_t domain_MODdim( index_t rank, index_t size, index_t n )
47     {
48     rank = size-rank-1;
49    
50     if( rank < n%size )
51     return (index_t)floor(n/size)+1;
52     return (index_t)floor(n/size);
53     }
54    
55    
56     /* Determines the number of nodes/elements etc along an axis which is numElementsGlobal long for domain rank */
57     /* A bit messy, but it only has to be done once... */
58     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 )
59     {
60     index_t i0;
61     dim_t numNodesGlobal = numElementsGlobal+1;
62    
63     (*numNodesLocal) = 2*domain_MODdim( rank, size, numNodesGlobal ) - 1;
64     if( rank<size-1 ) // add on node for right hand boundary
65     (*numNodesLocal) += 1;
66    
67     numElementsLocal[0] = domain_MODdim( rank, size, numNodesGlobal )+1;
68     periodicLocal[0] = periodicLocal[1] = FALSE;
69     nodesExternal[0] = nodesExternal[1] = 1;
70     if( periodic )
71     {
72     if( size==1 )
73     {
74     numElementsLocal[0] = numElementsGlobal;
75     nodesExternal[0] = nodesExternal[1] = 0;
76     periodicLocal[0] = periodicLocal[1] = TRUE;
77     }
78     else
79     {
80     if( rank==0 )
81     {
82     periodicLocal[0] = TRUE;
83     numNodesLocal[0]++;
84     }
85     if( rank==(size-1) )
86     {
87     periodicLocal[1] = TRUE;
88     numNodesLocal[0]--;
89     numElementsLocal[0]--;
90     }
91     }
92     }
93     else if( !periodic )
94     {
95     if( rank==0 ){
96     nodesExternal[0]--;
97     numElementsLocal[0]--;
98     }
99     if( rank==(size-1) )
100     {
101     nodesExternal[1]--;
102     numElementsLocal[0]--;
103     }
104     }
105     nodesExternal[0]*=2;
106     numNodesExternal[0] = nodesExternal[0]+nodesExternal[1];
107    
108     numElementsInternal[0] = numElementsLocal[0];
109     if( (rank==0) && (rank==size-1) );
110    
111     else if( !periodic && ( (rank==0) ^ (rank==size-1) ) )
112     numElementsInternal[0] -= 1;
113     else
114     numElementsInternal[0] -= 2;
115    
116     firstNode[0] = 0;
117     for( i0=0; i0<rank; i0++ )
118     firstNode[0] += domain_MODdim( i0, size, numNodesGlobal );
119     firstNode[0] *= 2;
120    
121     numDOFLocal[0] = numNodesLocal[0];
122     if( periodicLocal[0] )
123     numDOFLocal[0]--;
124    
125     DOFExternal[0] = nodesExternal[0];
126     DOFExternal[1] = nodesExternal[1];
127    
128     if( size>1 )
129     {
130     DOFBoundary[0] = periodicLocal[0]*1 + (rank>0)*1;
131     DOFBoundary[1] = periodicLocal[1]*2 + (rank<(size-1))*2;
132     }
133     else
134     {
135     DOFBoundary[0] = DOFBoundary[1] = 0;
136     }
137     /* some debugging printf statements */
138     //printf( "rank/size = %d/%d\nNodes : %d Local, %d External[%d %d], First = %d\nElements : %d Local\nDOF : %d Local, External [%d %d], Boundary [%d %d]\nperiodicLocal [%d %d]\n\n", rank, size, *numNodesLocal, *numNodesExternal, nodesExternal[0], nodesExternal[1], *firstNode, *numElementsLocal, *numDOFLocal, DOFExternal[0], DOFExternal[1], DOFBoundary[0], DOFBoundary[1], periodicLocal[0], periodicLocal[1] );
139     }
140     #endif
141    
142    
143     Finley_Mesh* Finley_RectangularMesh_Line3(dim_t* numElements,double* Length,bool_t* periodic,index_t order,bool_t useElementsOnFace)
144     #ifndef PASO_MPI
145     {
146 jgs 123 dim_t N0,NE0,i0,NDOF0,NFaceElements,NUMNODES;
147     index_t k,node0;
148 jgs 82 Finley_Mesh* out;
149     char name[50];
150     double time0=Finley_timer();
151     NE0=MAX(1,numElements[0]);
152     N0=2*NE0+1;
153     if (!periodic[0]) {
154     NDOF0=N0;
155     NFaceElements=2;
156     } else {
157     NDOF0=N0-1;
158     NFaceElements=0;
159     }
160    
161     /* allocate mesh: */
162    
163     sprintf(name,"Rectangular mesh with %d nodes",N0);
164 bcumming 730 /* TEMPFIX */
165 bcumming 751
166 jgs 82 out=Finley_Mesh_alloc(name,1,order);
167 jgs 150 if (! Finley_noError()) return NULL;
168 jgs 82
169     out->Elements=Finley_ElementFile_alloc(Line3,out->order);
170     if (useElementsOnFace) {
171     out->FaceElements=Finley_ElementFile_alloc(Line3Face,out->order);
172     out->ContactElements=Finley_ElementFile_alloc(Line3Face_Contact,out->order);
173     } else {
174     out->FaceElements=Finley_ElementFile_alloc(Point1,out->order);
175     out->ContactElements=Finley_ElementFile_alloc(Point1_Contact,out->order);
176     }
177     out->Points=Finley_ElementFile_alloc(Point1,out->order);
178 jgs 150 if (! Finley_noError()) {
179 jgs 82 Finley_Mesh_dealloc(out);
180     return NULL;
181     }
182    
183     /* allocate tables: */
184 bcumming 751
185 jgs 82 Finley_NodeFile_allocTable(out->Nodes,N0);
186 bcumming 751
187 jgs 82 Finley_ElementFile_allocTable(out->Elements,NE0);
188     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
189 jgs 150 if (! Finley_noError()) {
190 jgs 82 Finley_Mesh_dealloc(out);
191     return NULL;
192     }
193    
194     /* set nodes: */
195    
196     #pragma omp parallel for private(i0,k)
197     for (i0=0;i0<N0;i0++) {
198     k=i0;
199     out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(i0)/DBLE(N0-1)*Length[0];
200     out->Nodes->Id[k]=k;
201     out->Nodes->Tag[k]=0;
202     out->Nodes->degreeOfFreedom[k]=(i0%NDOF0);
203     }
204     if (!periodic[0]) {
205     out->Nodes->Tag[0]=1;
206     out->Nodes->Tag[N0-1]=2;
207     }
208    
209     /* set the elements: */
210    
211     #pragma omp parallel for private(i0,k,node0)
212     for (i0=0;i0<NE0;i0++) {
213     k=i0;
214     node0=2*i0;
215    
216     out->Elements->Id[k]=k;
217     out->Elements->Tag[k]=0;
218     out->Elements->Color[k]=COLOR_MOD(i0);
219    
220     out->Elements->Nodes[INDEX2(0,k,3)]=node0;
221     out->Elements->Nodes[INDEX2(1,k,3)]=node0+2;
222     out->Elements->Nodes[INDEX2(2,k,3)]=node0+1;
223     }
224 jgs 123 out->Elements->minColor=0;
225     out->Elements->maxColor=COLOR_MOD(0);
226 jgs 82
227     /* face elements: */
228     if (useElementsOnFace) {
229     NUMNODES=3;
230     } else {
231     NUMNODES=1;
232     }
233    
234     if (!periodic[0]) {
235     out->FaceElements->Id[0]=NE0;
236     out->FaceElements->Tag[0]=1;
237     out->FaceElements->Color[0]=0;
238     if (useElementsOnFace) {
239     out->FaceElements->Nodes[INDEX2(0,0,NUMNODES)]=0;
240     out->FaceElements->Nodes[INDEX2(1,0,NUMNODES)]=2;
241     out->FaceElements->Nodes[INDEX2(2,0,NUMNODES)]=1;
242     } else {
243     out->FaceElements->Nodes[INDEX2(0,0,NUMNODES)]=0;
244     }
245    
246     out->FaceElements->Id[1]=NE0+1;
247     out->FaceElements->Tag[1]=2;
248     out->FaceElements->Color[1]=1;
249     if (useElementsOnFace) {
250     out->FaceElements->Nodes[INDEX2(0,1,NUMNODES)]=N0-1;
251     out->FaceElements->Nodes[INDEX2(1,1,NUMNODES)]=N0-3;
252     out->FaceElements->Nodes[INDEX2(2,1,NUMNODES)]=N0-2;
253     } else {
254     out->FaceElements->Nodes[INDEX2(0,1,NUMNODES)]=N0-1;
255     }
256     }
257 jgs 123 out->FaceElements->minColor=0;
258     out->FaceElements->maxColor=1;
259 jgs 82
260     /* face elements done: */
261    
262     /* condense the nodes: */
263    
264     Finley_Mesh_resolveNodeIds(out);
265    
266     /* prepare mesh for further calculations:*/
267    
268     Finley_Mesh_prepare(out) ;
269    
270 jgs 149 #ifdef Finley_TRACE
271 jgs 82 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
272 jgs 149 #endif
273 jgs 82
274 jgs 150 if (! Finley_noError()) {
275 jgs 82 Finley_Mesh_dealloc(out);
276     return NULL;
277     }
278     return out;
279     }
280 bcumming 751 #else
281     /* MPI version */
282     {
283     dim_t N0, NE0, NE0_local, i0, NDOF0, NFaceElements, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal, DOFBoundary[2];
284     index_t *numForward=NULL, *numBackward=NULL;
285     index_t NUMNODES, node0,k, i,firstNode=0, DOFcount=0, forwardDOF[4], backwardDOF[4];
286     index_t targetDomain=-1;
287     bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE;
288     Finley_Mesh* out=NULL;
289     char name[50];
290     double time0=Finley_timer();
291     Paso_MPIInfo *mpi_info = NULL;
292 jgs 82
293 bcumming 751 /* dimensions the global domain */
294     NE0=MAX(1,numElements[0]);
295     N0=2*NE0+1;
296     NDOF0=N0;
297     if( periodic[0] )
298     NDOF0--;
299    
300     /* get MPI information */
301     mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
302     if (! Finley_noError()) {
303     Finley_Mesh_dealloc(out);
304     return NULL;
305     }
306    
307     if( mpi_info->rank==0 )
308     domLeft = TRUE;
309     if( mpi_info->rank==mpi_info->size-1 )
310     domRight = TRUE;
311     if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
312     domInternal = TRUE;
313    
314     /* dimensions of the local subdomain */
315     domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal, DOFBoundary );
316    
317     /* form face element if the local domain contains a periodic boundary */
318     NFaceElements = 0;
319     if( !periodic[0] )
320     {
321     if(domLeft)
322     NFaceElements++;
323     if(domRight)
324     NFaceElements++;
325     }
326    
327     /* allocate mesh: */
328     sprintf(name,"Rectangular mesh with %d nodes",N0);
329     out=Finley_Mesh_alloc(name,1,order,mpi_info);
330     if (! Finley_noError()) return NULL;
331    
332     out->Elements=Finley_ElementFile_alloc(Line3,out->order,mpi_info);
333     if (useElementsOnFace) {
334     out->FaceElements=Finley_ElementFile_alloc(Line3Face,out->order,mpi_info);
335     out->ContactElements=Finley_ElementFile_alloc(Line3Face_Contact,out->order,mpi_info);
336     } else {
337     out->FaceElements=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
338     out->ContactElements=Finley_ElementFile_alloc(Point1_Contact,out->order,mpi_info);
339     }
340     out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
341     if (! Finley_noError()) {
342     Finley_Mesh_dealloc(out);
343     return NULL;
344     }
345    
346     /* allocate tables: */
347     Finley_NodeFile_allocTable(out->Nodes,numNodesLocal+numNodesExternal);
348     Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal, DOFExternal[0]+DOFExternal[1], 0 );
349     Finley_ElementFile_allocTable(out->Elements,numElementsLocal);
350     if( NFaceElements )
351     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
352     if (! Finley_noError()) {
353     Finley_Mesh_dealloc(out);
354     return NULL;
355     }
356    
357     /* set nodes: */
358     #pragma omp parallel for private(i0,k)
359     /* local nodes */
360     for (i0=0;i0<numNodesLocal;i0++) {
361     k=i0;
362     out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(i0+firstNode)/DBLE(N0-1)*Length[0];
363     out->Nodes->Id[k]=k;
364     out->Nodes->Tag[k]=0;
365     out->Nodes->degreeOfFreedom[k]=i0%(numDOFLocal);
366     }
367    
368     /* external nodes - do left then right hand side */
369     /* the following only applies if more than one domain */
370    
371     /* this is messy, could be done cleaner */
372     if( mpi_info->size>1 )
373     {
374     DOFcount = numNodesLocal;
375     k=numNodesLocal;
376     if( mpi_info->rank!=0 || periodicLocal[0] )
377     {
378     /* left hand boundary is periodic -
379     add the nodes/DOF that define the element on the right hand boundary */
380     if( periodicLocal[0] )
381     {
382     k--;
383     out->Nodes->Coordinates[INDEX2(0,k,1)]=Length[0];
384     out->Nodes->Id[k]=k;
385     out->Nodes->Tag[k]=0;
386     out->Nodes->degreeOfFreedom[k]=0;
387     DOFcount--;
388     k++;
389     out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(N0-2)/DBLE(N0-1)*Length[0];
390     out->Nodes->Id[k]=k;
391     out->Nodes->Tag[k]=0;
392     out->Nodes->degreeOfFreedom[k]=DOFcount++;
393     k++;
394     out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(N0-3)/DBLE(N0-1)*Length[0];
395     out->Nodes->Id[k]=k;
396     out->Nodes->Tag[k]=0;
397     out->Nodes->degreeOfFreedom[k]=DOFcount++;
398     k++;
399     }
400     /* left hand boundary with another subdomain, need to add the nodes/DOFs that
401     defines the element that spans the boundary */
402     else
403     {
404     out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(firstNode-1)/DBLE(N0-1)*Length[0];
405     out->Nodes->Id[k]=k;
406     out->Nodes->Tag[k]=0;
407     out->Nodes->degreeOfFreedom[k]=DOFcount++;
408     k++;
409     out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(firstNode-2)/DBLE(N0-1)*Length[0];
410     out->Nodes->Id[k]=k;
411     out->Nodes->Tag[k]=0;
412     out->Nodes->degreeOfFreedom[k]=DOFcount++;
413     k++;
414     }
415     }
416     if( mpi_info->rank!=(mpi_info->size-1) || periodicLocal[1] )
417     {
418     /* right hand boundary is periodic - add the external reference to the distribution */
419     if( periodicLocal[1] )
420     {
421     out->Nodes->Coordinates[INDEX2(0,k,1)]=Length[0];
422     out->Nodes->Id[k]=k;
423     out->Nodes->Tag[k]=0;
424     out->Nodes->degreeOfFreedom[k]=k;
425     k++;
426     }
427     /* right hand boundary with another subdomain, need to add the nodes/DOFs that
428     defines the element that spans the boundary */
429     else
430     {
431     out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(firstNode+numNodesLocal-periodicLocal[0])/DBLE(N0-1)*Length[0];
432     out->Nodes->Id[k]=k;
433     out->Nodes->Tag[k]=0;
434     out->Nodes->degreeOfFreedom[k]=DOFcount;
435     k++;
436     }
437     }
438     /* setup boundary DOF data */
439     if( domInternal )
440     {
441     targetDomain = mpi_info->rank-1;
442     forwardDOF[0] = 0;
443     backwardDOF[0] = numNodesLocal+1; backwardDOF[1] = numNodesLocal;
444     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, forwardDOF );
445     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, backwardDOF );
446    
447     targetDomain = mpi_info->rank+1;
448     forwardDOF[0] = numNodesLocal-2; forwardDOF[1] = numNodesLocal-1;
449     backwardDOF[0] = numNodesLocal+2;
450     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, forwardDOF );
451     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, backwardDOF );
452     }
453     else if( mpi_info->size>2 || (mpi_info->size==2&&!periodic[0]) )
454     if( domLeft )
455     {
456     if( periodicLocal[0] )
457     {
458     targetDomain = mpi_info->size-1;
459     forwardDOF[0] = 0;
460     backwardDOF[0] = numNodesLocal; backwardDOF[1] = numNodesLocal-1;
461     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, forwardDOF );
462     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, backwardDOF );
463     }
464     targetDomain = mpi_info->rank+1;
465     forwardDOF[0] = numNodesLocal-2-periodicLocal[0];
466     forwardDOF[1] = numNodesLocal-1-periodicLocal[0];
467     backwardDOF[0] = numNodesLocal + periodicLocal[0];
468     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, forwardDOF );
469     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, backwardDOF );
470     }
471     else
472     {
473     targetDomain = mpi_info->rank-1;
474     forwardDOF[0] = 0;
475     backwardDOF[0] = numNodesLocal+1;
476     backwardDOF[1] = numNodesLocal;
477     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, forwardDOF );
478     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, backwardDOF );
479    
480     if( periodicLocal[1] )
481     {
482     targetDomain = 0;
483     forwardDOF[0] = numNodesLocal-1; forwardDOF[1] = numNodesLocal-1;
484     backwardDOF[0] = numNodesLocal+1;
485     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, forwardDOF );
486     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, backwardDOF );
487     }
488     }
489     else if( mpi_info->size==2 && periodic[0])
490     if( domLeft )
491     {
492     targetDomain = mpi_info->size-1;
493     forwardDOF[0] = 0;
494     forwardDOF[1] = numDOFLocal-2;
495     forwardDOF[2] = numDOFLocal-1;
496     backwardDOF[0] = numDOFLocal+1;
497     backwardDOF[1] = numDOFLocal;
498     backwardDOF[2] = numDOFLocal+2;
499     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 3, forwardDOF );
500     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 3, backwardDOF );
501     }
502     else
503     {
504     targetDomain = 0;
505     forwardDOF[0] = numDOFLocal-2;
506     forwardDOF[1] = numDOFLocal-1;
507     forwardDOF[2] = 0;
508     backwardDOF[0] = numDOFLocal+2;
509     backwardDOF[1] = numDOFLocal+1;
510     backwardDOF[2] = numDOFLocal;
511     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 3, forwardDOF );
512     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 3, backwardDOF );
513     }
514    
515    
516     if (! Finley_MPI_noError(mpi_info)) {
517     Finley_Mesh_dealloc(out);
518     return NULL;
519     }
520     out->Nodes->degreeOfFreedomDistribution->numBoundary = DOFBoundary[0] + DOFBoundary[1];
521     out->Nodes->degreeOfFreedomDistribution->numInternal = numDOFLocal - out->Nodes->degreeOfFreedomDistribution->numBoundary;
522    
523     /*printf( "\n============NODES=============\n" );
524     for( k=0; k<numNodesLocal; k++ )
525     printf( "\tI\tId %d\tDOF %d\tcoord [%g]\n", out->Nodes->Id[k], out->Nodes->degreeOfFreedom[k] , out->Nodes->Coordinates[INDEX2(0,k,1)] );
526     for( k=numNodesLocal; k<numNodesLocal+numNodesExternal; k++ )
527     printf( "\tE\tId %d\tDOF %d\tcoord [%g]\n", out->Nodes->Id[k], out->Nodes->degreeOfFreedom[k] , out->Nodes->Coordinates[INDEX2(0,k,1)] );
528    
529     for( k=0; k<out->Nodes->degreeOfFreedomDistribution->numNeighbours; k++ )
530     {
531     if( out->Nodes->degreeOfFreedomDistribution->neighbours[k]>=0 )
532     {
533     printf( "\t%d boundary DOF { ", out->Nodes->degreeOfFreedomDistribution->edges[k]->numForward );
534     for( i0=0; i0<out->Nodes->degreeOfFreedomDistribution->edges[k]->numForward; i0++ )
535     printf( "%d ", out->Nodes->degreeOfFreedomDistribution->edges[k]->indexForward[i0] );
536     printf("} to %d\n", out->Nodes->degreeOfFreedomDistribution->neighbours[k] );
537     printf( "\t%d boundary DOF { ", out->Nodes->degreeOfFreedomDistribution->edges[k]->numBackward );
538     for( i0=0; i0<out->Nodes->degreeOfFreedomDistribution->edges[k]->numBackward; i0++ )
539     printf( "%d ", out->Nodes->degreeOfFreedomDistribution->edges[k]->indexBackward[i0] );
540     printf("} from %d\n", out->Nodes->degreeOfFreedomDistribution->neighbours[k] );
541     }
542     }*/
543     }
544    
545     /* set the elements: */
546     /* form internal elements */
547     //#pragma omp parallel for private(i0,k)
548     for (i0=0;i0<numElementsInternal;i0++)
549     {
550     k=i0;
551     node0=2*i0;
552     out->Elements->Id[k]=k;
553     out->Elements->Tag[k]=0;
554     out->Elements->Color[k]=0;
555    
556     out->Elements->Nodes[INDEX2(0,k,3)]=node0;
557     out->Elements->Nodes[INDEX2(1,k,3)]=node0+2;
558     out->Elements->Nodes[INDEX2(2,k,3)]=node0+1;
559     }
560    
561     /* followed by boundary elements... */
562     i0 = numElementsInternal;
563     if( mpi_info->size>1 )
564     {
565     /* left hand boundary */
566     if( mpi_info->rank>0 ) /* left hand boundary is an internal boundary */
567     {
568     k=i0;
569     node0=numNodesLocal;
570     out->Elements->Id[k]=k;
571     out->Elements->Tag[k]=0;
572     out->Elements->Color[k]=0;
573    
574     out->Elements->Nodes[INDEX2(0,k,3)]=node0+1;
575     out->Elements->Nodes[INDEX2(1,k,3)]=0;
576     out->Elements->Nodes[INDEX2(2,k,3)]=node0;
577     i0++;
578     }
579     else if( periodicLocal[0] ) /* left hand boundary is a periodic boundary */
580     {
581     k=i0;
582     node0 = numNodesLocal;
583     out->Elements->Id[k]=k;
584     out->Elements->Tag[k]=0;
585     out->Elements->Color[k]=0;
586    
587     out->Elements->Nodes[INDEX2(0,k,3)]=node0+1;
588     out->Elements->Nodes[INDEX2(1,k,3)]=node0-1;
589     out->Elements->Nodes[INDEX2(2,k,3)]=node0;
590     i0++;
591     }
592    
593     /* right hand boundary */
594     if( mpi_info->rank<mpi_info->size-1 ) /* right hand boundary is an internal boundary */
595     {
596     k=i0;
597     out->Elements->Id[k]=k;
598     out->Elements->Tag[k]=0;
599     out->Elements->Color[k]=0;
600    
601     out->Elements->Nodes[INDEX2(0,k,3)]=numNodesLocal-2-periodicLocal[0];
602     out->Elements->Nodes[INDEX2(1,k,3)]=nodesExternal[0]+numNodesLocal;
603     out->Elements->Nodes[INDEX2(2,k,3)]=numNodesLocal-1-periodicLocal[0];
604     }
605     else if( periodicLocal[1] ) /* right hand boundary is a periodic boundary */
606     {
607     /* no need to do anything */;
608     k=i0;
609     out->Elements->Id[k]=k;
610     out->Elements->Tag[k]=0;
611     out->Elements->Color[k]=0;
612    
613     out->Elements->Nodes[INDEX2(0,k,3)]=numNodesLocal-2;
614     out->Elements->Nodes[INDEX2(1,k,3)]=numNodesLocal+2;
615     out->Elements->Nodes[INDEX2(2,k,3)]=numNodesLocal-1;
616     }
617     }
618     out->Elements->minColor=0;
619     out->Elements->maxColor=0;
620    
621     out->Elements->elementDistribution->numLocal = numElementsLocal;
622     out->Elements->elementDistribution->numInternal = numElementsInternal;
623     out->Elements->elementDistribution->numBoundary = numElementsLocal - numElementsInternal;
624    
625     /* face elements: */
626     if (useElementsOnFace) {
627     NUMNODES=3;
628     } else {
629     NUMNODES=1;
630     }
631     if ( domLeft && !periodicLocal[0] )
632     {
633     out->FaceElements->Id[0]=i0-1;
634     out->FaceElements->Tag[0]=1;
635     out->FaceElements->Color[0]=0;
636     if (useElementsOnFace) {
637     out->FaceElements->Nodes[INDEX2(0,0,NUMNODES)]=0;
638     out->FaceElements->Nodes[INDEX2(1,0,NUMNODES)]=2;
639     out->FaceElements->Nodes[INDEX2(2,0,NUMNODES)]=1;
640     } else {
641     out->FaceElements->Nodes[INDEX2(0,0,NUMNODES)]=0;
642     }
643     i0++;
644     }
645     if( domRight && !periodicLocal[1])
646     {
647     out->FaceElements->Id[domLeft]=i0;
648     out->FaceElements->Tag[domLeft]=2;
649     out->FaceElements->Color[domLeft]=1;
650     if (useElementsOnFace) {
651     out->FaceElements->Nodes[INDEX2(0,domLeft,NUMNODES)]=numNodesLocal-3;
652     out->FaceElements->Nodes[INDEX2(1,domLeft,NUMNODES)]=numNodesLocal-1;
653     out->FaceElements->Nodes[INDEX2(2,domLeft,NUMNODES)]=numNodesLocal-2;
654     } else {
655     out->FaceElements->Nodes[INDEX2(0,domLeft,NUMNODES)]=numNodesLocal-1;
656     }
657     }
658     out->FaceElements->maxColor=0;
659     out->FaceElements->minColor=0;
660     out->FaceElements->elementDistribution->numBoundary = 0;
661     if( domLeft || domRight && !periodic[0] )
662     out->FaceElements->elementDistribution->numLocal = out->FaceElements->elementDistribution->numInternal = domLeft + domRight;
663     else
664     out->FaceElements->elementDistribution->numLocal = out->FaceElements->elementDistribution->numInternal = 0;
665    
666     /* setup distribution info for other elements */
667     out->ContactElements->elementDistribution->numLocal = out->ContactElements->elementDistribution->numInternal = out->ContactElements->elementDistribution->numInternal = 0;
668     out->Points->elementDistribution->numLocal = out->Points->elementDistribution->numInternal = out->Points->elementDistribution->numInternal = 0;
669    
670     /*printf( "\n============ELEMENTS (%d)=============\n", out->Elements->numElements );
671     for( k=0; k<out->Elements->elementDistribution->numInternal; k++ )
672     {
673     printf( "I\tId %d : nodes [%d %d %d]->DOF [%d %d %d]\n", out->Elements->Id[k], out->Elements->Nodes[INDEX2(0,k,3)], out->Elements->Nodes[INDEX2(1,k,3)], out->Elements->Nodes[INDEX2(2,k,3)], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(0,k,3)]], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(1,k,3)]], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(2,k,3)]] );
674     }
675     for( k=out->Elements->elementDistribution->numInternal; k<out->Elements->elementDistribution->numLocal; k++ )
676     {
677     printf( "B\tId %d : nodes [%d %d %d]->DOF [%d %d %d]\n", out->Elements->Id[k], out->Elements->Nodes[INDEX2(0,k,3)], out->Elements->Nodes[INDEX2(1,k,3)], out->Elements->Nodes[INDEX2(2,k,3)], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(0,k,3)]], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(1,k,3)]], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(2,k,3)]] );
678     }
679     printf( "\n" );
680     for( k=0; k<out->FaceElements->numElements; k++ )
681     {
682     if( NUMNODES==3 )
683     printf( "F\tId %d : nodes [%d %d %d]->DOF [%d %d %d]\n", out->FaceElements->Id[k], out->FaceElements->Nodes[INDEX2(0,k,3)], out->FaceElements->Nodes[INDEX2(1,k,3)], out->FaceElements->Nodes[INDEX2(2,k,3)], out->Nodes->degreeOfFreedom[out->FaceElements->Nodes[INDEX2(0,k,3)]], out->Nodes->degreeOfFreedom[out->FaceElements->Nodes[INDEX2(1,k,3)]], out->Nodes->degreeOfFreedom[out->FaceElements->Nodes[INDEX2(2,k,3)]] );
684     else
685     printf( "F\tId %d : nodes [%d]->DOF [%d]\n", out->FaceElements->Id[k], out->FaceElements->Nodes[INDEX2(0,k,1)], out->Nodes->degreeOfFreedom[out->FaceElements->Nodes[INDEX2(0,k,1)]] );
686     }*/
687     /* face elements done: */
688    
689     /* condense the nodes: */
690    
691     Finley_Mesh_resolveNodeIds(out);
692    
693     /* setup the CommBuffer */
694     Finley_NodeDistribution_formCommBuffer( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
695     if ( !Finley_MPI_noError( mpi_info )) {
696     if( Finley_noError() )
697     Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );
698     Paso_MPIInfo_dealloc( mpi_info );
699     Finley_Mesh_dealloc(out);
700     return NULL;
701     }
702    
703     Finley_NodeDistribution_calculateIndexExternal( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
704    
705     /* prepare mesh for further calculatuions:*/
706     Finley_Mesh_prepare(out);
707    
708     if( !Finley_MPI_noError(mpi_info) )
709     {
710     if( Finley_noError() )
711     Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );
712     Paso_MPIInfo_dealloc( mpi_info );
713     Finley_Mesh_dealloc(out);
714     return NULL;
715     }
716    
717     /* free up memory */
718     Paso_MPIInfo_dealloc( mpi_info );
719    
720     #ifdef Finley_TRACE
721     printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
722     #endif
723    
724     return out;
725     }
726     #endif
727    
728 jgs 82 /*
729 jgs 149 * Revision 1.3 2005/09/01 03:31:36 jgs
730     * Merge of development branch dev-02 back to main trunk on 2005-09-01
731     *
732 jgs 150 * Revision 1.2.2.2 2005/09/07 06:26:19 gross
733     * the solver from finley are put into the standalone package paso now
734     *
735 jgs 149 * Revision 1.2.2.1 2005/08/24 02:02:18 gross
736     * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.
737     *
738 jgs 123 * Revision 1.2 2005/07/08 04:07:53 jgs
739     * Merge of development branch back to main trunk on 2005-07-08
740 jgs 82 *
741 jgs 123 * Revision 1.1.1.1.2.1 2005/06/29 02:34:52 gross
742     * some changes towards 64 integers in finley
743     *
744     * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
745     * initial import of project esys2
746     *
747 jgs 82 * Revision 1.1.1.1 2004/06/24 04:00:40 johng
748     * Initial version of eys using boost-python.
749     *
750     *
751     */
752    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26