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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26