/[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 787 - (show annotations)
Wed Jul 26 01:46:45 2006 UTC (13 years, 3 months ago) by bcumming
File MIME type: text/plain
File size: 27223 byte(s)
MPI update
Each element (normal elements, faceElements, ContactElements and points)
are now assigned a unique global id to streamline per-element
calculations and file IO of element data.



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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26