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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 750 - (show annotations)
Sun Jun 25 13:05:29 2006 UTC (16 years 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 /*
2 ************************************************************
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 */
12
13 /**************************************************************/
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 /* Author: gross@access.edu.au */
25 /* Version: $Id$ */
26
27 /**************************************************************/
28
29 #include "RectangularMesh.h"
30
31 /**************************************************************/
32
33
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 dim_t N0,NE0,i0,NDOF0,NFaceElements;
124 index_t NUMNODES,k;
125 Finley_Mesh* out;
126 char name[50];
127 double time0=Finley_timer();
128 NE0=MAX(1,numElements[0]);
129 N0=NE0+1;
130
131 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 if (! Finley_noError()) return NULL;
144
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 if (! Finley_noError()) {
155 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 if (! Finley_noError()) {
165 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 out->Elements->minColor=0;
197 out->Elements->maxColor=COLOR_MOD(0);
198
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 out->FaceElements->maxColor=1;
227 out->FaceElements->minColor=0;
228
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 #ifdef Finley_TRACE
241 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
242 #endif
243
244 if (! Finley_noError()) {
245 Finley_Mesh_dealloc(out);
246 return NULL;
247 }
248
249 return out;
250 }
251 #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
264 /* 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 /*
643 * 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 * 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 * 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 * Revision 1.2 2005/07/08 04:07:52 jgs
653 * Merge of development branch back to main trunk on 2005-07-08
654 *
655 * 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 * 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