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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26