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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26