/[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 749 - (show annotations)
Sun Jun 25 08:26:40 2006 UTC (16 years ago) by woo409
File MIME type: text/plain
File size: 23693 byte(s) 
+ Commiting some DEBUGGING code which will need to be removed later. I need to test on both win32 and altix. It appears something is being initialised incorrectly for Finley_ElementFile stuff on win32. I've modified Mesh_line2.c to assist in tracking this down (though it is common to anything that has useElementsOnFace true when constructed).
When I find it this is going to be one seriously stupid and expensive to find bug.
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 printf("\nDEBUG OUTPUT: %s(%d)\n: prediodic = %d", __FILE__,__LINE__,periodic[0]);
123 /* Serial/OpenMP version */
124 dim_t N0,NE0,i0,NDOF0,NFaceElements;
125 index_t NUMNODES,k;
126 Finley_Mesh* out;
127 char name[50];
128 double time0=Finley_timer();
129 NE0=MAX(1,numElements[0]);
130 N0=NE0+1;
131
132 if (!periodic[0]) {
133 NDOF0=N0;
134 NFaceElements=2;
135 } else {
136 NDOF0=N0-1;
137 NFaceElements=0;
138 }
139
140 /* allocate mesh: */
141
142 sprintf(name,"Rectangular mesh with %d nodes",N0);
143 out=Finley_Mesh_alloc(name,1,order);
144 if (! Finley_noError()) return NULL;
145
146 out->Elements=Finley_ElementFile_alloc(Line2,out->order);
147 if (useElementsOnFace) {
148 out->FaceElements=Finley_ElementFile_alloc(Line2Face,out->order);
149 out->ContactElements=Finley_ElementFile_alloc(Line2Face_Contact,out->order);
150 } else {
151 out->FaceElements=Finley_ElementFile_alloc(Point1,out->order);
152 out->ContactElements=Finley_ElementFile_alloc(Point1_Contact,out->order);
153 }
154 out->Points=Finley_ElementFile_alloc(Point1,out->order);
155 if (! Finley_noError()) {
156 Finley_Mesh_dealloc(out);
157 return NULL;
158 }
159
160 /* allocate tables: */
161
162 Finley_NodeFile_allocTable(out->Nodes,N0);
163 Finley_ElementFile_allocTable(out->Elements,NE0);
164 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
165 if (! Finley_noError()) {
166 Finley_Mesh_dealloc(out);
167 return NULL;
168 }
169
170 /* set nodes: */
171
172 #pragma omp parallel for private(i0,k)
173 for (i0=0;i0<N0;i0++) {
174 k=i0;
175 out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(i0)/DBLE(N0-1)*Length[0];
176 out->Nodes->Id[k]=k;
177 out->Nodes->Tag[k]=0;
178 out->Nodes->degreeOfFreedom[k]=(i0%NDOF0);
179 }
180 if (!periodic[0]) {
181 out->Nodes->Tag[0]=1;
182 out->Nodes->Tag[N0-1]=2;
183 }
184
185 /* set the elements: */
186
187 #pragma omp parallel for private(i0,k)
188 for (i0=0;i0<NE0;i0++) {
189 k=i0;
190 out->Elements->Id[k]=k;
191 out->Elements->Tag[k]=0;
192 out->Elements->Color[k]=COLOR_MOD(i0);
193
194 out->Elements->Nodes[INDEX2(0,k,2)]=i0;
195 out->Elements->Nodes[INDEX2(1,k,2)]=i0+1;
196 }
197 out->Elements->minColor=0;
198 out->Elements->maxColor=COLOR_MOD(0);
199
200 /* face elements: */
201 if (useElementsOnFace) {
202 NUMNODES=2;
203 } else {
204 NUMNODES=1;
205 }
206 if (!periodic[0]) {
207 out->FaceElements->Id[0]=NE0;
208 out->FaceElements->Tag[0]=1;
209 out->FaceElements->Color[0]=0;
210 if (useElementsOnFace) {
211 out->FaceElements->Nodes[INDEX2(0,0,NUMNODES)]=0;
212 out->FaceElements->Nodes[INDEX2(1,0,NUMNODES)]=1;
213 } else {
214 out->FaceElements->Nodes[INDEX2(0,0,NUMNODES)]=0;
215 }
216
217 out->FaceElements->Id[1]=NE0+1;
218 out->FaceElements->Tag[1]=2;
219 out->FaceElements->Color[1]=1;
220 if (useElementsOnFace) {
221 out->FaceElements->Nodes[INDEX2(0,1,NUMNODES)]=N0-1;
222 out->FaceElements->Nodes[INDEX2(1,1,NUMNODES)]=N0-2;
223 } else {
224 out->FaceElements->Nodes[INDEX2(0,1,NUMNODES)]=N0-1;
225 }
226 }
227 out->FaceElements->maxColor=1;
228 out->FaceElements->minColor=0;
229
230 /* face elements done: */
231
232
233 /* condense the nodes: */
234
235 Finley_Mesh_resolveNodeIds(out);
236
237 /* prepare mesh for further calculatuions:*/
238
239 Finley_Mesh_prepare(out) ;
240
241 #ifdef Finley_TRACE
242 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
243 #endif
244
245 if (! Finley_noError()) {
246 Finley_Mesh_dealloc(out);
247 return NULL;
248 }
249 // FIXME: dump out the mesh contents at this point to debug windows
250 printf("\nDEBUG OUTPUT: %s(%d)\nDumping mesh output:\n", __FILE__,__LINE__);
251 printf("Name: %s\n",out->Name);
252 printf("Reference Counter: %d\n",out->reference_counter);
253 printf("%s %d\n", out->FaceElements->ReferenceElement->Type->Name,out->FaceElements->numElements);
254 int NN=out->FaceElements->ReferenceElement->Type->numNodes;
255 for (int i=0;i<out->FaceElements->numElements;i++) {
256 printf("Id: %d Tag: %d Nodes:",out->FaceElements->Id[i],out->FaceElements->Tag[i]);
257 for (int j=0;j<NN;j++) printf(" %d",out->Nodes->Id[out->FaceElements->Nodes[INDEX2(j,i,NN)]]);
258 printf("\n");
259 }
260 //END FIXME
261
262 return out;
263 }
264 #else
265 {
266 /* Finley_Mesh* Finley_RectangularMesh_Line2(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace) { */
267 /* MPI version */
268 dim_t N0, NE0, NE0_local, i0, NDOF0, NFaceElements, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal;
269 index_t NUMNODES,k,firstNode=0, DOFcount=0, forwardDOF[2], backwardDOF[2];
270 index_t targetDomain=-1;
271 bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE;
272 Finley_Mesh* out=NULL;
273 char name[50];
274 double time0=Finley_timer();
275 Paso_MPIInfo *mpi_info = NULL;
276
277 /* dimensions the global domain */
278 NE0=MAX(1,numElements[0]);
279 N0=NE0+1;
280 NDOF0=N0;
281 if( periodic[0] )
282 NDOF0--;
283
284 /* get MPI information */
285 mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
286 if (! Finley_noError()) {
287 Finley_Mesh_dealloc(out);
288 return NULL;
289 }
290
291 if( mpi_info->rank==0 )
292 domLeft = TRUE;
293 if( mpi_info->rank==mpi_info->size-1 )
294 domRight = TRUE;
295 if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
296 domInternal = TRUE;
297
298 /* dimensions of the local subdomain */
299 domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal );
300
301 /* form face element if the local domain contains a periodic boundary */
302 NFaceElements = 0;
303 if( !periodic[0] )
304 {
305 if(domLeft)
306 NFaceElements++;
307 if(domRight)
308 NFaceElements++;
309 }
310
311
312 /* allocate mesh: */
313 sprintf(name,"Rectangular mesh with %d nodes",N0);
314 out=Finley_Mesh_alloc(name,1,order,mpi_info);
315 if (! Finley_noError()) return NULL;
316
317 out->Elements=Finley_ElementFile_alloc(Line2,out->order,mpi_info);
318 if (useElementsOnFace) {
319 out->FaceElements=Finley_ElementFile_alloc(Line2Face,out->order,mpi_info);
320 out->ContactElements=Finley_ElementFile_alloc(Line2Face_Contact,out->order,mpi_info);
321 } else {
322 out->FaceElements=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
323 out->ContactElements=Finley_ElementFile_alloc(Point1_Contact,out->order,mpi_info);
324 }
325 out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
326 if (! Finley_noError()) {
327 Finley_Mesh_dealloc(out);
328 return NULL;
329 }
330
331 /* allocate tables: */
332
333 Finley_NodeFile_allocTable(out->Nodes,numNodesLocal+numNodesExternal);
334 Finley_NodeDistibution_allocTable( out->Nodes->degreeOfFreedomDistribution, numNodesLocal, numNodesExternal, 0 );
335 Finley_ElementFile_allocTable(out->Elements,numElementsLocal);
336 if( NFaceElements )
337 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
338 if (! Finley_noError()) {
339 Finley_Mesh_dealloc(out);
340 return NULL;
341 }
342
343 /* set nodes: */
344 #pragma omp parallel for private(i0,k)
345 /* local nodes */
346 for (i0=0;i0<numNodesLocal;i0++) {
347 k=i0;
348 out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(i0+firstNode)/DBLE(N0-1)*Length[0];
349 out->Nodes->Id[k]=k;
350 out->Nodes->Tag[k]=0;
351 out->Nodes->degreeOfFreedom[k]=i0%(numDOFLocal);
352 }
353
354 /* external nodes - do left then right hand side */
355 /* the following only applies if more than one domain */
356 if( mpi_info->size>1 )
357 {
358 DOFcount = numNodesLocal;
359 k=numNodesLocal;
360 if( mpi_info->rank!=0 || periodicLocal[0] )
361 {
362 /* left hand boundary is periodic -
363 add the nodes/DOF that define the element on the right hand boundary */
364 if( periodicLocal[0] )
365 {
366 /* tag the left hand boundary appropriately */
367 out->Nodes->Tag[0]=1;
368 k--;
369 out->Nodes->Coordinates[INDEX2(0,k,1)]=Length[0];
370 out->Nodes->Id[k]=k;
371 out->Nodes->Tag[k]=2;
372 out->Nodes->degreeOfFreedom[k]=0;
373 DOFcount--;
374 k++;
375 out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(N0-2)/DBLE(N0-1)*Length[0];
376 out->Nodes->Id[k]=k;
377 out->Nodes->Tag[k]=0;
378 out->Nodes->degreeOfFreedom[k]=DOFcount++;
379 k++;
380 }
381 /* left hand boundary with another subdomain, need to add the node/DOF that
382 defines the element that spans the boundary */
383 else
384 {
385 out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(firstNode-1)/DBLE(N0-1)*Length[0];
386 out->Nodes->Id[k]=k;
387 out->Nodes->Tag[k]=0;
388 out->Nodes->degreeOfFreedom[k]=DOFcount++;
389 k++;
390 }
391 }
392 if( mpi_info->rank!=(mpi_info->size-1) || periodicLocal[1] )
393 {
394 /* right hand boundary is periodic - add the external reference to the distribution */
395 if( periodicLocal[1] )
396 {
397 out->Nodes->Coordinates[INDEX2(0,k,1)]=Length[0];
398 out->Nodes->Id[k]=k;
399 out->Nodes->Tag[k]=2;
400 out->Nodes->degreeOfFreedom[k]=k;
401 k++;
402 }
403 /* right hand boundary with another subdomain, need to add the node/DOF that
404 defines the element that spans the boundary */
405 else
406 {
407 out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(firstNode+numNodesLocal-periodicLocal[0])/DBLE(N0-1)*Length[0];
408 out->Nodes->Id[k]=k;
409 out->Nodes->Tag[k]=0;
410 out->Nodes->degreeOfFreedom[k]=DOFcount;
411 k++;
412 }
413
414 /* setup boundary DOF data */
415 if( domInternal )
416 {
417 targetDomain = mpi_info->rank-1;
418 forwardDOF[0] = 0;
419 backwardDOF[0] = numNodesLocal;
420 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, forwardDOF );
421 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, backwardDOF );
422
423 targetDomain = mpi_info->rank+1;
424 forwardDOF[0] = numNodesLocal-1;
425 backwardDOF[0] = numNodesLocal+1;
426 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, forwardDOF );
427 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, backwardDOF );
428 }
429 else if( domLeft )
430 {
431 if( periodicLocal[0] )
432 {
433 targetDomain = mpi_info->size-1;
434 forwardDOF[0] = 0;
435 backwardDOF[0] = numNodesLocal;
436 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, forwardDOF );
437 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, backwardDOF );
438 }
439 targetDomain = mpi_info->rank+1;
440 forwardDOF[0] = numNodesLocal-1-periodicLocal[0];
441 backwardDOF[0] = numNodesLocal + periodicLocal[0];
442 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, forwardDOF );
443 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, backwardDOF );
444 }
445 else
446 {
447 targetDomain = mpi_info->rank-1;
448 forwardDOF[0] = 0;
449 backwardDOF[0] = numNodesLocal;
450 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, forwardDOF );
451 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, backwardDOF );
452
453 if( periodicLocal[1] )
454 {
455 targetDomain = 0;
456 forwardDOF[0] = numNodesLocal-1;
457 backwardDOF[0] = numNodesLocal+1;
458 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, forwardDOF );
459 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, backwardDOF );
460 }
461 }
462
463 }
464 if (! Finley_MPI_noError(mpi_info)) {
465 Finley_Mesh_dealloc(out);
466 return NULL;
467 }
468
469 printf( "\n============NODES=============\n" );
470 for( k=0; k<numNodesLocal; k++ )
471 printf( "\tI\tId %d\tDOF %d\tcoord [%g]\n", out->Nodes->Id[k], out->Nodes->degreeOfFreedom[k] , out->Nodes->Coordinates[INDEX2(0,k,1)] );
472 for( k=numNodesLocal; k<numNodesLocal+numNodesExternal; k++ )
473 printf( "\tE\tId %d\tDOF %d\tcoord [%g]\n", out->Nodes->Id[k], out->Nodes->degreeOfFreedom[k] , out->Nodes->Coordinates[INDEX2(0,k,1)] );
474
475 for( k=0; k<out->Nodes->degreeOfFreedomDistribution->numNeighbours; k++ )
476 {
477 if( out->Nodes->degreeOfFreedomDistribution->neighbours[k]>=0 )
478 {
479 printf( "\t%d boundary DOF { ", out->Nodes->degreeOfFreedomDistribution->edges[k]->numForward );
480 for( i0=0; i0<out->Nodes->degreeOfFreedomDistribution->edges[k]->numForward; i0++ )
481 printf( "%d ", out->Nodes->degreeOfFreedomDistribution->edges[k]->indexForward[i0] );
482 printf("} to %d\n", out->Nodes->degreeOfFreedomDistribution->neighbours[k] );
483 printf( "\t%d boundary DOF { ", out->Nodes->degreeOfFreedomDistribution->edges[k]->numBackward );
484 for( i0=0; i0<out->Nodes->degreeOfFreedomDistribution->edges[k]->numBackward; i0++ )
485 printf( "%d ", out->Nodes->degreeOfFreedomDistribution->edges[k]->indexBackward[i0] );
486 printf("} from %d\n", out->Nodes->degreeOfFreedomDistribution->neighbours[k] );
487 }
488 }
489 }
490 else
491 {
492 printf( "\n============NODES=============\n" );
493 for( k=0; k<numNodesLocal; k++ )
494 printf( "\tI\tId %d\tDOF %d\tcoord [%g]\n", out->Nodes->Id[k], out->Nodes->degreeOfFreedom[k] , out->Nodes->Coordinates[INDEX2(0,k,1)] );
495 }
496
497 /* set the elements: */
498 /* form internal elements */
499 //#pragma omp parallel for private(i0,k)
500 for (i0=0;i0<numElementsInternal;i0++)
501 {
502 k=i0;
503 out->Elements->Id[k]=k;
504 out->Elements->Tag[k]=0;
505 out->Elements->Color[k]=0;
506
507 out->Elements->Nodes[INDEX2(0,k,2)]=i0;
508 out->Elements->Nodes[INDEX2(1,k,2)]=i0+1;
509 }
510
511 /* followed by boundary elements... */
512 i0 = numElementsInternal;
513 if( mpi_info->size>1 )
514 {
515 /* left hand boundary */
516 if( mpi_info->rank>0 ) /* left hand boundary is an internal boundary */
517 {
518 k=i0;
519 out->Elements->Id[k]=k;
520 out->Elements->Tag[k]=0;
521 out->Elements->Color[k]=0;
522
523 out->Elements->Nodes[INDEX2(0,k,2)]=i0+1;
524 out->Elements->Nodes[INDEX2(1,k,2)]=0;
525 i0++;
526 }
527 else if( periodicLocal[0] ) /* left hand boundary is a periodic 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)]=i0+2;
535 out->Elements->Nodes[INDEX2(1,k,2)]=i0+1;
536 i0++;
537 }
538
539 /* right hand boundary */
540 if( mpi_info->rank<mpi_info->size-1 ) /* right hand boundary is an internal boundary */
541 {
542 k=i0;
543 out->Elements->Id[k]=k;
544 out->Elements->Tag[k]=0;
545 out->Elements->Color[k]=0;
546
547 out->Elements->Nodes[INDEX2(0,k,2)]=numNodesLocal-1-periodicLocal[0];
548 out->Elements->Nodes[INDEX2(1,k,2)]=nodesExternal[0]+numNodesLocal;
549 i0++;
550 }
551 else if( periodicLocal[1] ) /* right hand boundary is a periodic boundary */
552 {
553 /* no need to do anything */;
554 k=i0;
555 out->Elements->Id[k]=k;
556 out->Elements->Tag[k]=0;
557 out->Elements->Color[k]=0;
558
559 out->Elements->Nodes[INDEX2(0,k,2)]=numNodesLocal-1;
560 out->Elements->Nodes[INDEX2(1,k,2)]=numNodesLocal+1;
561 }
562 }
563 out->Elements->minColor=0;
564 out->Elements->maxColor=0;
565
566 out->Elements->elementDistribution->numLocal = numElementsLocal;
567 out->Elements->elementDistribution->numInternal = numElementsInternal;
568 out->Elements->elementDistribution->numBoundary = numElementsLocal - numElementsInternal;
569
570 /* face elements: */
571 if (useElementsOnFace) {
572 NUMNODES=2;
573 } else {
574 NUMNODES=1;
575 }
576 if ( domLeft && !periodicLocal[0] )
577 {
578 out->FaceElements->Id[0]=i0-1;
579 out->FaceElements->Tag[0]=1;
580 out->FaceElements->Color[0]=0;
581 if (useElementsOnFace) {
582 out->FaceElements->Nodes[INDEX2(0,0,NUMNODES)]=0;
583 out->FaceElements->Nodes[INDEX2(1,0,NUMNODES)]=1;
584 } else {
585 out->FaceElements->Nodes[INDEX2(0,0,NUMNODES)]=0;
586 }
587 i0++;
588 }
589 if( domRight && !periodicLocal[1])
590 {
591 out->FaceElements->Id[domLeft]=i0;
592 out->FaceElements->Tag[domLeft]=2;
593 out->FaceElements->Color[domLeft]=1;
594 /* TODO */
595 /* check that the correct indices have been used */
596 if (useElementsOnFace) {
597 out->FaceElements->Nodes[INDEX2(0,domLeft,NUMNODES)]=numNodesLocal-2;
598 out->FaceElements->Nodes[INDEX2(1,domLeft,NUMNODES)]=numNodesLocal-1;
599 } else {
600 out->FaceElements->Nodes[INDEX2(0,domLeft,NUMNODES)]=numNodesLocal-1;
601 }
602 }
603 out->FaceElements->maxColor=0;
604 out->FaceElements->minColor=0;
605 if( domLeft || domRight && !periodic[0] )
606 out->FaceElements->elementDistribution->numLocal = out->FaceElements->elementDistribution->numInternal = domLeft + domRight;
607
608 printf( "\n============ELEMENTS (%d)=============\n", out->Elements->numElements );
609 for( k=0; k<out->Elements->elementDistribution->numInternal; k++ )
610 {
611 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)]] );
612 }
613 for( k=out->Elements->elementDistribution->numInternal; k<out->Elements->elementDistribution->numLocal; k++ )
614 {
615 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)]] );
616 }
617 for( k=0; k<out->FaceElements->numElements; k++ )
618 {
619 if( NUMNODES==2 )
620 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)]] );
621 else
622 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)]] );
623 }
624 /* face elements done: */
625
626 /* condense the nodes: */
627 Finley_Mesh_resolveNodeIds(out);
628
629 /* prepare mesh for further calculatuions:*/
630
631 /* TEMPFIX */
632 /* this isn't needed for a mesh generated this way */
633 //Finley_Mesh_prepare(out);
634
635 if (! Finley_MPI_noError( mpi_info )) {
636 Paso_MPIInfo_dealloc( mpi_info );
637 Finley_Mesh_dealloc(out);
638 return NULL;
639 }
640
641 /* free up memory */
642 Paso_MPIInfo_dealloc( mpi_info );
643
644 #ifdef Finley_TRACE
645 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
646 #endif
647
648 return out;
649 }
650 #endif
651
652
653
654
655 /*
656 * Revision 1.3 2005/09/01 03:31:36 jgs
657 * Merge of development branch dev-02 back to main trunk on 2005-09-01
658 *
659 * Revision 1.2.2.2 2005/09/07 06:26:19 gross
660 * the solver from finley are put into the standalone package paso now
661 *
662 * Revision 1.2.2.1 2005/08/24 02:02:18 gross
663 * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.
664 *
665 * Revision 1.2 2005/07/08 04:07:52 jgs
666 * Merge of development branch back to main trunk on 2005-07-08
667 *
668 * Revision 1.1.1.1.2.2 2005/06/30 01:53:56 gross
669 * a bug in coloring fixed
670 *
671 * Revision 1.1.1.1.2.1 2005/06/29 02:34:52 gross
672 * some changes towards 64 integers in finley
673 *
674 * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
675 * initial import of project esys2
676 *
677 * Revision 1.1.1.1 2004/06/24 04:00:40 johng
678 * Initial version of eys using boost-python.
679 *
680 *
681 */
682

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
  ViewVC Help
Powered by ViewVC 1.1.26