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

Contents of /trunk/finley/src/Mesh_rec4.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: 34191 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 a numElements[0] x numElements[1] mesh with first order elements (Rec4) in the rectangle */
18 /* [0,Length[0]] x [0,Length[1]]. order is the desired accuracy of the integration scheme. */
19
20
21 /**************************************************************/
22
23 /* Author: gross@access.edu.au */
24 /* Version: $Id$ */
25
26 /**************************************************************/
27
28 #include "RectangularMesh.h"
29
30 /**************************************************************/
31
32 #ifdef PASO_MPI
33 /* get the number of nodes/elements for domain with rank=rank, of size processors
34 where n is the total number of nodes/elements in the global domain */
35 static index_t domain_MODdim( index_t rank, index_t size, index_t n )
36 {
37 rank = size-rank-1;
38
39 if( rank < n%size )
40 return (index_t)floor(n/size)+1;
41 return (index_t)floor(n/size);
42 }
43
44
45 /* Determines the number of nodes/elements etc along an axis which is numElementsGlobal long for domain rank */
46 /* A bit messy, but it only has to be done once... */
47 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 *numDOFInternal )
48 {
49 index_t i0;
50 dim_t numNodesGlobal = numElementsGlobal+1;
51
52 (*numNodesLocal) = domain_MODdim( rank, size, numNodesGlobal );
53
54 numElementsLocal[0] = numNodesLocal[0]+1;
55 periodicLocal[0] = periodicLocal[1] = FALSE;
56 nodesExternal[0] = nodesExternal[1] = 1;
57 if( periodic )
58 {
59 if( size==1 )
60 {
61 numElementsLocal[0] = numElementsGlobal;
62 nodesExternal[0] = nodesExternal[1] = 0;
63 periodicLocal[0] = periodicLocal[1] = TRUE;
64 }
65 else
66 {
67 if( rank==0 )
68 {
69 periodicLocal[0] = TRUE;
70 numNodesLocal[0]++;
71 }
72 if( rank==(size-1) )
73 {
74 periodicLocal[1] = TRUE;
75 numNodesLocal[0]--;
76 numElementsLocal[0]--;
77 }
78 }
79 }
80 else if( !periodic )
81 {
82 if( rank==0 ){
83 nodesExternal[0]--;
84 numElementsLocal[0]--;
85 }
86 if( rank==(size-1) )
87 {
88 nodesExternal[1]--;
89 numElementsLocal[0]--;
90 }
91 }
92 numNodesExternal[0] = nodesExternal[0]+nodesExternal[1];
93
94 numElementsInternal[0] = numElementsLocal[0];
95 if( (rank==0) && (rank==size-1) );
96
97 else if( !periodic && ( (rank==0) ^ (rank==size-1) ) )
98 numElementsInternal[0] -= 1;
99 else
100 numElementsInternal[0] -= 2;
101
102 firstNode[0] = 0;
103 for( i0=0; i0<rank; i0++ )
104 firstNode[0] += domain_MODdim( i0, size, numNodesGlobal );
105
106 numDOFLocal[0] = numNodesLocal[0];
107 if( periodicLocal[0] )
108 numDOFLocal[0]--;
109
110 numDOFInternal[0] = numDOFLocal[0];
111 if( size>1 )
112 {
113 if( rank>0 || periodic ){
114 numDOFInternal[0] -= 1;
115 }
116 if( rank<(size-1) || periodic ){
117 numDOFInternal[0] -= 1;
118 }
119 }
120 DOFExternal[0] = nodesExternal[0];
121 DOFExternal[1] = nodesExternal[1];
122 }
123 #endif
124
125 Finley_Mesh* Finley_RectangularMesh_Rec4(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
126 #ifndef PASO_MPI
127 {
128 dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES,M0,M1;
129 index_t k,node0;
130 Finley_Mesh* out;
131 char name[50];
132 double time0=Finley_timer();
133 NE0=MAX(1,numElements[0]);
134 NE1=MAX(1,numElements[1]);
135 N0=NE0+1;
136 N1=NE1+1;
137
138 if (N0<=N1) {
139 M0=1;
140 M1=N0;
141 } else {
142 M0=N1;
143 M1=1;
144 }
145
146 NFaceElements=0;
147 if (!periodic[0]) {
148 NDOF0=N0;
149 NFaceElements+=2*NE1;
150 } else {
151 NDOF0=N0-1;
152 }
153 if (!periodic[1]) {
154 NDOF1=N1;
155 NFaceElements+=2*NE0;
156 } else {
157 NDOF1=N1-1;
158 }
159
160 /* allocate mesh: */
161
162 sprintf(name,"Rectangular %d x %d mesh",N0,N1);
163 out=Finley_Mesh_alloc(name,2,order);
164
165 if (! Finley_noError()) return NULL;
166
167 out->Elements=Finley_ElementFile_alloc(Rec4,out->order);
168 if (useElementsOnFace) {
169 out->FaceElements=Finley_ElementFile_alloc(Rec4Face,out->order);
170 out->ContactElements=Finley_ElementFile_alloc(Rec4Face_Contact,out->order);
171 } else {
172 out->FaceElements=Finley_ElementFile_alloc(Line2,out->order);
173 out->ContactElements=Finley_ElementFile_alloc(Line2_Contact,out->order);
174 }
175 out->Points=Finley_ElementFile_alloc(Point1,out->order);
176
177 if (! Finley_noError()) {
178 Finley_Mesh_dealloc(out);
179 return NULL;
180 }
181
182 /* allocate tables: */
183 Finley_NodeFile_allocTable(out->Nodes,N0*N1);
184 Finley_ElementFile_allocTable(out->Elements,NE0*NE1);
185 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
186
187 if (! Finley_noError()) {
188 Finley_Mesh_dealloc(out);
189 return NULL;
190 }
191 /* set nodes: */
192
193 #pragma omp parallel for private(i0,i1,k)
194 for (i1=0;i1<N1;i1++) {
195 for (i0=0;i0<N0;i0++) {
196 k=M0*i0+M1*i1;
197 out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(i0)/DBLE(N0-1)*Length[0];
198 out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
199 out->Nodes->Id[k]=i0+N0*i1;
200 out->Nodes->Tag[k]=0;
201 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1);
202 }
203 }
204 /* tags for the faces: */
205 if (!periodic[1]) {
206 for (i0=0;i0<N0;i0++) {
207 out->Nodes->Tag[M0*i0+M1*0]+=10;
208 out->Nodes->Tag[M0*i0+M1*(N1-1)]+=20;
209 }
210 }
211 if (!periodic[0]) {
212 for (i1=0;i1<N1;i1++) {
213 out->Nodes->Tag[M0*0+M1*i1]+=1;
214 out->Nodes->Tag[M0*(N0-1)+M1*i1]+=2;
215 }
216 }
217
218 /* set the elements: */
219
220 #pragma omp parallel for private(i0,i1,k,node0)
221 for (i1=0;i1<NE1;i1++) {
222 for (i0=0;i0<NE0;i0++) {
223 k=i0+NE0*i1;
224 node0=i0+i1*N0;
225
226 out->Elements->Id[k]=k;
227 out->Elements->Tag[k]=0;
228 out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);
229
230 out->Elements->Nodes[INDEX2(0,k,4)]=node0;
231 out->Elements->Nodes[INDEX2(1,k,4)]=node0+1;
232 out->Elements->Nodes[INDEX2(2,k,4)]=node0+N0+1;
233 out->Elements->Nodes[INDEX2(3,k,4)]=node0+N0;
234
235 }
236 }
237 out->Elements->minColor=0;
238 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
239
240 /* face elements: */
241
242 if (useElementsOnFace) {
243 NUMNODES=4;
244 } else {
245 NUMNODES=2;
246 }
247 totalNECount=NE0*NE1;
248 faceNECount=0;
249
250 /* elements on boundary 010 (x2=0): */
251 if (!periodic[1]) {
252 #pragma omp parallel for private(i0,k,node0)
253 for (i0=0;i0<NE0;i0++) {
254 k=i0+faceNECount;
255 node0=i0;
256
257 out->FaceElements->Id[k]=i0+totalNECount;
258 out->FaceElements->Tag[k]=10;
259 out->FaceElements->Color[k]=i0%2;
260
261 if (useElementsOnFace) {
262 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
263 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
264 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
265 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
266 } else {
267 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
268 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
269 }
270 }
271 totalNECount+=NE0;
272 faceNECount+=NE0;
273
274 /* ** elements on boundary 020 (x2=1): */
275
276 #pragma omp parallel for private(i0,k,node0)
277 for (i0=0;i0<NE0;i0++) {
278 k=i0+faceNECount;
279 node0=i0+(NE1-1)*N0;
280
281 out->FaceElements->Id[k]=i0+totalNECount;
282 out->FaceElements->Tag[k]=20;
283 out->FaceElements->Color[k]=i0%2+2;
284
285 if (useElementsOnFace) {
286 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0+1;
287 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
288 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;
289 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
290 } else {
291 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0+1;
292 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
293 }
294 }
295 totalNECount+=NE0;
296 faceNECount+=NE0;
297 }
298 if (!periodic[0]) {
299 /* ** elements on boundary 001 (x1=0): */
300 #pragma omp parallel for private(i1,k,node0)
301 for (i1=0;i1<NE1;i1++) {
302 k=i1+faceNECount;
303 node0=i1*N0;
304
305 out->FaceElements->Id[k]=i1+totalNECount;
306 out->FaceElements->Tag[k]=1;
307 out->FaceElements->Color[k]=i1%2+4;
308
309 if (useElementsOnFace) {
310 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
311 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
312 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
313 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
314 } else {
315 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
316 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
317
318 }
319 }
320 totalNECount+=NE1;
321 faceNECount+=NE1;
322
323 /* ** elements on boundary 002 (x1=1): */
324
325 #pragma omp parallel for private(i1,k,node0)
326 for (i1=0;i1<NE1;i1++) {
327 k=i1+faceNECount;
328 node0=(NE0-1)+i1*N0;
329
330 out->FaceElements->Id[k]=i1+totalNECount;
331 out->FaceElements->Tag[k]=2;
332 out->FaceElements->Color[k]=i1%2+6;
333
334 if (useElementsOnFace) {
335 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
336 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
337 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0;
338 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
339 } else {
340 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
341 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
342 }
343 }
344 totalNECount+=NE1;
345 faceNECount+=NE1;
346 }
347 out->FaceElements->minColor=0;
348 out->FaceElements->maxColor=7;
349
350
351 /* face elements done: */
352
353 /* condense the nodes: */
354
355 Finley_Mesh_resolveNodeIds(out);
356
357 /* prepare mesh for further calculatuions:*/
358
359 Finley_Mesh_prepare(out) ;
360
361 #ifdef Finley_TRACE
362 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
363 #endif
364
365 if (! Finley_noError()) {
366 Finley_Mesh_dealloc(out);
367 return NULL;
368 }
369 return out;
370 }
371
372 /*----------------------------------------------------------------------------
373 MPI VERSION
374
375
376 ASSUMPTIONS
377 ===========
378
379 - the domain dimension is large enough in the NE0 direction for each domain
380 to be 2 internal nodes wide in that direction. Face element calculation gets
381 buggered up otherwise.
382 - if dividing into two domains with periodic[0]=TRUE , then the global domain
383 must be at least 4 elements along the NE0 direction, otherwise redundant
384 nodes are generated.
385
386 These problems can be avoided by ensuring that numElements[0]/mpiSize >= 2
387
388 if domains are small enough to cause these problems, you probably don't
389 need to solve them in parallel. If you really do, rotate the domain so that the
390 long axis is aligned with NE0.
391 ------------------------------------------------------------------------------*/
392
393
394 #else
395 {
396 dim_t N0,N1,NE0,NE1,i0,i1, NE0_local, NDOF0,NDOF1, NFaceElements, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal, faceNECount, totalNECount,M0,M1,numDOFInternal;
397 index_t innerLoop=0;
398 index_t NUMNODES,k,firstNode=0, DOFcount=0, forwardDOF[2], backwardDOF[2], node0, node1;
399 index_t targetDomain=-1;
400 index_t *indexForward=NULL;
401 bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE;
402 Finley_Mesh* out;
403 char name[50];
404 double time0=Finley_timer();
405
406 NE0=MAX(1,numElements[0]);
407 NE1=MAX(1,numElements[1]);
408 N0=NE0+1;
409 N1=NE1+1;
410 Paso_MPIInfo *mpi_info = NULL;
411
412 /* get MPI information */
413 mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
414 if (! Finley_noError())
415 return NULL;
416
417 if( mpi_info->rank==0 )
418 domLeft = TRUE;
419 if( mpi_info->rank==mpi_info->size-1 )
420 domRight = TRUE;
421 if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
422 domInternal = TRUE;
423
424 /* dimensions of the local subdomain */
425 domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal, &numDOFInternal );
426
427 if (N0<=N1) {
428 M0=1;
429 M1=N0;
430 } else {
431 M0=N1;
432 M1=1;
433 }
434
435 NFaceElements=0;
436 if (!periodic[0])
437 {
438 if( domLeft )
439 NFaceElements+=NE1;
440 if( domRight )
441 NFaceElements+=NE1;
442 }
443 if (!periodic[1]) {
444 NDOF1=N1;
445 NFaceElements+=2*numElementsLocal;
446 } else {
447 NDOF1=N1-1;
448 }
449
450 /* allocate mesh: */
451
452 sprintf(name,"Rectangular %d x %d mesh",N0,N1);
453 out=Finley_Mesh_alloc( name, 2, order, mpi_info );
454
455 if (! Finley_noError()) return NULL;
456
457 out->Elements=Finley_ElementFile_alloc( Rec4, out->order, mpi_info );
458 if (useElementsOnFace) {
459 out->FaceElements=Finley_ElementFile_alloc(Rec4Face,out->order, mpi_info );
460 out->ContactElements=Finley_ElementFile_alloc(Rec4Face_Contact,out->order, mpi_info );
461 } else {
462 out->FaceElements=Finley_ElementFile_alloc(Line2,out->order, mpi_info );
463 out->ContactElements=Finley_ElementFile_alloc(Line2_Contact,out->order, mpi_info );
464 }
465 out->Points=Finley_ElementFile_alloc(Point1,out->order, mpi_info );
466 if (! Finley_noError()) {
467 Finley_Mesh_dealloc(out);
468 return NULL;
469 }
470
471 /* allocate tables: */
472 Finley_NodeFile_allocTable( out->Nodes, (numNodesLocal + numNodesExternal)*N1 );
473 Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal*NDOF1, (DOFExternal[0]+DOFExternal[1])*NDOF1, 0 );
474 Finley_ElementFile_allocTable(out->Elements,numElementsLocal*NE1);
475 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
476
477 if (! Finley_noError()) {
478 Finley_Mesh_dealloc(out);
479 return NULL;
480 }
481 /* set nodes: */
482
483 //#pragma omp parallel for private(i0,i1,k)
484 if( mpi_info->size==1 )
485 innerLoop = numNodesLocal;
486 else
487 innerLoop = numNodesLocal-periodicLocal[0];
488
489 for (k=0,i1=0;i1<N1;i1++) {
490 for ( i0=0;i0<innerLoop;i0++,k++) {
491 out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(i0+firstNode)/DBLE(N0-1)*Length[0];
492 out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
493 out->Nodes->Id[k]=i0+innerLoop*i1;
494 out->Nodes->Tag[k]=0;
495 out->Nodes->degreeOfFreedom[k]=i0%numDOFLocal + numDOFLocal*(i1%NDOF1);
496 }
497 }
498
499 /* tag 'em */
500 for( i0=0; i0<numNodesLocal; i0++ )
501 {
502 out->Nodes->Tag[i0] += 10; // bottom
503 out->Nodes->Tag[i0+numNodesLocal*(N1-1)] += 20; // top
504 }
505 if( domLeft && !periodicLocal[0] )
506 for( i1=0; i1<N1; i1++ )
507 out->Nodes->Tag[i1*numNodesLocal] += 1; // left
508
509 if( domRight && !periodicLocal[0] )
510 for( i1=0; i1<N1; i1++ )
511 out->Nodes->Tag[(i1+1)*numNodesLocal-1] += 2; // right
512
513 /* external nodes */
514 k = innerLoop*N1;
515 DOFcount = out->Nodes->degreeOfFreedom[k-1]+1;
516 if( mpi_info->size>1 )
517 {
518 /* add forward communication information for boundary nodes */
519 indexForward = TMPMEMALLOC( NDOF1, index_t );
520 if( domInternal || domRight || periodicLocal[0] )
521 {
522 for( int i=0; i<NDOF1; i++ )
523 indexForward[i] = i*numDOFLocal;
524 targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;
525 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );
526 }
527 if( domInternal || domLeft || periodicLocal[1] )
528 {
529 for( int i=0; i<NDOF1; i++ )
530 indexForward[i] = (i+1)*numDOFLocal - 1;
531 targetDomain = (mpi_info->rank+1) % mpi_info->size;
532 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );
533 }
534 TMPMEMFREE( indexForward );
535
536 /* LHS */
537 if( periodicLocal[0] )
538 {
539 for (i1=0; i1<N1; i1++, k++)
540 {
541 out->Nodes->Coordinates[INDEX2(0,k,2)]=Length[0];
542 out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
543 out->Nodes->Id[k]=k;
544 out->Nodes->Tag[k]=0;
545 out->Nodes->degreeOfFreedom[k] = (i1 % NDOF1)*numDOFLocal;
546 }
547 out->Nodes->Tag[k-N1] = 10;
548 out->Nodes->Tag[k-1] = 20;
549 for (i1=0; i1<N1; i1++, k++)
550 {
551 out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(NE0-1)/DBLE(NE0)*Length[0];
552 out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
553 out->Nodes->Id[k]=k;
554 out->Nodes->Tag[k]=0;
555 out->Nodes->degreeOfFreedom[k] = DOFcount + i1%NDOF1;
556 }
557 out->Nodes->Tag[k-N1] = 10;
558 out->Nodes->Tag[k-1] = 20;
559 DOFcount = out->Nodes->degreeOfFreedom[k-1]+1;
560
561 targetDomain = mpi_info->size-1;
562 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, out->Nodes->degreeOfFreedom+(numNodesLocal*N1) );
563 }
564 if( mpi_info->rank>0 )
565 {
566 for (i1=0; i1<N1; i1++, k++)
567 {
568 out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(firstNode-1)/DBLE(N0-1)*Length[0];
569 out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
570 out->Nodes->Id[k]=k;
571 out->Nodes->Tag[k]=0;
572 out->Nodes->degreeOfFreedom[k] = DOFcount + i1%NDOF1;
573 }
574 out->Nodes->Tag[k-N1] = 10;
575 out->Nodes->Tag[k-1] = 20;
576 DOFcount = out->Nodes->degreeOfFreedom[k-1]+1;
577
578 targetDomain = mpi_info->rank-1;
579 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, out->Nodes->degreeOfFreedom+(k-N1) );
580 }
581
582 /* RHS */
583 if( periodicLocal[1] || (mpi_info->rank < mpi_info->size-1) )
584 {
585 for (i1=0; i1<N1; i1++, k++)
586 {
587 out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(firstNode+numNodesLocal-periodicLocal[0])/DBLE(N0-1)*Length[0];
588 out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
589 out->Nodes->Id[k]=k;
590 out->Nodes->Tag[k]=0;
591 out->Nodes->degreeOfFreedom[k] = DOFcount + i1%NDOF1;
592 }
593 out->Nodes->Tag[k-N1] = 10;
594 out->Nodes->Tag[k-1] = 20;
595 DOFcount = out->Nodes->degreeOfFreedom[k-1]+1;
596
597 targetDomain = (mpi_info->rank+1) % mpi_info->size;
598 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, out->Nodes->degreeOfFreedom+(k-N1) );
599 }
600 }
601 out->Nodes->degreeOfFreedomDistribution->numInternal = NDOF1*numDOFInternal;
602 out->Nodes->degreeOfFreedomDistribution->numBoundary = out->Nodes->degreeOfFreedomDistribution->numLocal- out->Nodes->degreeOfFreedomDistribution->numInternal;
603 /* set the elements: */
604
605 /* internal */
606 //#pragma omp parallel for private(i0,i1,k,node0)
607 for ( i1=0, k=0; i1<NE1; i1++)
608 {
609 for ( i0=0; i0<numElementsInternal; i0++, k++)
610 {
611 node0=i0+i1*innerLoop;
612
613 out->Elements->Id[k]=k;
614 out->Elements->Tag[k]=0;
615 out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
616
617 out->Elements->Nodes[INDEX2(0,k,4)]=node0;
618 out->Elements->Nodes[INDEX2(1,k,4)]=node0+1;
619 out->Elements->Nodes[INDEX2(2,k,4)]=node0+innerLoop+1;
620 out->Elements->Nodes[INDEX2(3,k,4)]=node0+innerLoop;
621
622 }
623 }
624 /* boundary */
625 if( mpi_info->size>1 )
626 {
627 if( domInternal )
628 {
629 /* left */
630 for( i1=0; i1<NE1; i1++, k++ )
631 {
632 node1=numNodesLocal*N1 + i1;
633 node0=i1*numNodesLocal;
634
635 out->Elements->Id[k]=k;
636 out->Elements->Tag[k]=0;
637 out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
638
639 out->Elements->Nodes[INDEX2(0,k,4)]=node1;
640 out->Elements->Nodes[INDEX2(1,k,4)]=node0;
641 out->Elements->Nodes[INDEX2(2,k,4)]=node0+numNodesLocal;
642 out->Elements->Nodes[INDEX2(3,k,4)]=node1+1;
643 }
644
645 /* right */
646 for( i1=0; i1<NE1; i1++, k++ )
647 {
648 node0 = (i1+1)*numNodesLocal-1;
649 node1 = (numNodesLocal+1)*N1 + i1;
650
651 out->Elements->Id[k]=k;
652 out->Elements->Tag[k]=0;
653 out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
654
655 out->Elements->Nodes[INDEX2(0,k,4)]=node0;
656 out->Elements->Nodes[INDEX2(1,k,4)]=node1;
657 out->Elements->Nodes[INDEX2(2,k,4)]=node1+1;
658 out->Elements->Nodes[INDEX2(3,k,4)]=node0+numNodesLocal;
659 }
660 }
661 else if( domLeft )
662 {
663 /* left */
664 if( periodicLocal[0] )
665 {
666 node1 = numNodesLocal*N1;
667 node0 = innerLoop*N1;
668 for( i1=0; i1<NE1; i1++, k++, node1++, node0++ )
669 {
670 out->Elements->Id[k]=k;
671 out->Elements->Tag[k]=0;
672 out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
673
674 out->Elements->Nodes[INDEX2(0,k,4)]=node1;
675 out->Elements->Nodes[INDEX2(1,k,4)]=node0;
676 out->Elements->Nodes[INDEX2(2,k,4)]=node0+1;
677 out->Elements->Nodes[INDEX2(3,k,4)]=node1+1;
678 }
679 }
680 /* right */
681 for( i1=0; i1<NE1; i1++, k++ )
682 {
683 node0 = (i1+1)*innerLoop-1;
684 node1 = (numNodesLocal+periodicLocal[0])*N1 + i1;
685
686 out->Elements->Id[k]=k;
687 out->Elements->Tag[k]=0;
688 out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
689
690 out->Elements->Nodes[INDEX2(0,k,4)]=node0;
691 out->Elements->Nodes[INDEX2(1,k,4)]=node1;
692 out->Elements->Nodes[INDEX2(2,k,4)]=node1+1;
693 out->Elements->Nodes[INDEX2(3,k,4)]=node0+numNodesLocal-periodicLocal[0];
694 }
695 }
696 else if( domRight )
697 {
698 /* left */
699 for( i1=0; i1<NE1; i1++, k++ )
700 {
701 node1=numNodesLocal*N1 + i1;
702 node0=i1*numNodesLocal;
703
704 out->Elements->Id[k]=k;
705 out->Elements->Tag[k]=0;
706 out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
707
708 out->Elements->Nodes[INDEX2(0,k,4)]=node1;
709 out->Elements->Nodes[INDEX2(1,k,4)]=node0;
710 out->Elements->Nodes[INDEX2(2,k,4)]=node0+numNodesLocal;
711 out->Elements->Nodes[INDEX2(3,k,4)]=node1+1;
712 }
713 /* right */
714 if( periodicLocal[1] )
715 {
716 for( i1=0; i1<NE1; i1++, k++ )
717 {
718 node0=(i1+1)*numNodesLocal-1;
719 node1 = (numNodesLocal+1)*N1 + i1;
720
721 out->Elements->Id[k]=k;
722 out->Elements->Tag[k]=0;
723 out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
724
725 out->Elements->Nodes[INDEX2(0,k,4)]=node0;
726 out->Elements->Nodes[INDEX2(1,k,4)]=node1;
727 out->Elements->Nodes[INDEX2(2,k,4)]=node1+1;
728 out->Elements->Nodes[INDEX2(3,k,4)]=node0+numNodesLocal;
729 }
730 }
731 }
732 }
733 out->Elements->minColor=0;
734 out->Elements->maxColor=0; //COLOR_MOD(0)+3*COLOR_MOD(0);
735
736 if( domInternal )
737 {
738 out->Elements->elementDistribution->numInternal = NE1*(numElementsLocal-2);
739 out->Elements->elementDistribution->numBoundary = NE1*2;
740 }
741 else if( mpi_info->size==1 )
742 {
743 out->Elements->elementDistribution->numInternal = NE1*numElementsLocal;
744 }
745 else
746 {
747 out->Elements->elementDistribution->numInternal += NE1*(numElementsLocal-1-periodic[0]);
748 out->Elements->elementDistribution->numBoundary += NE1*(1 + periodic[0]);
749 }
750
751 out->FaceElements->elementDistribution->numLocal = out->FaceElements->elementDistribution->numInternal + out->FaceElements->elementDistribution->numBoundary;
752
753 /* face elements: */
754
755 if (useElementsOnFace) {
756 NUMNODES=4;
757 } else {
758 NUMNODES=2;
759 }
760 totalNECount=numElementsLocal*NE1;
761 faceNECount=0;
762
763 /* do the internal elements first */
764 if (!periodic[1])
765 {
766 /* elements on boundary 010 (x1=0): */
767 #pragma omp parallel for private(i0,k,node0)
768 for (i0=0;i0<numElementsLocal-numNodesExternal;i0++) {
769 k=i0+faceNECount;
770 node0=i0;
771
772 out->FaceElements->Id[k]=i0+totalNECount;
773 out->FaceElements->Tag[k] = 10;
774 out->FaceElements->Color[k] = 0; // i0%2;
775 /* WARNING */
776 /* should be using numDOFLocal instead of numNodesLocal for the case of left hand domain with periodic[0]=true */
777 if (useElementsOnFace) {
778 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
779 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
780 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal+1;
781 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal;
782 } else {
783 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
784 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
785 }
786 }
787 totalNECount += numElementsLocal-numNodesExternal;
788 faceNECount += numElementsLocal-numNodesExternal;
789
790 /* ** elements on boundary 020 (x1=1): */
791
792 #pragma omp parallel for private(i0,k,node0)
793 for (i0=0;i0<numElementsLocal-numNodesExternal;i0++) {
794 k=i0+faceNECount;
795 node0=i0+(NE1-1)*numNodesLocal;
796
797 out->FaceElements->Id[k]=i0+totalNECount;
798 out->FaceElements->Tag[k] = 20;
799 out->FaceElements->Color[k] = 0; // i0%2+2;
800
801 if (useElementsOnFace) {
802 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numNodesLocal+1;
803 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal;
804 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;
805 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
806 } else {
807 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numNodesLocal+1;
808 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal;
809 }
810 }
811 totalNECount += numElementsLocal-numNodesExternal;
812 faceNECount += numElementsLocal-numNodesExternal;
813 }
814 if ( !periodic[0])
815 {
816 /* ** elements on boundary 001 (x0=0): */
817 if( domLeft )
818 {
819 #pragma omp parallel for private(i1,k,node0)
820 for (i1=0;i1<NE1;i1++) {
821 k=i1+faceNECount;
822 node0=i1*numNodesLocal;
823
824 out->FaceElements->Id[k]=i1+totalNECount;
825 out->FaceElements->Tag[k] = 1;
826 out->FaceElements->Color[k] = 0; //i1%2+4;
827
828 if (useElementsOnFace) {
829 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numNodesLocal;
830 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
831 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
832 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal+1;
833 } else {
834 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numNodesLocal;
835 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
836
837 }
838 }
839 totalNECount+=NE1;
840 faceNECount+=NE1;
841 }
842 /* ** elements on boundary 002 (x0=1): */
843 if( domRight )
844 {
845 #pragma omp parallel for private(i1,k,node0)
846 for (i1=0;i1<NE1;i1++) {
847 k=i1+faceNECount;
848 node0=(numNodesLocal-2)+i1*numNodesLocal;
849
850 out->FaceElements->Id[k]=i1+totalNECount;
851 out->FaceElements->Tag[k] = 2;
852 out->FaceElements->Color[k] = 0; //i1%2+6;
853
854 if (useElementsOnFace) {
855 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
856 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal+1;
857 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal;
858 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
859 } else {
860 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
861 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal+1;
862 }
863 }
864 totalNECount+=NE1;
865 faceNECount+=NE1;
866 }
867 }
868
869 if( mpi_info->size>1 )
870 {
871 if( !periodic[1] && (domInternal || domRight) )
872 {
873 // bottom left
874 k = faceNECount;
875 node0 = numNodesLocal*N1;
876
877 out->FaceElements->Id[k]=totalNECount;
878 out->FaceElements->Tag[k] = 10;
879 out->FaceElements->Color[k] = 0; //i1%2+6;
880
881 if (useElementsOnFace) {
882 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
883 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=0;
884 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=numNodesLocal;
885 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
886 } else {
887 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
888 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=0;
889 }
890 totalNECount++;
891 faceNECount++;
892
893 // top left
894 k = faceNECount;
895 node0 = (numNodesLocal+1)*N1 - 2;
896
897 out->FaceElements->Id[k]=totalNECount;
898 out->FaceElements->Tag[k] = 20;
899 out->FaceElements->Color[k] = 0; //i1%2+6;
900
901 if (useElementsOnFace) {
902 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=numNodesLocal*(N1-1);
903 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
904 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;
905 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=numNodesLocal*(N1-2);
906 } else {
907 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=numNodesLocal*(N1-1);
908 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
909 }
910 totalNECount++;
911 faceNECount++;
912 }
913 if( !periodic[1] && (domInternal || domLeft) )
914 {
915 // bottom right
916 k = faceNECount;
917 node0 = (numNodesLocal+nodesExternal[0])*N1;
918
919 out->FaceElements->Id[k]=totalNECount;
920 out->FaceElements->Tag[k] = 10;
921 out->FaceElements->Color[k] = 0; //i1%2+6;
922
923 if (useElementsOnFace) {
924 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=numNodesLocal-1;
925 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
926 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
927 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=2*numNodesLocal-1;
928 } else {
929 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=numNodesLocal-1;
930 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
931 }
932 totalNECount++;
933 faceNECount++;
934
935 // top right
936 k = faceNECount;
937 node0 = (numNodesLocal+1+nodesExternal[0])*N1-2;
938
939 out->FaceElements->Id[k]=totalNECount;
940 out->FaceElements->Tag[k] = 20;
941 out->FaceElements->Color[k] = 0; //i1%2+6;
942
943 if (useElementsOnFace) {
944 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
945 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=numNodesLocal*N1-1;
946 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=numNodesLocal*(N1-1)-1;
947 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
948 } else {
949 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
950 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=numNodesLocal*N1-1;
951 }
952 totalNECount++;
953 faceNECount++;
954 }
955 }
956 out->FaceElements->minColor = 0;
957 out->FaceElements->maxColor = 0; //7;
958
959 out->FaceElements->elementDistribution->numInternal = 0;
960 out->FaceElements->elementDistribution->numBoundary = 0;
961 if( domInternal )
962 {
963 if( !periodic[1] )
964 {
965 out->FaceElements->elementDistribution->numInternal = 2*(numElementsLocal-2);
966 out->FaceElements->elementDistribution->numBoundary = 4;
967 }
968 }
969 else if( mpi_info->size==1 )
970 {
971 if( !periodic[0] )
972 out->FaceElements->elementDistribution->numInternal += NE1*2;
973 if( !periodic[1] )
974 out->FaceElements->elementDistribution->numInternal += numElementsLocal*2;
975 }
976 else
977 {
978 if( !periodic[0] )
979 out->FaceElements->elementDistribution->numInternal += NE1;
980 if( !periodic[1] )
981 {
982 out->FaceElements->elementDistribution->numInternal += 2*(numElementsLocal-1-periodic[0]);
983 out->FaceElements->elementDistribution->numBoundary += 2*(1 + periodic[0]);
984 }
985 }
986
987 out->FaceElements->elementDistribution->numLocal = faceNECount;
988
989 /* face elements done: */
990
991 /* setup distribution info for other elements */
992 out->ContactElements->elementDistribution->numLocal = out->ContactElements->elementDistribution->numInternal = out->ContactElements->elementDistribution->numInternal = 0;
993 out->Points->elementDistribution->numLocal = out->Points->elementDistribution->numInternal = out->Points->elementDistribution->numInternal = 0;
994
995 /* condense the nodes: */
996 Finley_Mesh_resolveNodeIds( out );
997
998 /* setup the CommBuffer */
999 Finley_NodeDistribution_formCommBuffer( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
1000 if ( !Finley_MPI_noError( mpi_info )) {
1001 if( Finley_noError() )
1002 Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );
1003 Paso_MPIInfo_dealloc( mpi_info );
1004 Finley_Mesh_dealloc(out);
1005 return NULL;
1006 }
1007
1008 Finley_NodeDistribution_calculateIndexExternal( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
1009
1010 /* prepare mesh for further calculatuions:*/
1011 Finley_Mesh_prepare(out) ;
1012
1013 #ifdef Finley_TRACE
1014 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
1015 #endif
1016
1017 if( !Finley_MPI_noError(mpi_info) )
1018 {
1019 if( Finley_noError() )
1020 Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );
1021 Paso_MPIInfo_dealloc( mpi_info );
1022 Finley_Mesh_dealloc(out);
1023 return NULL;
1024 }
1025
1026 /* free up memory */
1027 Paso_MPIInfo_dealloc( mpi_info );
1028
1029 return out;
1030 }
1031 #endif
1032
1033
1034

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26