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

Contents of /trunk/finley/src/Mesh_line3.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 751 - (show annotations)
Mon Jun 26 01:46:34 2006 UTC (13 years, 2 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 /*
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 second order elements (Line3) in the interval */
18 /* [0,Length[0]]. order is the desired accuracy of the integration scheme. */
19
20 /**************************************************************/
21
22 /* Author: gross@access.edu.au */
23 /* Version: $Id$ */
24
25 /**************************************************************/
26
27 #include "RectangularMesh.h"
28
29
30
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 dim_t N0,NE0,i0,NDOF0,NFaceElements,NUMNODES;
147 index_t k,node0;
148 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 /* TEMPFIX */
165
166 out=Finley_Mesh_alloc(name,1,order);
167 if (! Finley_noError()) return NULL;
168
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 if (! Finley_noError()) {
179 Finley_Mesh_dealloc(out);
180 return NULL;
181 }
182
183 /* allocate tables: */
184
185 Finley_NodeFile_allocTable(out->Nodes,N0);
186
187 Finley_ElementFile_allocTable(out->Elements,NE0);
188 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
189 if (! Finley_noError()) {
190 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 out->Elements->minColor=0;
225 out->Elements->maxColor=COLOR_MOD(0);
226
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 out->FaceElements->minColor=0;
258 out->FaceElements->maxColor=1;
259
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 #ifdef Finley_TRACE
271 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
272 #endif
273
274 if (! Finley_noError()) {
275 Finley_Mesh_dealloc(out);
276 return NULL;
277 }
278 return out;
279 }
280 #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
293 /* 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 /*
729 * 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 * 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 * 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 * Revision 1.2 2005/07/08 04:07:53 jgs
739 * Merge of development branch back to main trunk on 2005-07-08
740 *
741 * 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 * 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