/[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 969 - (show annotations)
Tue Feb 13 23:02:23 2007 UTC (12 years, 9 months ago) by ksteube
File MIME type: text/plain
File size: 29207 byte(s)
Parallelization using MPI for solution of implicit problems.

Parallelization for explicit problems has already been accomplished in
the main SVN branch.

This is incomplete and is not ready for use.


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 if ( ! Finley_Mesh_isPrepared(out)) {
452 Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");
453 }
454 } else {
455 Finley_Mesh_dealloc(out);
456 }
457 return out;
458 }
459
460 #ifdef PASO_MPI
461 Finley_Mesh* Finley_RectangularMesh_Rec8(int* numElements,double* Length,int* periodic,int order,int useElementsOnFace) {
462 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];
463 dim_t N0t, NDOF0t;
464 index_t *numForward=NULL, *numBackward=NULL;
465 index_t kk, k,i,j,node0,firstNodeConstruct,firstNode=0, DOFcount=0, targetDomain=-1,faceNEBoundary;
466 Finley_Mesh* out=NULL;
467 char name[50];
468 index_t *indexForward=NULL, *indexBackward=NULL;
469 double time0=Finley_timer();
470 bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft=FALSE, boundaryRight=FALSE;
471 Paso_MPIInfo *mpi_info = NULL;
472 index_t faceStart, numFaceLeft;
473
474 /* these are used in constructing the face elements, and give the permutation
475 of the node order of a normal element for each face */
476 index_t *facePerm;
477 index_t face0[]={0, 1, 2, 3, 4, 5, 6, 7};
478 index_t face1[]={2, 3, 0, 1, 6, 7, 4, 5};
479 index_t face2[]={3, 0, 1, 2, 7, 4, 5, 6};
480 index_t face3[]={1, 2, 3, 0, 5, 6, 7, 4};
481 index_t face0nc[]={0, 1, 4};
482 index_t face1nc[]={2, 3, 6};
483 index_t face2nc[]={3, 0, 7};
484 index_t face3nc[]={1, 2, 5};
485
486 NE0=MAX(1,numElements[0]);
487 NE1=MAX(1,numElements[1]);
488 N0=2*NE0+1;
489 N1=2*NE1+1;
490
491 /* get MPI information */
492 mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
493 if (! Finley_noError()) {
494 Finley_Mesh_dealloc(out);
495 return NULL;
496 }
497
498 // use the serial method for the single CPU case
499 if( mpi_info->size==1 ){
500 out = Finley_RectangularMesh_Rec8_singleCPU( numElements,Length,periodic,order,useElementsOnFace, mpi_info);
501 return out;
502 }
503
504 if( mpi_info->rank==0 )
505 domLeft = TRUE;
506 if( mpi_info->rank==mpi_info->size-1 )
507 domRight = TRUE;
508 if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
509 domInternal = TRUE;
510
511 /* dimensions of the local subdomain */
512 domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal, DOFBoundary );
513
514 NFaceElements=0;
515 if (!periodic[0]) {
516 NDOF0=N0;
517 NFaceElements+=(domLeft + domRight)*NE1;
518 } else {
519 NDOF0=N0-1;
520 }
521 if (!periodic[1]) {
522 NDOF1=N1;
523 NFaceElements+=2*numElementsLocal;
524 } else {
525 NDOF1=N1-1;
526 }
527
528 boundaryLeft = !domLeft || periodicLocal[0];
529 boundaryRight = !domRight || periodicLocal[1];
530 N0t = numNodesLocal + boundaryRight + boundaryLeft*2;
531 NDOF0t = numDOFLocal + boundaryRight + boundaryLeft*2;
532 firstNodeConstruct = firstNode - 2*boundaryLeft;
533 firstNodeConstruct = firstNodeConstruct<0 ? N0+firstNodeConstruct-1 : firstNodeConstruct;
534
535 /* allocate mesh: */
536
537 sprintf(name,"Rectangular %d x %d mesh",N0,N1);
538 out=Finley_Mesh_alloc(name,2,order,mpi_info);
539
540 if (! Finley_noError()) return NULL;
541
542 out->Elements=Finley_ElementFile_alloc(Rec8,out->order,mpi_info);
543 if (useElementsOnFace) {
544 out->FaceElements=Finley_ElementFile_alloc(Rec8Face,out->order,mpi_info);
545 out->ContactElements=Finley_ElementFile_alloc(Rec8Face_Contact,out->order,mpi_info);
546 } else {
547 out->FaceElements=Finley_ElementFile_alloc(Line3,out->order,mpi_info);
548 out->ContactElements=Finley_ElementFile_alloc(Line3_Contact,out->order,mpi_info);
549 }
550 out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
551
552 if (! Finley_noError()) {
553 Finley_Mesh_dealloc(out);
554 return NULL;
555 }
556
557 /* allocate tables: */
558 Finley_NodeFile_allocTable(out->Nodes,N0t*N1);
559 Finley_ElementFile_allocTable(out->Elements,numElementsLocal*NE1);
560 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
561
562 Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0t*NDOF1, NDOF1*3, 0 );
563 Finley_ElementDistribution_allocTable( out->Elements->elementDistribution, numElementsLocal*NE1, NE1*(numElementsLocal-boundaryRight*(!periodic[1])) );
564 Finley_ElementDistribution_allocTable( out->FaceElements->elementDistribution, NFaceElements, NFaceElements-2*boundaryRight*(!periodic[1]) );
565 if (! Finley_noError()) {
566 Finley_Mesh_dealloc(out);
567 return NULL;
568 }
569
570 /* set nodes: */
571 #pragma omp parallel for private(i0,i1,k)
572 for (i1=0;i1<N1;i1++) {
573 for (i0=0;i0<N0t;i0++, k++) {
574 k=i0+i1*N0t;
575 out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE((firstNodeConstruct + i0) % N0)/DBLE(N0-1)*Length[0];
576 out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
577 out->Nodes->Id[k]=i0 + i1*N0t;
578 out->Nodes->Tag[k]=0;
579 out->Nodes->Dom[k]=NODE_INTERNAL;
580 out->Nodes->degreeOfFreedom[k]=i0 + (i1%NDOF1)*N0t;
581 }
582 /* take care to make the DOF reflect the internal/external position of the DOF */
583 if( boundaryLeft ){
584 out->Nodes->Dom[i1*N0t] = NODE_EXTERNAL;
585 out->Nodes->Dom[i1*N0t+1] = NODE_EXTERNAL;
586 out->Nodes->Dom[i1*N0t+2] = NODE_BOUNDARY;
587 if( periodicLocal[0] ){
588 out->Nodes->degreeOfFreedom[i1*N0t+3] = out->Nodes->degreeOfFreedom[i1*N0t+2];
589 out->Nodes->Dom[i1*N0t+3] = NODE_BOUNDARY;
590 }
591 }
592 if( boundaryRight ){
593 out->Nodes->Dom[(i1+1)*N0t-1] = NODE_EXTERNAL;
594 out->Nodes->Dom[(i1+1)*N0t-2] = NODE_BOUNDARY;
595 out->Nodes->Dom[(i1+1)*N0t-3] = NODE_BOUNDARY;
596 }
597 }
598 /* tags for the faces: */
599 if (!periodic[1]) {
600 for (i0=0;i0<N0t;i0++) {
601 out->Nodes->Tag[i0]+=10;
602 out->Nodes->Tag[i0+N0t*(N1-1)]+=20;
603 }
604 }
605 if (domLeft && !periodic[0]) {
606 for (i1=0;i1<N1;i1++) {
607 out->Nodes->Tag[i1*N0t]+=1;
608 }
609 }
610 if (domRight && !periodic[0]){
611 for (i1=0;i1<N1;i1++) {
612 out->Nodes->Tag[(i1+1)*N0t-1]+=2;
613 }
614 }
615
616 /* setup the receive information */
617 indexBackward = TMPMEMALLOC( NDOF1*2, index_t );
618 indexForward = TMPMEMALLOC( NDOF1*2, index_t );
619 if( !(mpi_info->size==2 && periodicLocal[0])){
620 /* Backward (receive) DOF */
621 if( boundaryLeft ){
622 targetDomain = mpi_info->rank-1<0 ? mpi_info->size-1 : mpi_info->rank-1;
623
624 for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){
625 indexBackward[i1] = out->Nodes->degreeOfFreedom[i*N0t];
626 indexBackward[i1+1] = out->Nodes->degreeOfFreedom[i*N0t+1];
627 indexForward[i] = out->Nodes->degreeOfFreedom[i*N0t+2];
628 }
629
630 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );
631 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexBackward );
632 }
633
634 if( boundaryRight ){
635 targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;
636
637 for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){
638 indexBackward[i] = out->Nodes->degreeOfFreedom[(i+1)*N0t-1];
639 indexForward[i1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-3];
640 indexForward[i1+1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-2];
641 }
642
643 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexForward );
644 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexBackward );
645 }
646 }
647 else{
648 targetDomain = 1;
649
650 for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){
651 indexBackward[i] = out->Nodes->degreeOfFreedom[(i+1)*N0t-1];
652 indexForward[i1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-3];
653 indexForward[i1+1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-2];
654 }
655
656 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexForward );
657 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexBackward );
658
659 for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){
660 indexBackward[i1] = out->Nodes->degreeOfFreedom[i*N0t];
661 indexBackward[i1+1] = out->Nodes->degreeOfFreedom[i*N0t+1];
662 indexForward[i] = out->Nodes->degreeOfFreedom[i*N0t+2];
663 }
664
665 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );
666 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexBackward );
667 }
668 TMPMEMFREE( indexBackward );
669 TMPMEMFREE( indexForward );
670
671 /* set the elements: */
672 #pragma omp parallel for private(i0,i1,k,node0)
673 for (i1=0;i1<NE1;i1++) {
674 for (i0=0;i0<numElementsLocal;i0++) {
675 k=i0+numElementsLocal*i1;
676 node0 = (periodicLocal[0] && !i0) ? 2*i1*N0t : 2*i0+2*i1*N0t+periodicLocal[0] ;
677
678 out->Elements->Id[k]=((firstNodeConstruct/2+i0)%NE0) * NE1 + i1;
679 out->Elements->Tag[k]=0;
680 out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);
681 out->Elements->Dom[k] = ELEMENT_INTERNAL;
682
683 out->Elements->Nodes[INDEX2(0,k,8)]=node0;
684 out->Elements->Nodes[INDEX2(1,k,8)]=node0+2;
685 out->Elements->Nodes[INDEX2(2,k,8)]=node0+2*N0t+2;
686 out->Elements->Nodes[INDEX2(3,k,8)]=node0+2*N0t;
687 out->Elements->Nodes[INDEX2(4,k,8)]=node0+1;
688 out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0t+2;
689 out->Elements->Nodes[INDEX2(6,k,8)]=node0+2*N0t+1;
690 out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0t;
691 }
692 if( boundaryLeft )
693 out->Elements->Dom[i1*numElementsLocal] = ELEMENT_BOUNDARY;
694 if( boundaryRight )
695 out->Elements->Dom[(i1+1)*numElementsLocal-1] = ELEMENT_BOUNDARY;
696 }
697 out->Elements->minColor=0;
698 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
699
700 /* set the distribution information for the elements */
701 Finley_ElementFile_setDomainFlags( out->Elements );
702
703 /* face elements: */
704 if (useElementsOnFace) {
705 NUMNODES=8;
706 } else {
707 NUMNODES=3;
708 }
709
710 Finley_ElementFile_setDomainFlags( out->FaceElements );
711 totalNECount=numElementsLocal*NE1;
712 faceNECount=0;
713 faceNEBoundary = totalNECount + numElementsInternal*(!periodic[1])*2 + NE1*(!periodic[0])*(domLeft+domRight);
714
715 numFaceLeft = domLeft*(!periodic[0])*NE1;
716 faceStart = out->FaceElements->elementDistribution->vtxdist[mpi_info->rank];
717 if (!periodic[1]) {
718 /* ** elements on boundary 010 (x1=0): */
719
720 #pragma omp parallel for private(i0,k,kk,facePerm)
721 for (i0=0;i0<numElementsLocal; i0++) {
722 k=i0+faceNECount;
723 kk=i0;
724 facePerm = useElementsOnFace ? face0 : face0nc;
725
726 out->FaceElements->Id[k]=faceStart + numFaceLeft + i0*2;
727 out->FaceElements->Tag[k]=10;
728 out->FaceElements->Color[k]=i0%2;
729 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
730
731 for( j=0; j<NUMNODES; j++ )
732 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
733 }
734 if( boundaryLeft ){
735 out->FaceElements->Dom[faceNECount]=ELEMENT_BOUNDARY;
736 if( periodicLocal[0] )
737 out->FaceElements->Dom[faceNECount+1]=ELEMENT_BOUNDARY;
738 }
739 if( boundaryRight ){
740 out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;
741 out->FaceElements->Id[faceNECount+numElementsLocal-1]=out->FaceElements->elementDistribution->vtxdist[mpi_info->rank+1];
742 }
743 totalNECount+=numElementsLocal;
744 faceNECount+=numElementsLocal;
745
746 /* ** elements on boundary 020 (x1=1): */
747
748 #pragma omp parallel for private(i0,k,kk,facePerm)
749 for (i0=0;i0<numElementsLocal;i0++) {
750 k=i0+faceNECount;
751 kk=i0 + numElementsLocal*(NE1-1);
752 facePerm = useElementsOnFace ? face1 : face1nc;
753
754 out->FaceElements->Id[k]=faceStart + numFaceLeft + i0*2+1;
755 out->FaceElements->Tag[k]=20;
756 out->FaceElements->Color[k]=i0%2+2;
757 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
758
759 for( j=0; j<NUMNODES; j++ )
760 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
761 }
762 if( boundaryLeft ){
763 out->FaceElements->Dom[faceNECount]=ELEMENT_BOUNDARY;
764 if( periodicLocal[0] )
765 out->FaceElements->Dom[faceNECount+1]=ELEMENT_BOUNDARY;
766 }
767 if( boundaryRight ){
768 out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;
769 out->FaceElements->Id[faceNECount+numElementsLocal-1]=out->FaceElements->elementDistribution->vtxdist[mpi_info->rank+1]+1;
770 }
771 totalNECount+=numElementsLocal;
772 faceNECount+=numElementsLocal;
773 }
774 if (domLeft && !periodicLocal[0]) {
775 /* ** elements on boundary 001 (x0=0): */
776
777 #pragma omp parallel for private(i1,k,kk,facePerm)
778 for (i1=0;i1<NE1;i1++) {
779 k=i1+faceNECount;
780 kk=i1*numElementsLocal;
781 facePerm = useElementsOnFace ? face2 : face2nc;
782
783 out->FaceElements->Id[k]=faceStart + i1;
784 out->FaceElements->Tag[k]=1;
785 out->FaceElements->Color[k]=(i1%2)+4;
786 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
787
788 for( j=0; j<NUMNODES; j++ )
789 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
790 }
791 totalNECount+=NE1;
792 faceNECount+=NE1;
793
794 /* ** elements on boundary 002 (x0=1): */
795 }
796 if (domRight && !periodicLocal[1]) {
797 #pragma omp parallel for private(i1,k,kk,facePerm)
798 for (i1=0;i1<NE1;i1++) {
799 k=i1+faceNECount;
800 kk=(i1+1)*numElementsLocal-1;
801 facePerm = useElementsOnFace ? face3 : face3nc;
802
803 out->FaceElements->Id[k]=faceStart + numElementsLocal*2*(!periodic[1])+i1;
804 out->FaceElements->Tag[k]=2;
805 out->FaceElements->Color[k]=(i1%2)+6;
806 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
807
808 for( j=0; j<NUMNODES; j++ )
809 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
810 }
811 totalNECount+=NE1;
812 faceNECount+=NE1;
813 }
814 out->FaceElements->minColor=0;
815 out->FaceElements->maxColor=7;
816 out->FaceElements->numElements = faceNECount;
817
818 /* setup distribution info for other elements */
819 Finley_ElementFile_setDomainFlags( out->ContactElements );
820 Finley_ElementFile_setDomainFlags( out->Points );
821
822 Finley_Mesh_prepareElementDistribution( out );
823
824 /* reorder the degrees of freedom */
825 Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
826
827 /* condense the nodes: */
828 Finley_Mesh_resolveNodeIds(out);
829 if ( !Finley_MPI_noError( mpi_info )) {
830 Paso_MPIInfo_dealloc( mpi_info );
831 Finley_Mesh_dealloc(out);
832 return NULL;
833 }
834
835 /* prepare mesh for further calculatuions:*/
836 Finley_Mesh_prepare(out) ;
837
838 #ifdef Finley_TRACE
839 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
840 #endif
841
842 if ( !Finley_MPI_noError( mpi_info )) {
843 Paso_MPIInfo_dealloc( mpi_info );
844 Finley_Mesh_dealloc(out);
845 return NULL;
846 }
847 if (Finley_noError()) {
848 if (! Finley_Mesh_isPrepared(out)) {
849 Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");
850 }
851 }
852 return out;
853 }
854 #endif

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26