/[escript]/branches/intelc_win32/finley/src/Mesh_line2.c
ViewVC logotype

Annotation of /branches/intelc_win32/finley/src/Mesh_line2.c

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26