/[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 1062 - (show annotations)
Mon Mar 26 06:17:53 2007 UTC (12 years, 8 months ago) by gross
File MIME type: text/plain
File size: 29733 byte(s)
reduced integration schemes are implemented now for grad, integrate, etc. Tests still to be added.
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, index_t reduced_order, int useElementsOnFace, Paso_MPIInfo *mpi_info)
131 #else
132 Finley_Mesh* Finley_RectangularMesh_Rec8(int* numElements,double* Length,int* periodic,int order, index_t reduced_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, reduced_order,mpi_info);
177 #else
178 out=Finley_Mesh_alloc(name,2,order, reduced_order);
179 #endif
180 if (! Finley_noError()) return NULL;
181
182 #ifdef PASO_MPI
183 out->Elements=Finley_ElementFile_alloc(Rec8,out->order, out->reduced_order,mpi_info);
184 if (useElementsOnFace) {
185 out->FaceElements=Finley_ElementFile_alloc(Rec8Face,out->order, out->reduced_order,mpi_info);
186 out->ContactElements=Finley_ElementFile_alloc(Rec8Face_Contact,out->order, out->reduced_order,mpi_info);
187 } else {
188 out->FaceElements=Finley_ElementFile_alloc(Line3,out->order, out->reduced_order,mpi_info);
189 out->ContactElements=Finley_ElementFile_alloc(Line3_Contact,out->order, out->reduced_order,mpi_info);
190 }
191 out->Points=Finley_ElementFile_alloc(Point1,out->order, out->reduced_order,mpi_info);
192 #else
193 out->Elements=Finley_ElementFile_alloc(Rec8,out->order, out->reduced_order);
194 if (useElementsOnFace) {
195 out->FaceElements=Finley_ElementFile_alloc(Rec8Face,out->order, out->reduced_order);
196 out->ContactElements=Finley_ElementFile_alloc(Rec8Face_Contact,out->order, out->reduced_order);
197 } else {
198 out->FaceElements=Finley_ElementFile_alloc(Line3,out->order, out->reduced_order);
199 out->ContactElements=Finley_ElementFile_alloc(Line3_Contact,out->order, out->reduced_order);
200 }
201 out->Points=Finley_ElementFile_alloc(Point1,out->order, out->reduced_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, index_t reduced_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, forwardDOF=NULL, backwardDOF=NULL, 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, reduced_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,reduced_order,mpi_info);
539
540 if (! Finley_noError()) return NULL;
541
542 out->Elements=Finley_ElementFile_alloc(Rec8,out->order, out->reduced_order,mpi_info);
543 if (useElementsOnFace) {
544 out->FaceElements=Finley_ElementFile_alloc(Rec8Face,out->order, out->reduced_order,mpi_info);
545 out->ContactElements=Finley_ElementFile_alloc(Rec8Face_Contact,out->order, out->reduced_order,mpi_info);
546 } else {
547 out->FaceElements=Finley_ElementFile_alloc(Line3,out->order, out->reduced_order,mpi_info);
548 out->ContactElements=Finley_ElementFile_alloc(Line3_Contact,out->order, out->reduced_order,mpi_info);
549 }
550 out->Points=Finley_ElementFile_alloc(Point1,out->order, out->reduced_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