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

Annotation of /temp_trunk_copy/finley/src/Mesh_rec4.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 964 - (hide annotations)
Tue Feb 13 05:10:26 2007 UTC (12 years, 9 months ago) by gross
Original Path: trunk/finley/src/Mesh_rec4.c
File MIME type: text/plain
File size: 27562 byte(s)
The set/getRefVal functions of Data objects have been removed (mainly to avoid later problems with MPI).
Moreover, a faster access to the reference id of samples has been introduced. I don't think that anybody will
profit form this at this stage but it will allow a faster dump of data objects.


1 jgs 150 /*
2 elspeth 616 ************************************************************
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 jgs 150 */
12    
13 jgs 82 /**************************************************************/
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 bcumming 730 #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 bcumming 751 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 bcumming 730 {
48     index_t i0;
49     dim_t numNodesGlobal = numElementsGlobal+1;
50    
51     (*numNodesLocal) = domain_MODdim( rank, size, numNodesGlobal );
52 bcumming 751
53 bcumming 730 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 bcumming 751
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 bcumming 730 DOFExternal[0] = nodesExternal[0];
120     DOFExternal[1] = nodesExternal[1];
121     }
122 bcumming 751 #endif
123 bcumming 730
124 bcumming 782 #ifdef PASO_MPI
125     Finley_Mesh* Finley_RectangularMesh_Rec4_singleCPU(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace,Paso_MPIInfo *mpi_info)
126     #else
127 bcumming 730 Finley_Mesh* Finley_RectangularMesh_Rec4(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
128 bcumming 782 #endif
129 bcumming 730 {
130 jgs 153 dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES,M0,M1;
131 jgs 123 index_t k,node0;
132 jgs 82 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 jgs 153 if (N0<=N1) {
141     M0=1;
142     M1=N0;
143     } else {
144     M0=N1;
145     M1=1;
146     }
147    
148 jgs 82 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 bcumming 782 #ifdef PASO_MPI
166     out=Finley_Mesh_alloc(name,2,order,mpi_info);
167    
168     if (! Finley_noError()) return NULL;
169    
170     out->Elements=Finley_ElementFile_alloc(Rec4,out->order,mpi_info);
171     if (useElementsOnFace) {
172     out->FaceElements=Finley_ElementFile_alloc(Rec4Face,out->order,mpi_info);
173     out->ContactElements=Finley_ElementFile_alloc(Rec4Face_Contact,out->order,mpi_info);
174     } else {
175     out->FaceElements=Finley_ElementFile_alloc(Line2,out->order,mpi_info);
176     out->ContactElements=Finley_ElementFile_alloc(Line2_Contact,out->order,mpi_info);
177     }
178     out->Points=Finley_ElementFile_alloc(Point1,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 bcumming 787
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 bcumming 782 #else
194 jgs 82 out=Finley_Mesh_alloc(name,2,order);
195 bcumming 730
196 jgs 150 if (! Finley_noError()) return NULL;
197 jgs 82
198     out->Elements=Finley_ElementFile_alloc(Rec4,out->order);
199     if (useElementsOnFace) {
200     out->FaceElements=Finley_ElementFile_alloc(Rec4Face,out->order);
201     out->ContactElements=Finley_ElementFile_alloc(Rec4Face_Contact,out->order);
202     } else {
203     out->FaceElements=Finley_ElementFile_alloc(Line2,out->order);
204     out->ContactElements=Finley_ElementFile_alloc(Line2_Contact,out->order);
205     }
206     out->Points=Finley_ElementFile_alloc(Point1,out->order);
207 bcumming 730
208 jgs 150 if (! Finley_noError()) {
209 jgs 82 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 bcumming 782 #endif
218 jgs 150 if (! Finley_noError()) {
219 jgs 82 Finley_Mesh_dealloc(out);
220     return NULL;
221     }
222     /* set nodes: */
223 jgs 153 #pragma omp parallel for private(i0,i1,k)
224 jgs 82 for (i1=0;i1<N1;i1++) {
225     for (i0=0;i0<N0;i0++) {
226 jgs 153 k=M0*i0+M1*i1;
227 jgs 82 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 jgs 153 out->Nodes->Id[k]=i0+N0*i1;
230 jgs 82 out->Nodes->Tag[k]=0;
231 jgs 153 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1);
232 bcumming 782 #ifdef PASO_MPI
233     out->Nodes->Dom[k]=NODE_INTERNAL;
234     #endif
235 jgs 82 }
236     }
237     /* tags for the faces: */
238     if (!periodic[1]) {
239     for (i0=0;i0<N0;i0++) {
240 jgs 153 out->Nodes->Tag[M0*i0+M1*0]+=10;
241     out->Nodes->Tag[M0*i0+M1*(N1-1)]+=20;
242 jgs 82 }
243     }
244     if (!periodic[0]) {
245     for (i1=0;i1<N1;i1++) {
246 jgs 153 out->Nodes->Tag[M0*0+M1*i1]+=1;
247     out->Nodes->Tag[M0*(N0-1)+M1*i1]+=2;
248 jgs 82 }
249     }
250 jgs 153
251 jgs 82 /* 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 bcumming 782 #ifdef PASO_MPI
263     out->Elements->Dom[k]=ELEMENT_INTERNAL;
264     #endif
265 jgs 82
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 jgs 123 out->Elements->minColor=0;
274     out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
275 jgs 82
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 bcumming 782 #ifdef PASO_MPI
297     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
298     #endif
299 jgs 82
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 bcumming 782 #ifdef PASO_MPI
324     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
325     #endif
326 jgs 82
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 bcumming 782 #ifdef PASO_MPI
351     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
352     #endif
353 jgs 82
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 bcumming 782 #ifdef PASO_MPI
379     out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
380     #endif
381 jgs 82
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 jgs 123 out->FaceElements->minColor=0;
396     out->FaceElements->maxColor=7;
397 bcumming 782
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 jgs 82
404 bcumming 782 /* reorder the degrees of freedom */
405     Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
406     #endif
407    
408 jgs 82 /* condense the nodes: */
409     Finley_Mesh_resolveNodeIds(out);
410    
411     /* prepare mesh for further calculatuions:*/
412     Finley_Mesh_prepare(out) ;
413    
414 jgs 149 #ifdef Finley_TRACE
415 jgs 82 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
416 jgs 149 #endif
417 gross 964 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 jgs 82 Finley_Mesh_dealloc(out);
423     }
424     return out;
425     }
426 bcumming 782 #ifdef PASO_MPI
427     Finley_Mesh* Finley_RectangularMesh_Rec4(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
428 bcumming 730 {
429 bcumming 782 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 bcumming 730 Finley_Mesh* out;
435     char name[50];
436     double time0=Finley_timer();
437 bcumming 787 index_t faceStart, numFaceLeft, i;
438 bcumming 730
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 bcumming 782 /* 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,useElementsOnFace,mpi_info);
453     return out;
454     }
455    
456 bcumming 730 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 bcumming 751 domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal, &numDOFInternal );
465 bcumming 730
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 bcumming 751 NFaceElements+=2*numElementsLocal;
485 bcumming 730 } else {
486     NDOF1=N1-1;
487     }
488 bcumming 782
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 bcumming 730
496     sprintf(name,"Rectangular %d x %d mesh",N0,N1);
497     out=Finley_Mesh_alloc( name, 2, order, mpi_info );
498    
499     if (! Finley_noError()) return NULL;
500    
501     out->Elements=Finley_ElementFile_alloc( Rec4, out->order, mpi_info );
502     if (useElementsOnFace) {
503     out->FaceElements=Finley_ElementFile_alloc(Rec4Face,out->order, mpi_info );
504     out->ContactElements=Finley_ElementFile_alloc(Rec4Face_Contact,out->order, mpi_info );
505     } else {
506     out->FaceElements=Finley_ElementFile_alloc(Line2,out->order, mpi_info );
507     out->ContactElements=Finley_ElementFile_alloc(Line2_Contact,out->order, mpi_info );
508     }
509     out->Points=Finley_ElementFile_alloc(Point1,out->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 bcumming 787 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 bcumming 730 if (! Finley_noError()) {
525     Finley_Mesh_dealloc(out);
526     return NULL;
527     }
528 bcumming 782
529 bcumming 730 /* set nodes: */
530 bcumming 782 #pragma omp parallel for private(i0,i1,k)
531 bcumming 787 for (i1=0;i1<N1;i1++) {
532 bcumming 782 for ( i0=0;i0<N0t;i0++,k++) {
533 bcumming 787 k=i0+i1*N0t;
534 bcumming 782 out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE((i0+firstNodeConstruct) % N0)/DBLE(N0-1)*Length[0];
535 bcumming 730 out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
536 bcumming 782 out->Nodes->Id[k]=k;
537 bcumming 730 out->Nodes->Tag[k]=0;
538 bcumming 782 out->Nodes->degreeOfFreedom[k]=i0 + (i1%NDOF1)*N0t;
539     out->Nodes->Dom[k]=NODE_INTERNAL;
540 bcumming 730 }
541     }
542 bcumming 782 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 bcumming 730 /* tag 'em */
562 bcumming 782 for(i0=0; i0<N0t; i0++ )
563 bcumming 730 {
564     out->Nodes->Tag[i0] += 10; // bottom
565 bcumming 782 out->Nodes->Tag[i0+N0t*(N1-1)] += 20; // top
566 bcumming 730 }
567     if( domLeft && !periodicLocal[0] )
568     for( i1=0; i1<N1; i1++ )
569 bcumming 782 out->Nodes->Tag[i1*N0t] += 1; // left
570 bcumming 730
571     if( domRight && !periodicLocal[0] )
572     for( i1=0; i1<N1; i1++ )
573 bcumming 782 out->Nodes->Tag[(i1+1)*N0t-1] += 2; // right
574 bcumming 730
575 bcumming 782 /* 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 bcumming 730
608 bcumming 782 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 bcumming 730
618 bcumming 782 /* elements: */
619 bcumming 730
620 bcumming 782 #pragma omp parallel for private(i0,i1,k,node0)
621 bcumming 787 for ( i1=0; i1<NE1; i1++) {
622 bcumming 782 for ( i0=0; i0<numElementsLocal; i0++, k++) {
623 bcumming 787 k=i1*numElementsLocal + i0;
624 bcumming 782 node0 = (periodicLocal[0] && !i0) ? i1*N0t : i1*N0t + i0 + periodicLocal[0];
625 bcumming 730
626 bcumming 787 out->Elements->Id[k]=((firstNodeConstruct+i0)%NE0) * NE1 + i1;
627 bcumming 730 out->Elements->Tag[k]=0;
628 bcumming 782 out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);
629     out->Elements->Dom[k]=ELEMENT_INTERNAL;
630 bcumming 730
631     out->Elements->Nodes[INDEX2(0,k,4)]=node0;
632     out->Elements->Nodes[INDEX2(1,k,4)]=node0+1;
633 bcumming 782 out->Elements->Nodes[INDEX2(2,k,4)]=node0+N0t+1;
634     out->Elements->Nodes[INDEX2(3,k,4)]=node0+N0t;
635 bcumming 730 }
636     }
637     out->Elements->minColor=0;
638 bcumming 782 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 bcumming 730
648     /* face elements: */
649     if (useElementsOnFace) {
650     NUMNODES=4;
651     } else {
652     NUMNODES=2;
653     }
654     totalNECount=numElementsLocal*NE1;
655     faceNECount=0;
656    
657 bcumming 787 numFaceLeft = domLeft*(!periodic[0])*NE1;
658     faceStart = out->FaceElements->elementDistribution->vtxdist[mpi_info->rank];
659 bcumming 782 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 bcumming 787 out->FaceElements->Id[k]=faceStart + numFaceLeft + i0*2;
668 bcumming 782 out->FaceElements->Tag[k] = 10;
669     out->FaceElements->Color[k] = i0%2;
670     out->FaceElements->Dom[k] = ELEMENT_INTERNAL;
671 bcumming 730
672 bcumming 782 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 bcumming 787 if( boundaryRight ){
688 bcumming 782 out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;
689 bcumming 787 out->FaceElements->Id[faceNECount+numElementsLocal-1]=out->FaceElements->elementDistribution->vtxdist[mpi_info->rank+1];
690     }
691 bcumming 730
692 bcumming 782 totalNECount += numElementsLocal;
693     faceNECount += numElementsLocal;
694    
695     /* ** elements on boundary 020 (x1=1): */
696 bcumming 730
697 bcumming 782 #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 bcumming 730
702 bcumming 787 out->FaceElements->Id[k]=faceStart + numFaceLeft + i0*2+1;
703 bcumming 751 out->FaceElements->Tag[k] = 20;
704 bcumming 782 out->FaceElements->Color[k] = i0%2+2;
705     out->FaceElements->Dom[k] = ELEMENT_INTERNAL;
706 bcumming 730
707 bcumming 751 if (useElementsOnFace) {
708 bcumming 782 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 bcumming 751 } else {
713 bcumming 782 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0t+1;
714     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0t;
715 bcumming 751 }
716 bcumming 782 }
717     if( boundaryLeft ){
718     out->FaceElements->Dom[faceNECount]=ELEMENT_BOUNDARY;
719     if( periodicLocal[0] )
720     out->FaceElements->Dom[faceNECount+1]=ELEMENT_BOUNDARY;
721 bcumming 751 }
722 bcumming 787 if( boundaryRight ){
723 bcumming 782 out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;
724 bcumming 787 out->FaceElements->Id[faceNECount+numElementsLocal-1]=out->FaceElements->elementDistribution->vtxdist[mpi_info->rank+1]+1;
725     }
726 bcumming 730
727 bcumming 782 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 bcumming 730
738 bcumming 787 out->FaceElements->Id[k]=faceStart + i1;
739 bcumming 782 out->FaceElements->Tag[k] = 1;
740     out->FaceElements->Color[k] = i1%2+4;
741     out->FaceElements->Dom[k] = ELEMENT_INTERNAL;
742    
743 bcumming 751 if (useElementsOnFace) {
744 bcumming 782 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0t;
745 bcumming 751 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
746     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
747 bcumming 782 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0t+1;
748 bcumming 751 } else {
749 bcumming 782 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0t;
750 bcumming 751 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
751 bcumming 782
752 bcumming 751 }
753 bcumming 782 }
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 bcumming 730
765 bcumming 787 out->FaceElements->Id[k]=faceStart + numElementsLocal*2*(!periodic[1])+i1;
766 bcumming 782 out->FaceElements->Tag[k] = 2;
767     out->FaceElements->Color[k] = i1%2+6;
768     out->FaceElements->Dom[k] = ELEMENT_INTERNAL;
769 bcumming 730
770 bcumming 751 if (useElementsOnFace) {
771     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
772 bcumming 782 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0t+1;
773     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0t;
774 bcumming 751 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
775     } else {
776     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
777 bcumming 782 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0t+1;
778 bcumming 751 }
779     }
780 bcumming 782 totalNECount+=NE1;
781     faceNECount+=NE1;
782 bcumming 751 }
783 bcumming 730 out->FaceElements->minColor = 0;
784 bcumming 782 out->FaceElements->maxColor = 7;
785 bcumming 730
786 bcumming 782 out->FaceElements->numElements = faceNECount;
787     Finley_ElementFile_setDomainFlags( out->FaceElements );
788 bcumming 730
789 bcumming 782 /* setup distribution info for other elements */
790     Finley_ElementFile_setDomainFlags( out->ContactElements );
791     Finley_ElementFile_setDomainFlags( out->Points );
792    
793 bcumming 787 Finley_Mesh_prepareElementDistribution( out );
794 bcumming 782
795     /* reorder the degrees of freedom */
796     Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
797 bcumming 751
798 bcumming 730 /* condense the nodes: */
799 bcumming 782 Finley_Mesh_resolveNodeIds(out);
800     if( !Finley_MPI_noError(mpi_info) )
801     {
802 bcumming 751 Paso_MPIInfo_dealloc( mpi_info );
803     Finley_Mesh_dealloc(out);
804     return NULL;
805 bcumming 782 }
806 bcumming 751
807 bcumming 730 /* prepare mesh for further calculatuions:*/
808 bcumming 782 Finley_Mesh_prepare(out);
809 bcumming 751 if( !Finley_MPI_noError(mpi_info) )
810     {
811     Paso_MPIInfo_dealloc( mpi_info );
812     Finley_Mesh_dealloc(out);
813     return NULL;
814     }
815 bcumming 782
816     /* free up memory */
817 bcumming 730 Paso_MPIInfo_dealloc( mpi_info );
818    
819 bcumming 782 //print_mesh_statistics( out, TRUE );
820    
821     #ifdef Finley_TRACE
822     printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
823     #endif
824 gross 964 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 bcumming 782 return out;
830 bcumming 730 }
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