/[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 782 - (show annotations)
Tue Jul 18 00:47:47 2006 UTC (12 years, 9 months ago) by bcumming
File MIME type: text/plain
File size: 26100 byte(s)
Large number of changes to Finley for meshing in MPI.

- optimisation and neatening up of rectcanglular mesh generation code
- first and second order 1D, 2D and 3D rectangular meshes are now
  available in finley and escript using MPI.
- reduced meshes now generated in MPI, and interpolation to and from 
  reduced data types now supported.  

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 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 )
47 {
48 index_t i0;
49 dim_t numNodesGlobal = numElementsGlobal+1;
50
51 (*numNodesLocal) = domain_MODdim( rank, size, numNodesGlobal );
52
53 numElementsLocal[0] = numNodesLocal[0]+1;
54 periodicLocal[0] = periodicLocal[1] = FALSE;
55 nodesExternal[0] = nodesExternal[1] = 1;
56 if( periodic )
57 {
58 if( size==1 )
59 {
60 numElementsLocal[0] = numElementsGlobal;
61 nodesExternal[0] = nodesExternal[1] = 0;
62 periodicLocal[0] = periodicLocal[1] = TRUE;
63 }
64 else
65 {
66 if( rank==0 )
67 {
68 periodicLocal[0] = TRUE;
69 numNodesLocal[0]++;
70 }
71 if( rank==(size-1) )
72 {
73 periodicLocal[1] = TRUE;
74 numNodesLocal[0]--;
75 numElementsLocal[0]--;
76 }
77 }
78 }
79 else if( !periodic )
80 {
81 if( rank==0 ){
82 nodesExternal[0]--;
83 numElementsLocal[0]--;
84 }
85 if( rank==(size-1) )
86 {
87 nodesExternal[1]--;
88 numElementsLocal[0]--;
89 }
90 }
91 numNodesExternal[0] = nodesExternal[0]+nodesExternal[1];
92
93 numElementsInternal[0] = numElementsLocal[0];
94 if( (rank==0) && (rank==size-1) );
95
96 else if( !periodic && ( (rank==0) ^ (rank==size-1) ) )
97 numElementsInternal[0] -= 1;
98 else
99 numElementsInternal[0] -= 2;
100
101 firstNode[0] = 0;
102 for( i0=0; i0<rank; i0++ )
103 firstNode[0] += domain_MODdim( i0, size, numNodesGlobal );
104
105 numDOFLocal[0] = numNodesLocal[0];
106 if( periodicLocal[0] )
107 numDOFLocal[0]--;
108
109 numDOFInternal[0] = numDOFLocal[0];
110 if( size>1 )
111 {
112 if( rank>0 || periodic ){
113 numDOFInternal[0] -= 1;
114 }
115 if( rank<(size-1) || periodic ){
116 numDOFInternal[0] -= 1;
117 }
118 }
119 DOFExternal[0] = nodesExternal[0];
120 DOFExternal[1] = nodesExternal[1];
121 }
122 #endif
123
124 #ifdef PASO_MPI
125 Finley_Mesh* Finley_RectangularMesh_Rec4_singleCPU(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace,Paso_MPIInfo *mpi_info)
126 #else
127 Finley_Mesh* Finley_RectangularMesh_Rec4(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
128 #endif
129 {
130 dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES,M0,M1;
131 index_t k,node0;
132 Finley_Mesh* out;
133 char name[50];
134 double time0=Finley_timer();
135 NE0=MAX(1,numElements[0]);
136 NE1=MAX(1,numElements[1]);
137 N0=NE0+1;
138 N1=NE1+1;
139
140 if (N0<=N1) {
141 M0=1;
142 M1=N0;
143 } else {
144 M0=N1;
145 M1=1;
146 }
147
148 NFaceElements=0;
149 if (!periodic[0]) {
150 NDOF0=N0;
151 NFaceElements+=2*NE1;
152 } else {
153 NDOF0=N0-1;
154 }
155 if (!periodic[1]) {
156 NDOF1=N1;
157 NFaceElements+=2*NE0;
158 } else {
159 NDOF1=N1-1;
160 }
161
162 /* allocate mesh: */
163
164 sprintf(name,"Rectangular %d x %d mesh",N0,N1);
165 #ifdef PASO_MPI
166 out=Finley_Mesh_alloc(name,2,order,mpi_info);
167
168 if (! Finley_noError()) return NULL;
169
170 out->Elements=Finley_ElementFile_alloc(Rec4,out->order,mpi_info);
171 if (useElementsOnFace) {
172 out->FaceElements=Finley_ElementFile_alloc(Rec4Face,out->order,mpi_info);
173 out->ContactElements=Finley_ElementFile_alloc(Rec4Face_Contact,out->order,mpi_info);
174 } else {
175 out->FaceElements=Finley_ElementFile_alloc(Line2,out->order,mpi_info);
176 out->ContactElements=Finley_ElementFile_alloc(Line2_Contact,out->order,mpi_info);
177 }
178 out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
179
180 if (! Finley_noError()) {
181 Finley_Mesh_dealloc(out);
182 return NULL;
183 }
184
185 /* allocate tables: */
186 Finley_NodeFile_allocTable(out->Nodes,N0*N1);
187 Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, N0*N1, 0, 0);
188 Finley_ElementFile_allocTable(out->Elements,NE0*NE1);
189 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
190 #else
191 out=Finley_Mesh_alloc(name,2,order);
192
193 if (! Finley_noError()) return NULL;
194
195 out->Elements=Finley_ElementFile_alloc(Rec4,out->order);
196 if (useElementsOnFace) {
197 out->FaceElements=Finley_ElementFile_alloc(Rec4Face,out->order);
198 out->ContactElements=Finley_ElementFile_alloc(Rec4Face_Contact,out->order);
199 } else {
200 out->FaceElements=Finley_ElementFile_alloc(Line2,out->order);
201 out->ContactElements=Finley_ElementFile_alloc(Line2_Contact,out->order);
202 }
203 out->Points=Finley_ElementFile_alloc(Point1,out->order);
204
205 if (! Finley_noError()) {
206 Finley_Mesh_dealloc(out);
207 return NULL;
208 }
209
210 /* allocate tables: */
211 Finley_NodeFile_allocTable(out->Nodes,N0*N1);
212 Finley_ElementFile_allocTable(out->Elements,NE0*NE1);
213 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
214 #endif
215 if (! Finley_noError()) {
216 Finley_Mesh_dealloc(out);
217 return NULL;
218 }
219 /* set nodes: */
220 #pragma omp parallel for private(i0,i1,k)
221 for (i1=0;i1<N1;i1++) {
222 for (i0=0;i0<N0;i0++) {
223 k=M0*i0+M1*i1;
224 out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(i0)/DBLE(N0-1)*Length[0];
225 out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
226 out->Nodes->Id[k]=i0+N0*i1;
227 out->Nodes->Tag[k]=0;
228 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1);
229 #ifdef PASO_MPI
230 out->Nodes->Dom[k]=NODE_INTERNAL;
231 #endif
232 }
233 }
234 /* tags for the faces: */
235 if (!periodic[1]) {
236 for (i0=0;i0<N0;i0++) {
237 out->Nodes->Tag[M0*i0+M1*0]+=10;
238 out->Nodes->Tag[M0*i0+M1*(N1-1)]+=20;
239 }
240 }
241 if (!periodic[0]) {
242 for (i1=0;i1<N1;i1++) {
243 out->Nodes->Tag[M0*0+M1*i1]+=1;
244 out->Nodes->Tag[M0*(N0-1)+M1*i1]+=2;
245 }
246 }
247
248 /* set the elements: */
249
250 #pragma omp parallel for private(i0,i1,k,node0)
251 for (i1=0;i1<NE1;i1++) {
252 for (i0=0;i0<NE0;i0++) {
253 k=i0+NE0*i1;
254 node0=i0+i1*N0;
255
256 out->Elements->Id[k]=k;
257 out->Elements->Tag[k]=0;
258 out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);
259 #ifdef PASO_MPI
260 out->Elements->Dom[k]=ELEMENT_INTERNAL;
261 #endif
262
263 out->Elements->Nodes[INDEX2(0,k,4)]=node0;
264 out->Elements->Nodes[INDEX2(1,k,4)]=node0+1;
265 out->Elements->Nodes[INDEX2(2,k,4)]=node0+N0+1;
266 out->Elements->Nodes[INDEX2(3,k,4)]=node0+N0;
267
268 }
269 }
270 out->Elements->minColor=0;
271 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
272
273 /* face elements: */
274
275 if (useElementsOnFace) {
276 NUMNODES=4;
277 } else {
278 NUMNODES=2;
279 }
280 totalNECount=NE0*NE1;
281 faceNECount=0;
282
283 /* elements on boundary 010 (x2=0): */
284 if (!periodic[1]) {
285 #pragma omp parallel for private(i0,k,node0)
286 for (i0=0;i0<NE0;i0++) {
287 k=i0+faceNECount;
288 node0=i0;
289
290 out->FaceElements->Id[k]=i0+totalNECount;
291 out->FaceElements->Tag[k]=10;
292 out->FaceElements->Color[k]=i0%2;
293 #ifdef PASO_MPI
294 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
295 #endif
296
297 if (useElementsOnFace) {
298 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
299 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
300 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
301 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
302 } else {
303 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
304 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
305 }
306 }
307 totalNECount+=NE0;
308 faceNECount+=NE0;
309
310 /* ** elements on boundary 020 (x2=1): */
311
312 #pragma omp parallel for private(i0,k,node0)
313 for (i0=0;i0<NE0;i0++) {
314 k=i0+faceNECount;
315 node0=i0+(NE1-1)*N0;
316
317 out->FaceElements->Id[k]=i0+totalNECount;
318 out->FaceElements->Tag[k]=20;
319 out->FaceElements->Color[k]=i0%2+2;
320 #ifdef PASO_MPI
321 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
322 #endif
323
324 if (useElementsOnFace) {
325 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0+1;
326 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
327 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;
328 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
329 } else {
330 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0+1;
331 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
332 }
333 }
334 totalNECount+=NE0;
335 faceNECount+=NE0;
336 }
337 if (!periodic[0]) {
338 /* ** elements on boundary 001 (x1=0): */
339 #pragma omp parallel for private(i1,k,node0)
340 for (i1=0;i1<NE1;i1++) {
341 k=i1+faceNECount;
342 node0=i1*N0;
343
344 out->FaceElements->Id[k]=i1+totalNECount;
345 out->FaceElements->Tag[k]=1;
346 out->FaceElements->Color[k]=i1%2+4;
347 #ifdef PASO_MPI
348 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
349 #endif
350
351 if (useElementsOnFace) {
352 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
353 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
354 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
355 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
356 } else {
357 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
358 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
359
360 }
361 }
362 totalNECount+=NE1;
363 faceNECount+=NE1;
364
365 /* ** elements on boundary 002 (x1=1): */
366
367 #pragma omp parallel for private(i1,k,node0)
368 for (i1=0;i1<NE1;i1++) {
369 k=i1+faceNECount;
370 node0=(NE0-1)+i1*N0;
371
372 out->FaceElements->Id[k]=i1+totalNECount;
373 out->FaceElements->Tag[k]=2;
374 out->FaceElements->Color[k]=i1%2+6;
375 #ifdef PASO_MPI
376 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
377 #endif
378
379 if (useElementsOnFace) {
380 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
381 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
382 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0;
383 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
384 } else {
385 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
386 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
387 }
388 }
389 totalNECount+=NE1;
390 faceNECount+=NE1;
391 }
392 out->FaceElements->minColor=0;
393 out->FaceElements->maxColor=7;
394
395 #ifdef PASO_MPI
396 Finley_ElementFile_setDomainFlags( out->Elements );
397 Finley_ElementFile_setDomainFlags( out->FaceElements );
398 Finley_ElementFile_setDomainFlags( out->ContactElements );
399 Finley_ElementFile_setDomainFlags( out->Points );
400
401 /* reorder the degrees of freedom */
402 Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
403 #endif
404
405 /* condense the nodes: */
406 Finley_Mesh_resolveNodeIds(out);
407
408 /* prepare mesh for further calculatuions:*/
409 Finley_Mesh_prepare(out) ;
410
411 #ifdef Finley_TRACE
412 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
413 #endif
414
415 if (! Finley_noError()) {
416 Finley_Mesh_dealloc(out);
417 return NULL;
418 }
419 return out;
420 }
421 #ifdef PASO_MPI
422 Finley_Mesh* Finley_RectangularMesh_Rec4(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
423 {
424 dim_t N0,N1,NE0,NE1,i0,i1,N0t, NDOF0t, NE0_local, NDOF0,NDOF1, NFaceElements, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal, faceNECount, totalNECount,M0,M1,numDOFInternal;
425 index_t NUMNODES,k,firstNode=0, DOFcount=0, node0, node1;
426 index_t targetDomain=-1, firstNodeConstruct;
427 index_t *forwardDOF=NULL, *backwardDOF=NULL;
428 bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft=FALSE, boundaryRight=FALSE;
429 Finley_Mesh* out;
430 char name[50];
431 double time0=Finley_timer();
432
433 NE0=MAX(1,numElements[0]);
434 NE1=MAX(1,numElements[1]);
435 N0=NE0+1;
436 N1=NE1+1;
437 Paso_MPIInfo *mpi_info = NULL;
438
439 /* get MPI information */
440 mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
441 if (! Finley_noError())
442 return NULL;
443
444 /* use the serial code to generate the mesh in the 1-CPU case */
445 if( mpi_info->size==1 ){
446 out = Finley_RectangularMesh_Rec4_singleCPU(numElements,Length,periodic,order,useElementsOnFace,mpi_info);
447 return out;
448 }
449
450 if( mpi_info->rank==0 )
451 domLeft = TRUE;
452 if( mpi_info->rank==mpi_info->size-1 )
453 domRight = TRUE;
454 if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
455 domInternal = TRUE;
456
457 /* dimensions of the local subdomain */
458 domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal, &numDOFInternal );
459
460 if (N0<=N1) {
461 M0=1;
462 M1=N0;
463 } else {
464 M0=N1;
465 M1=1;
466 }
467
468 NFaceElements=0;
469 if (!periodic[0])
470 {
471 if( domLeft )
472 NFaceElements+=NE1;
473 if( domRight )
474 NFaceElements+=NE1;
475 }
476 if (!periodic[1]) {
477 NDOF1=N1;
478 NFaceElements+=2*numElementsLocal;
479 } else {
480 NDOF1=N1-1;
481 }
482
483 boundaryLeft = !domLeft || periodicLocal[0];
484 boundaryRight = !domRight || periodicLocal[1];
485 N0t = numNodesLocal + boundaryRight + boundaryLeft;
486 NDOF0t = numDOFLocal + boundaryRight + boundaryLeft;
487 firstNodeConstruct = firstNode - boundaryLeft;
488 firstNodeConstruct = firstNodeConstruct<0 ? N0-2 : firstNodeConstruct;
489
490 sprintf(name,"Rectangular %d x %d mesh",N0,N1);
491 out=Finley_Mesh_alloc( name, 2, order, mpi_info );
492
493 if (! Finley_noError()) return NULL;
494
495 out->Elements=Finley_ElementFile_alloc( Rec4, out->order, mpi_info );
496 if (useElementsOnFace) {
497 out->FaceElements=Finley_ElementFile_alloc(Rec4Face,out->order, mpi_info );
498 out->ContactElements=Finley_ElementFile_alloc(Rec4Face_Contact,out->order, mpi_info );
499 } else {
500 out->FaceElements=Finley_ElementFile_alloc(Line2,out->order, mpi_info );
501 out->ContactElements=Finley_ElementFile_alloc(Line2_Contact,out->order, mpi_info );
502 }
503 out->Points=Finley_ElementFile_alloc(Point1,out->order, mpi_info );
504 if (! Finley_noError()) {
505 Finley_Mesh_dealloc(out);
506 return NULL;
507 }
508
509 /* allocate tables: */
510 Finley_NodeFile_allocTable( out->Nodes, (numNodesLocal + numNodesExternal)*N1 );
511 Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal*NDOF1, 2*NDOF1, 0 );
512 Finley_ElementFile_allocTable(out->Elements,numElementsLocal*NE1);
513 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
514
515 if (! Finley_noError()) {
516 Finley_Mesh_dealloc(out);
517 return NULL;
518 }
519
520 /* set nodes: */
521 #pragma omp parallel for private(i0,i1,k)
522 for (k=0,i1=0;i1<N1;i1++) {
523 for ( i0=0;i0<N0t;i0++,k++) {
524 out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE((i0+firstNodeConstruct) % N0)/DBLE(N0-1)*Length[0];
525 out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
526 out->Nodes->Id[k]=k;
527 out->Nodes->Tag[k]=0;
528 out->Nodes->degreeOfFreedom[k]=i0 + (i1%NDOF1)*N0t;
529 out->Nodes->Dom[k]=NODE_INTERNAL;
530 }
531 }
532 if( boundaryLeft ){
533 for( i1=0; i1<N1; i1++ ){
534 out->Nodes->Dom[N0t*i1] = NODE_EXTERNAL;
535 out->Nodes->Dom[N0t*i1+1] = NODE_BOUNDARY;
536 }
537 }
538 if( boundaryRight ){
539 for( i1=0; i1<N1; i1++ ){
540 out->Nodes->Dom[N0t*(i1+1)-1] = NODE_EXTERNAL;
541 out->Nodes->Dom[N0t*(i1+1)-2] = NODE_BOUNDARY;
542 }
543 }
544 if( periodicLocal[0] ){
545 for( i1=0; i1<N1; i1++ ){
546 out->Nodes->degreeOfFreedom[i1*N0t+2] = out->Nodes->degreeOfFreedom[i1*N0t+1];
547 out->Nodes->Dom[N0t*i1+2] = NODE_BOUNDARY;
548 }
549 }
550
551 /* tag 'em */
552 for(i0=0; i0<N0t; i0++ )
553 {
554 out->Nodes->Tag[i0] += 10; // bottom
555 out->Nodes->Tag[i0+N0t*(N1-1)] += 20; // top
556 }
557 if( domLeft && !periodicLocal[0] )
558 for( i1=0; i1<N1; i1++ )
559 out->Nodes->Tag[i1*N0t] += 1; // left
560
561 if( domRight && !periodicLocal[0] )
562 for( i1=0; i1<N1; i1++ )
563 out->Nodes->Tag[(i1+1)*N0t-1] += 2; // right
564
565 /* form the boudary communication information */
566 forwardDOF = MEMALLOC(NDOF1,index_t);
567 backwardDOF = MEMALLOC(NDOF1,index_t);
568 if( !(mpi_info->size==2 && periodicLocal[0])){
569 if( boundaryLeft ) {
570 targetDomain = mpi_info->rank-1 < 0 ? mpi_info->size-1 : mpi_info->rank-1;
571 for( i1=0; i1<NDOF1; i1++ ){
572 forwardDOF[i1] = out->Nodes->degreeOfFreedom[i1*N0t+1];
573 backwardDOF[i1] = out->Nodes->degreeOfFreedom[i1*N0t];
574 }
575 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, forwardDOF );
576 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, backwardDOF );
577 }
578 if( boundaryRight ) {
579 targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;
580 for( i1=0; i1<NDOF1; i1++ ){
581 forwardDOF[i1] = out->Nodes->degreeOfFreedom[(i1+1)*N0t-2];
582 backwardDOF[i1] = out->Nodes->degreeOfFreedom[(i1+1)*N0t-1];
583 }
584 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, forwardDOF );
585 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, backwardDOF );
586 }
587 } else{
588 /* periodic boundary conditions with 2 domains, need to change the order in which domain 0 passes boundary data */
589 targetDomain = 1;
590
591 for( i1=0; i1<NDOF1; i1++ ){
592 forwardDOF[i1] = out->Nodes->degreeOfFreedom[(i1+1)*N0t-2];
593 backwardDOF[i1] = out->Nodes->degreeOfFreedom[(i1+1)*N0t-1];
594 }
595 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, forwardDOF );
596 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, backwardDOF );
597
598 for( i1=0; i1<NDOF1; i1++ ){
599 forwardDOF[i1] = out->Nodes->degreeOfFreedom[i1*N0t+1];
600 backwardDOF[i1] = out->Nodes->degreeOfFreedom[i1*N0t];
601 }
602 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, forwardDOF );
603 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, backwardDOF );
604 }
605 MEMFREE( forwardDOF );
606 MEMFREE( backwardDOF );
607
608 /* elements: */
609
610 #pragma omp parallel for private(i0,i1,k,node0)
611 for ( i1=0, k=0; i1<NE1; i1++) {
612 for ( i0=0; i0<numElementsLocal; i0++, k++) {
613 node0 = (periodicLocal[0] && !i0) ? i1*N0t : i1*N0t + i0 + periodicLocal[0];
614
615 out->Elements->Id[k]=k;
616 out->Elements->Tag[k]=0;
617 out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);
618 out->Elements->Dom[k]=ELEMENT_INTERNAL;
619
620 out->Elements->Nodes[INDEX2(0,k,4)]=node0;
621 out->Elements->Nodes[INDEX2(1,k,4)]=node0+1;
622 out->Elements->Nodes[INDEX2(2,k,4)]=node0+N0t+1;
623 out->Elements->Nodes[INDEX2(3,k,4)]=node0+N0t;
624 }
625 }
626 out->Elements->minColor=0;
627 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
628 if( boundaryLeft )
629 for( i1=0; i1<NE1; i1++ )
630 out->Elements->Dom[i1*numElementsLocal]=ELEMENT_BOUNDARY;
631 if( boundaryRight )
632 for( i1=0; i1<NE1; i1++ )
633 out->Elements->Dom[(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
634
635 Finley_ElementFile_setDomainFlags( out->Elements );
636
637 /* face elements: */
638 if (useElementsOnFace) {
639 NUMNODES=4;
640 } else {
641 NUMNODES=2;
642 }
643 totalNECount=numElementsLocal*NE1;
644 faceNECount=0;
645
646 if( !periodic[1] ){
647
648 /* elements on boundary 010 (x1=0): */
649 #pragma omp parallel for private(i0,k,node0)
650 for (i0=0;i0<numElementsLocal;i0++) {
651 k=i0+faceNECount;
652 node0 = (periodicLocal[0] && !i0) ? 0 : i0 + periodicLocal[0];
653
654 out->FaceElements->Id[k]=i0+totalNECount;
655 out->FaceElements->Tag[k] = 10;
656 out->FaceElements->Color[k] = i0%2;
657 out->FaceElements->Dom[k] = ELEMENT_INTERNAL;
658
659 if (useElementsOnFace) {
660 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
661 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
662 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0t+1;
663 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0t;
664 } else {
665 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
666 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
667 }
668 }
669 if( boundaryLeft ){
670 out->FaceElements->Dom[faceNECount]=ELEMENT_BOUNDARY;
671 if( periodicLocal[0] )
672 out->FaceElements->Dom[faceNECount+1]=ELEMENT_BOUNDARY;
673 }
674 if( boundaryRight )
675 out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;
676
677 totalNECount += numElementsLocal;
678 faceNECount += numElementsLocal;
679
680 /* ** elements on boundary 020 (x1=1): */
681
682 #pragma omp parallel for private(i0,k,node0)
683 for (i0=0;i0<numElementsLocal;i0++) {
684 k=i0+faceNECount;
685 node0 = (periodicLocal[0] && !i0) ? (NE1-1)*N0t : i0 + periodicLocal[0] + (NE1-1)*N0t;
686
687 out->FaceElements->Id[k]=i0+totalNECount;
688 out->FaceElements->Tag[k] = 20;
689 out->FaceElements->Color[k] = i0%2+2;
690 out->FaceElements->Dom[k] = ELEMENT_INTERNAL;
691
692 if (useElementsOnFace) {
693 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0t+1;
694 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0t;
695 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;
696 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
697 } else {
698 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0t+1;
699 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0t;
700 }
701 }
702 if( boundaryLeft ){
703 out->FaceElements->Dom[faceNECount]=ELEMENT_BOUNDARY;
704 if( periodicLocal[0] )
705 out->FaceElements->Dom[faceNECount+1]=ELEMENT_BOUNDARY;
706 }
707 if( boundaryRight )
708 out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;
709
710 totalNECount += numElementsLocal;
711 faceNECount += numElementsLocal;
712 }
713 /* ** elements on boundary 001 (x0=0): */
714 if( domLeft && !periodic[0] )
715 {
716 #pragma omp parallel for private(i1,k,node0)
717 for (i1=0;i1<NE1;i1++) {
718 k=i1+faceNECount;
719 node0=i1*N0t;
720
721 out->FaceElements->Id[k]=i1+totalNECount;
722 out->FaceElements->Tag[k] = 1;
723 out->FaceElements->Color[k] = i1%2+4;
724 out->FaceElements->Dom[k] = ELEMENT_INTERNAL;
725
726 if (useElementsOnFace) {
727 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0t;
728 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
729 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
730 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0t+1;
731 } else {
732 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0t;
733 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
734
735 }
736 }
737 totalNECount+=NE1;
738 faceNECount+=NE1;
739 }
740 /* ** elements on boundary 002 (x0=1): */
741 if( domRight && !periodic[0])
742 {
743 #pragma omp parallel for private(i1,k,node0)
744 for (i1=0;i1<NE1;i1++) {
745 k=i1+faceNECount;
746 node0=(N0t-2)+i1*N0t;
747
748 out->FaceElements->Id[k]=i1+totalNECount;
749 out->FaceElements->Tag[k] = 2;
750 out->FaceElements->Color[k] = i1%2+6;
751 out->FaceElements->Dom[k] = ELEMENT_INTERNAL;
752
753 if (useElementsOnFace) {
754 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
755 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0t+1;
756 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0t;
757 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
758 } else {
759 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
760 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0t+1;
761 }
762 }
763 totalNECount+=NE1;
764 faceNECount+=NE1;
765 }
766 out->FaceElements->minColor = 0;
767 out->FaceElements->maxColor = 7;
768
769 out->FaceElements->numElements = faceNECount;
770 Finley_ElementFile_setDomainFlags( out->FaceElements );
771
772 /* setup distribution info for other elements */
773 Finley_ElementFile_setDomainFlags( out->ContactElements );
774 Finley_ElementFile_setDomainFlags( out->Points );
775
776
777 /* reorder the degrees of freedom */
778 Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
779
780 /* condense the nodes: */
781 Finley_Mesh_resolveNodeIds(out);
782 if( !Finley_MPI_noError(mpi_info) )
783 {
784 Paso_MPIInfo_dealloc( mpi_info );
785 Finley_Mesh_dealloc(out);
786 return NULL;
787 }
788
789 /* prepare mesh for further calculatuions:*/
790 Finley_Mesh_prepare(out);
791 if( !Finley_MPI_noError(mpi_info) )
792 {
793 Paso_MPIInfo_dealloc( mpi_info );
794 Finley_Mesh_dealloc(out);
795 return NULL;
796 }
797
798 /* free up memory */
799 Paso_MPIInfo_dealloc( mpi_info );
800
801 //print_mesh_statistics( out, TRUE );
802
803 #ifdef Finley_TRACE
804 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
805 #endif
806
807 return out;
808 }
809 #endif
810
811
812

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26