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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 730 - (hide annotations)
Mon May 15 04:03:49 2006 UTC (13 years, 5 months ago) by bcumming
File MIME type: text/plain
File size: 37208 byte(s)


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     static void Finley_Mesh_printDistributed( Finley_Mesh *in )
34     {
35     dim_t k, i0, NUMNODES;
36     Finley_Mesh *out=in;
37    
38     NUMNODES = in->FaceElements->ReferenceElement->Type->numNodes;
39    
40     printf( "\n============NODES=============\n" );
41     for( k=0; k<in->Nodes->numNodes; k++ )
42     printf( "\tId %d\tDOF %d\tcoord [%3g%3g]\n", out->Nodes->Id[k], out->Nodes->degreeOfFreedom[k] , out->Nodes->Coordinates[INDEX2(0,k,2)], out->Nodes->Coordinates[INDEX2(1,k,2)] );
43     for( k=0; k<out->Nodes->degreeOfFreedomDistribution->numNeighbours; k++ )
44     {
45     if( out->Nodes->degreeOfFreedomDistribution->neighbours[k]>=0 )
46     {
47     printf( "\t%d boundary DOF { ", out->Nodes->degreeOfFreedomDistribution->edges[k]->numForward );
48     for( i0=0; i0<out->Nodes->degreeOfFreedomDistribution->edges[k]->numForward; i0++ )
49     printf( "%d ", out->Nodes->degreeOfFreedomDistribution->edges[k]->indexForward[i0] );
50     printf("} to %d\n", out->Nodes->degreeOfFreedomDistribution->neighbours[k] );
51     printf( "\t%d external DOF { ", out->Nodes->degreeOfFreedomDistribution->edges[k]->numBackward );
52     for( i0=0; i0<out->Nodes->degreeOfFreedomDistribution->edges[k]->numBackward; i0++ )
53     printf( "%d ", out->Nodes->degreeOfFreedomDistribution->edges[k]->indexBackward[i0] );
54     printf("} from %d\n", out->Nodes->degreeOfFreedomDistribution->neighbours[k] );
55     }
56     }
57    
58     printf( "\n============ELEMENTS=============\n");
59     for( k=0; k<in->Elements->elementDistribution->numInternal; k++ )
60     {
61     printf( "I\tId %d : nodes [%d %d %d %d]->DOF [%d %d %d %d]\n", out->Elements->Id[k], out->Elements->Nodes[INDEX2(0,k,4)], out->Elements->Nodes[INDEX2(1,k,4)], out->Elements->Nodes[INDEX2(2,k,4)], out->Elements->Nodes[INDEX2(3,k,4)], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(0,k,4)]], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(1,k,4)]], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(2,k,4)]], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(3,k,4)]] );
62     }
63    
64     for( k=in->Elements->elementDistribution->numInternal; k<in->Elements->elementDistribution->numInternal+in->Elements->elementDistribution->numBoundary; k++ )
65     {
66     printf( "B\tId %d : nodes [%d %d %d %d]->DOF [%d %d %d %d]\n", out->Elements->Id[k], out->Elements->Nodes[INDEX2(0,k,4)], out->Elements->Nodes[INDEX2(1,k,4)], out->Elements->Nodes[INDEX2(2,k,4)], out->Elements->Nodes[INDEX2(3,k,4)], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(0,k,4)]], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(1,k,4)]], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(2,k,4)]], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(3,k,4)]] );
67     }
68    
69     for( k=in->FaceElements->elementDistribution->numInternal; k<in->FaceElements->elementDistribution->numInternal+in->FaceElements->elementDistribution->numBoundary; k++ )
70     {
71     if( NUMNODES==4 )
72     printf( "F\tId %d : nodes [%d %d %d %d]->DOF [%d %d %d %d]\n", out->FaceElements->Id[k], out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)], out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)], out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)], out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)], out->Nodes->degreeOfFreedom[out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]], out->Nodes->degreeOfFreedom[out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]], out->Nodes->degreeOfFreedom[out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]], out->Nodes->degreeOfFreedom[out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]] );
73     else
74     printf( "F\tId %d : nodes [%d %d]->DOF [%d %d]\n", out->FaceElements->Id[k], out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)], out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)], out->Nodes->degreeOfFreedom[out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]], out->Nodes->degreeOfFreedom[out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]] );
75     }
76    
77     }
78     #endif
79    
80    
81    
82     /* get the number of nodes/elements for domain with rank=rank, of size processors
83     where n is the total number of nodes/elements in the global domain */
84     static index_t domain_MODdim( index_t rank, index_t size, index_t n )
85     {
86     rank = size-rank-1;
87    
88     if( rank < n%size )
89     return (index_t)floor(n/size)+1;
90     return (index_t)floor(n/size);
91     }
92    
93    
94     /* Determines the number of nodes/elements etc along an axis which is numElementsGlobal long for domain rank */
95     /* A bit messy, but it only has to be done once... */
96     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 )
97     {
98     index_t i0;
99     dim_t numNodesGlobal = numElementsGlobal+1;
100    
101     (*numNodesLocal) = domain_MODdim( rank, size, numNodesGlobal );
102    
103     numElementsLocal[0] = numNodesLocal[0]+1;
104     periodicLocal[0] = periodicLocal[1] = FALSE;
105     nodesExternal[0] = nodesExternal[1] = 1;
106     if( periodic )
107     {
108     if( size==1 )
109     {
110     numElementsLocal[0] = numElementsGlobal;
111     nodesExternal[0] = nodesExternal[1] = 0;
112     periodicLocal[0] = periodicLocal[1] = TRUE;
113     }
114     else
115     {
116     if( rank==0 )
117     {
118     periodicLocal[0] = TRUE;
119     numNodesLocal[0]++;
120     }
121     if( rank==(size-1) )
122     {
123     periodicLocal[1] = TRUE;
124     numNodesLocal[0]--;
125     numElementsLocal[0]--;
126     }
127     }
128     }
129     else if( !periodic )
130     {
131     if( rank==0 ){
132     nodesExternal[0]--;
133     numElementsLocal[0]--;
134     }
135     if( rank==(size-1) )
136     {
137     nodesExternal[1]--;
138     numElementsLocal[0]--;
139     }
140     }
141     numNodesExternal[0] = nodesExternal[0]+nodesExternal[1];
142    
143     numElementsInternal[0] = numElementsLocal[0];
144     if( (rank==0) && (rank==size-1) );
145    
146     else if( !periodic && ( (rank==0) ^ (rank==size-1) ) )
147     numElementsInternal[0] -= 1;
148     else
149     numElementsInternal[0] -= 2;
150    
151     firstNode[0] = 0;
152     for( i0=0; i0<rank; i0++ )
153     firstNode[0] += domain_MODdim( i0, size, numNodesGlobal );
154    
155     numDOFLocal[0] = numNodesLocal[0];
156     if( periodicLocal[0] )
157     {
158     numDOFLocal[0]--;
159     }
160     DOFExternal[0] = nodesExternal[0];
161     DOFExternal[1] = nodesExternal[1];
162    
163     /* some debugging printf statements */
164     //printf( "rank/size = %d/%d\nNodes : %d Local, %d External[%d %d], First = %d\nElements : %d Local\nDOF : %d Local, External [%d %d]\nperiodicLocal [%d %d]\n\n", rank, size, *numNodesLocal, *numNodesExternal, nodesExternal[0], nodesExternal[1], *firstNode, *numElementsLocal, *numDOFLocal, DOFExternal[0], DOFExternal[1], periodicLocal[0], periodicLocal[1] );
165     }
166    
167    
168     Finley_Mesh* Finley_RectangularMesh_Rec4(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
169     #ifndef PASO_MPI
170     {
171 jgs 153 dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES,M0,M1;
172 jgs 123 index_t k,node0;
173 jgs 82 Finley_Mesh* out;
174     char name[50];
175     double time0=Finley_timer();
176     NE0=MAX(1,numElements[0]);
177     NE1=MAX(1,numElements[1]);
178     N0=NE0+1;
179     N1=NE1+1;
180    
181 jgs 153 if (N0<=N1) {
182     M0=1;
183     M1=N0;
184     } else {
185     M0=N1;
186     M1=1;
187     }
188    
189 jgs 82 NFaceElements=0;
190     if (!periodic[0]) {
191     NDOF0=N0;
192     NFaceElements+=2*NE1;
193     } else {
194     NDOF0=N0-1;
195     }
196     if (!periodic[1]) {
197     NDOF1=N1;
198     NFaceElements+=2*NE0;
199     } else {
200     NDOF1=N1-1;
201     }
202    
203     /* allocate mesh: */
204    
205     sprintf(name,"Rectangular %d x %d mesh",N0,N1);
206     out=Finley_Mesh_alloc(name,2,order);
207 bcumming 730
208 jgs 150 if (! Finley_noError()) return NULL;
209 jgs 82
210     out->Elements=Finley_ElementFile_alloc(Rec4,out->order);
211     if (useElementsOnFace) {
212     out->FaceElements=Finley_ElementFile_alloc(Rec4Face,out->order);
213     out->ContactElements=Finley_ElementFile_alloc(Rec4Face_Contact,out->order);
214     } else {
215     out->FaceElements=Finley_ElementFile_alloc(Line2,out->order);
216     out->ContactElements=Finley_ElementFile_alloc(Line2_Contact,out->order);
217     }
218     out->Points=Finley_ElementFile_alloc(Point1,out->order);
219 bcumming 730
220 jgs 150 if (! Finley_noError()) {
221 jgs 82 Finley_Mesh_dealloc(out);
222     return NULL;
223     }
224    
225     /* allocate tables: */
226     Finley_NodeFile_allocTable(out->Nodes,N0*N1);
227     Finley_ElementFile_allocTable(out->Elements,NE0*NE1);
228     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
229 bcumming 730
230 jgs 150 if (! Finley_noError()) {
231 jgs 82 Finley_Mesh_dealloc(out);
232     return NULL;
233     }
234     /* set nodes: */
235 jgs 153
236     #pragma omp parallel for private(i0,i1,k)
237 jgs 82 for (i1=0;i1<N1;i1++) {
238     for (i0=0;i0<N0;i0++) {
239 jgs 153 k=M0*i0+M1*i1;
240 jgs 82 out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(i0)/DBLE(N0-1)*Length[0];
241     out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
242 jgs 153 out->Nodes->Id[k]=i0+N0*i1;
243 jgs 82 out->Nodes->Tag[k]=0;
244 jgs 153 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1);
245 jgs 82 }
246     }
247     /* tags for the faces: */
248     if (!periodic[1]) {
249     for (i0=0;i0<N0;i0++) {
250 jgs 153 out->Nodes->Tag[M0*i0+M1*0]+=10;
251     out->Nodes->Tag[M0*i0+M1*(N1-1)]+=20;
252 jgs 82 }
253     }
254     if (!periodic[0]) {
255     for (i1=0;i1<N1;i1++) {
256 jgs 153 out->Nodes->Tag[M0*0+M1*i1]+=1;
257     out->Nodes->Tag[M0*(N0-1)+M1*i1]+=2;
258 jgs 82 }
259     }
260 jgs 153
261 jgs 82 /* set the elements: */
262    
263     #pragma omp parallel for private(i0,i1,k,node0)
264     for (i1=0;i1<NE1;i1++) {
265     for (i0=0;i0<NE0;i0++) {
266     k=i0+NE0*i1;
267     node0=i0+i1*N0;
268    
269     out->Elements->Id[k]=k;
270     out->Elements->Tag[k]=0;
271     out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);
272    
273     out->Elements->Nodes[INDEX2(0,k,4)]=node0;
274     out->Elements->Nodes[INDEX2(1,k,4)]=node0+1;
275     out->Elements->Nodes[INDEX2(2,k,4)]=node0+N0+1;
276     out->Elements->Nodes[INDEX2(3,k,4)]=node0+N0;
277    
278     }
279     }
280 jgs 123 out->Elements->minColor=0;
281     out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
282 jgs 82
283     /* face elements: */
284    
285     if (useElementsOnFace) {
286     NUMNODES=4;
287     } else {
288     NUMNODES=2;
289     }
290     totalNECount=NE0*NE1;
291     faceNECount=0;
292    
293     /* elements on boundary 010 (x2=0): */
294     if (!periodic[1]) {
295     #pragma omp parallel for private(i0,k,node0)
296     for (i0=0;i0<NE0;i0++) {
297     k=i0+faceNECount;
298     node0=i0;
299    
300     out->FaceElements->Id[k]=i0+totalNECount;
301     out->FaceElements->Tag[k]=10;
302     out->FaceElements->Color[k]=i0%2;
303    
304     if (useElementsOnFace) {
305     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
306     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
307     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
308     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
309     } else {
310     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
311     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
312     }
313     }
314     totalNECount+=NE0;
315     faceNECount+=NE0;
316    
317     /* ** elements on boundary 020 (x2=1): */
318    
319     #pragma omp parallel for private(i0,k,node0)
320     for (i0=0;i0<NE0;i0++) {
321     k=i0+faceNECount;
322     node0=i0+(NE1-1)*N0;
323    
324     out->FaceElements->Id[k]=i0+totalNECount;
325     out->FaceElements->Tag[k]=20;
326     out->FaceElements->Color[k]=i0%2+2;
327    
328     if (useElementsOnFace) {
329     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0+1;
330     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
331     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;
332     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
333     } else {
334     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0+1;
335     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
336     }
337     }
338     totalNECount+=NE0;
339     faceNECount+=NE0;
340     }
341     if (!periodic[0]) {
342     /* ** elements on boundary 001 (x1=0): */
343     #pragma omp parallel for private(i1,k,node0)
344     for (i1=0;i1<NE1;i1++) {
345     k=i1+faceNECount;
346     node0=i1*N0;
347    
348     out->FaceElements->Id[k]=i1+totalNECount;
349     out->FaceElements->Tag[k]=1;
350     out->FaceElements->Color[k]=i1%2+4;
351    
352     if (useElementsOnFace) {
353     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
354     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
355     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
356     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
357     } else {
358     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
359     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
360    
361     }
362     }
363     totalNECount+=NE1;
364     faceNECount+=NE1;
365    
366     /* ** elements on boundary 002 (x1=1): */
367    
368     #pragma omp parallel for private(i1,k,node0)
369     for (i1=0;i1<NE1;i1++) {
370     k=i1+faceNECount;
371     node0=(NE0-1)+i1*N0;
372    
373     out->FaceElements->Id[k]=i1+totalNECount;
374     out->FaceElements->Tag[k]=2;
375     out->FaceElements->Color[k]=i1%2+6;
376    
377     if (useElementsOnFace) {
378     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
379     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
380     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0;
381     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
382     } else {
383     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
384     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
385     }
386     }
387     totalNECount+=NE1;
388     faceNECount+=NE1;
389     }
390 jgs 123 out->FaceElements->minColor=0;
391     out->FaceElements->maxColor=7;
392 jgs 82
393 bcumming 730
394 jgs 82 /* face elements done: */
395    
396     /* condense the nodes: */
397    
398     Finley_Mesh_resolveNodeIds(out);
399    
400     /* prepare mesh for further calculatuions:*/
401    
402     Finley_Mesh_prepare(out) ;
403    
404 jgs 149 #ifdef Finley_TRACE
405 jgs 82 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
406 jgs 149 #endif
407 jgs 82
408 jgs 150 if (! Finley_noError()) {
409 jgs 82 Finley_Mesh_dealloc(out);
410     return NULL;
411     }
412     return out;
413     }
414 bcumming 730
415     /*----------------------------------------------------------------------------
416     MPI VERSION
417    
418    
419     ASSUMPTIONS
420     ===========
421    
422     - the domain dimension is large enough in the NE0 direction for each domain
423     to be 2 internal nodes wide in that direction. Face element calculation gets
424     buggered up otherwise.
425     - if dividing into two domains with periodic[0]=TRUE , then the global domain
426     must be at least 4 elements along the NE0 direction, otherwise redundant
427     nodes are generated.
428    
429     These problems can be avoided by ensuring that numElements[0]/mpiSize >= 2
430    
431     if domains are small enough to cause these problems, you probably don't
432     need to solve them in parallel. If you really do, rotate the domain so that the
433     long axis is aligned with NE0.
434     ------------------------------------------------------------------------------*/
435    
436    
437     #else
438     {
439     dim_t N0,N1,NE0,NE1,i0,i1, NE0_local, NDOF0,NDOF1, NFaceElements, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal, faceNECount, totalNECount,M0,M1;
440     index_t innerLoop=0;
441     index_t NUMNODES,k,firstNode=0, DOFcount=0, forwardDOF[2], backwardDOF[2], node0, node1;
442     index_t targetDomain=-1;
443     index_t *indexForward=NULL;
444     bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE;
445     Finley_Mesh* out;
446     char name[50];
447     double time0=Finley_timer();
448    
449     NE0=MAX(1,numElements[0]);
450     NE1=MAX(1,numElements[1]);
451     N0=NE0+1;
452     N1=NE1+1;
453     Paso_MPIInfo *mpi_info = NULL;
454    
455     /* get MPI information */
456     mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
457     if (! Finley_noError())
458     return NULL;
459    
460     if( mpi_info->rank==0 )
461     domLeft = TRUE;
462     if( mpi_info->rank==mpi_info->size-1 )
463     domRight = TRUE;
464     if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
465     domInternal = TRUE;
466    
467     /* dimensions of the local subdomain */
468     domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal );
469    
470     if (N0<=N1) {
471     M0=1;
472     M1=N0;
473     } else {
474     M0=N1;
475     M1=1;
476     }
477    
478     NFaceElements=0;
479     if (!periodic[0])
480     {
481     if( domLeft )
482     NFaceElements+=NE1;
483     if( domRight )
484     NFaceElements+=NE1;
485     }
486     if (!periodic[1]) {
487     NDOF1=N1;
488     NFaceElements+=2*NE0;
489     } else {
490     NDOF1=N1-1;
491     }
492    
493     /* allocate mesh: */
494    
495     sprintf(name,"Rectangular %d x %d mesh",N0,N1);
496     out=Finley_Mesh_alloc( name, 2, order, mpi_info );
497    
498     if (! Finley_noError()) return NULL;
499    
500     out->Elements=Finley_ElementFile_alloc( Rec4, out->order, mpi_info );
501     if (useElementsOnFace) {
502     out->FaceElements=Finley_ElementFile_alloc(Rec4Face,out->order, mpi_info );
503     out->ContactElements=Finley_ElementFile_alloc(Rec4Face_Contact,out->order, mpi_info );
504     } else {
505     out->FaceElements=Finley_ElementFile_alloc(Line2,out->order, mpi_info );
506     out->ContactElements=Finley_ElementFile_alloc(Line2_Contact,out->order, mpi_info );
507     }
508     out->Points=Finley_ElementFile_alloc(Point1,out->order, mpi_info );
509     if (! Finley_noError()) {
510     Finley_Mesh_dealloc(out);
511     return NULL;
512     }
513    
514     /* allocate tables: */
515     Finley_NodeFile_allocTable( out->Nodes, (numNodesLocal + numNodesExternal)*N1 );
516     Finley_NodeDistibution_allocTable( out->Nodes->degreeOfFreedomDistribution, numNodesLocal*NDOF1, numNodesExternal*NDOF1, 0 );
517     Finley_ElementFile_allocTable(out->Elements,numElementsLocal*NE1);
518     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
519    
520     if (! Finley_noError()) {
521     Finley_Mesh_dealloc(out);
522     return NULL;
523     }
524     /* set nodes: */
525    
526     //#pragma omp parallel for private(i0,i1,k)
527     if( mpi_info->size==1 )
528     innerLoop = numNodesLocal;
529     else
530     innerLoop = numNodesLocal-periodicLocal[0];
531    
532     for (k=0,i1=0;i1<N1;i1++) {
533     for ( i0=0;i0<innerLoop;i0++,k++) {
534     out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(i0+firstNode)/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]=i0+innerLoop*i1;
537     out->Nodes->Tag[k]=0;
538     out->Nodes->degreeOfFreedom[k]=i0%numDOFLocal + numDOFLocal*(i1%NDOF1);
539     }
540     }
541    
542     /* tag 'em */
543     for( i0=0; i0<numNodesLocal; i0++ )
544     {
545     out->Nodes->Tag[i0] += 10; // bottom
546     out->Nodes->Tag[i0+numNodesLocal*(N1-1)] += 20; // top
547     }
548     if( domLeft && !periodicLocal[0] )
549     for( i1=0; i1<N1; i1++ )
550     out->Nodes->Tag[i1*numNodesLocal] += 1; // left
551    
552     if( domRight && !periodicLocal[0] )
553     for( i1=0; i1<N1; i1++ )
554     out->Nodes->Tag[(i1+1)*numNodesLocal-1] += 2; // right
555    
556     /* external nodes */
557     k = innerLoop*N1;
558     DOFcount = out->Nodes->degreeOfFreedom[k-1]+1;
559     if( mpi_info->size>1 )
560     {
561     /* add forward communication information for boundary nodes */
562     indexForward = TMPMEMALLOC( NDOF1, index_t );
563     if( domInternal || domRight || periodicLocal[0] )
564     {
565     for( int i=0; i<NDOF1; i++ )
566     indexForward[i] = i*numDOFLocal;
567     targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;
568     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );
569     }
570     if( domInternal || domLeft || periodicLocal[1] )
571     {
572     for( int i=0; i<NDOF1; i++ )
573     indexForward[i] = (i+1)*numDOFLocal - 1;
574     targetDomain = (mpi_info->rank+1) % mpi_info->size;
575     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );
576     }
577     TMPMEMFREE( indexForward );
578    
579     /* LHS */
580     if( periodicLocal[0] )
581     {
582     for (i1=0; i1<N1; i1++, k++)
583     {
584     out->Nodes->Coordinates[INDEX2(0,k,2)]=Length[0];
585     out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
586     out->Nodes->Id[k]=k;
587     out->Nodes->Tag[k]=0;
588     out->Nodes->degreeOfFreedom[k] = (i1 % NDOF1)*numDOFLocal;
589     }
590     out->Nodes->Tag[k-N1] = 10;
591     out->Nodes->Tag[k-1] = 20;
592     for (i1=0; i1<N1; i1++, k++)
593     {
594     out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(NE0-1)/DBLE(NE0)*Length[0];
595     out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
596     out->Nodes->Id[k]=k;
597     out->Nodes->Tag[k]=0;
598     out->Nodes->degreeOfFreedom[k] = DOFcount + i1%NDOF1;
599     }
600     out->Nodes->Tag[k-N1] = 10;
601     out->Nodes->Tag[k-1] = 20;
602     DOFcount = out->Nodes->degreeOfFreedom[k-1]+1;
603    
604     targetDomain = mpi_info->size-1;
605     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, out->Nodes->degreeOfFreedom+(numNodesLocal*N1) );
606     }
607     if( mpi_info->rank>0 )
608     {
609     for (i1=0; i1<N1; i1++, k++)
610     {
611     out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(firstNode-1)/DBLE(N0-1)*Length[0];
612     out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
613     out->Nodes->Id[k]=k;
614     out->Nodes->Tag[k]=0;
615     out->Nodes->degreeOfFreedom[k] = DOFcount + i1%NDOF1;
616     }
617     out->Nodes->Tag[k-N1] = 10;
618     out->Nodes->Tag[k-1] = 20;
619     DOFcount = out->Nodes->degreeOfFreedom[k-1]+1;
620    
621     targetDomain = mpi_info->rank-1;
622     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, out->Nodes->degreeOfFreedom+(k-N1) );
623     }
624    
625     /* RHS */
626     if( periodicLocal[1] || (mpi_info->rank < mpi_info->size-1) )
627     {
628     for (i1=0; i1<N1; i1++, k++)
629     {
630     out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(firstNode+numNodesLocal-periodicLocal[0])/DBLE(N0-1)*Length[0];
631     out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
632     out->Nodes->Id[k]=k;
633     out->Nodes->Tag[k]=0;
634     out->Nodes->degreeOfFreedom[k] = DOFcount + i1%NDOF1;
635     }
636     out->Nodes->Tag[k-N1] = 10;
637     out->Nodes->Tag[k-1] = 20;
638     DOFcount = out->Nodes->degreeOfFreedom[k-1]+1;
639    
640     targetDomain = (mpi_info->rank+1) % mpi_info->size;
641     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, out->Nodes->degreeOfFreedom+(k-N1) );
642     }
643     }
644    
645     /* set the elements: */
646    
647     /* internal */
648     //#pragma omp parallel for private(i0,i1,k,node0)
649     for ( i1=0, k=0; i1<NE1; i1++)
650     {
651     for ( i0=0; i0<numElementsInternal; i0++, k++)
652     {
653     node0=i0+i1*innerLoop;
654    
655     out->Elements->Id[k]=k;
656     out->Elements->Tag[k]=0;
657     out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
658    
659     out->Elements->Nodes[INDEX2(0,k,4)]=node0;
660     out->Elements->Nodes[INDEX2(1,k,4)]=node0+1;
661     out->Elements->Nodes[INDEX2(2,k,4)]=node0+innerLoop+1;
662     out->Elements->Nodes[INDEX2(3,k,4)]=node0+innerLoop;
663    
664     }
665     }
666     /* boundary */
667     if( mpi_info->size>1 )
668     {
669     if( domInternal )
670     {
671     /* left */
672     for( i1=0; i1<NE1; i1++, k++ )
673     {
674     node1=numNodesLocal*N1 + i1;
675     node0=i1*numNodesLocal;
676    
677     out->Elements->Id[k]=k;
678     out->Elements->Tag[k]=0;
679     out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
680    
681     out->Elements->Nodes[INDEX2(0,k,4)]=node1;
682     out->Elements->Nodes[INDEX2(1,k,4)]=node0;
683     out->Elements->Nodes[INDEX2(2,k,4)]=node0+numNodesLocal;
684     out->Elements->Nodes[INDEX2(3,k,4)]=node1+1;
685     }
686    
687     /* right */
688     for( i1=0; i1<NE1; i1++, k++ )
689     {
690     node0 = (i1+1)*numNodesLocal-1;
691     node1 = (numNodesLocal+1)*N1 + i1;
692    
693     out->Elements->Id[k]=k;
694     out->Elements->Tag[k]=0;
695     out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
696    
697     out->Elements->Nodes[INDEX2(0,k,4)]=node0;
698     out->Elements->Nodes[INDEX2(1,k,4)]=node1;
699     out->Elements->Nodes[INDEX2(2,k,4)]=node1+1;
700     out->Elements->Nodes[INDEX2(3,k,4)]=node0+numNodesLocal;
701     }
702     }
703     else if( domLeft )
704     {
705     /* left */
706     if( periodicLocal[0] )
707     {
708     node1 = numNodesLocal*N1;
709     node0 = innerLoop*N1;
710     for( i1=0; i1<NE1; i1++, k++, node1++, node0++ )
711     {
712     out->Elements->Id[k]=k;
713     out->Elements->Tag[k]=0;
714     out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
715    
716     out->Elements->Nodes[INDEX2(0,k,4)]=node1;
717     out->Elements->Nodes[INDEX2(1,k,4)]=node0;
718     out->Elements->Nodes[INDEX2(2,k,4)]=node0+1;
719     out->Elements->Nodes[INDEX2(3,k,4)]=node1+1;
720     }
721     }
722     /* right */
723     for( i1=0; i1<NE1; i1++, k++ )
724     {
725     node0 = (i1+1)*innerLoop-1;
726     node1 = (numNodesLocal+periodicLocal[0])*N1 + i1;
727    
728     out->Elements->Id[k]=k;
729     out->Elements->Tag[k]=0;
730     out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
731    
732     out->Elements->Nodes[INDEX2(0,k,4)]=node0;
733     out->Elements->Nodes[INDEX2(1,k,4)]=node1;
734     out->Elements->Nodes[INDEX2(2,k,4)]=node1+1;
735     out->Elements->Nodes[INDEX2(3,k,4)]=node0+numNodesLocal-periodicLocal[0];
736     }
737     }
738     else if( domRight )
739     {
740     /* left */
741     for( i1=0; i1<NE1; i1++, k++ )
742     {
743     node1=numNodesLocal*N1 + i1;
744     node0=i1*numNodesLocal;
745    
746     out->Elements->Id[k]=k;
747     out->Elements->Tag[k]=0;
748     out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
749    
750     out->Elements->Nodes[INDEX2(0,k,4)]=node1;
751     out->Elements->Nodes[INDEX2(1,k,4)]=node0;
752     out->Elements->Nodes[INDEX2(2,k,4)]=node0+numNodesLocal;
753     out->Elements->Nodes[INDEX2(3,k,4)]=node1+1;
754     }
755     /* right */
756     if( periodicLocal[1] )
757     {
758     for( i1=0; i1<NE1; i1++, k++ )
759     {
760     node0=(i1+1)*numNodesLocal-1;
761     node1 = (numNodesLocal+1)*N1 + i1;
762    
763     out->Elements->Id[k]=k;
764     out->Elements->Tag[k]=0;
765     out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
766    
767     out->Elements->Nodes[INDEX2(0,k,4)]=node0;
768     out->Elements->Nodes[INDEX2(1,k,4)]=node1;
769     out->Elements->Nodes[INDEX2(2,k,4)]=node1+1;
770     out->Elements->Nodes[INDEX2(3,k,4)]=node0+numNodesLocal;
771     }
772     }
773     }
774     }
775     out->Elements->minColor=0;
776     out->Elements->maxColor=0; //COLOR_MOD(0)+3*COLOR_MOD(0);
777    
778     if( domInternal )
779     {
780     out->Elements->elementDistribution->numInternal = NE1*(numElementsLocal-2);
781     out->Elements->elementDistribution->numBoundary = NE1*2;
782     }
783     else if( mpi_info->size==1 )
784     {
785     out->Elements->elementDistribution->numInternal = NE1*numElementsLocal;
786     }
787     else
788     {
789     out->Elements->elementDistribution->numInternal += NE1*(numElementsLocal-1-periodic[0]);
790     out->Elements->elementDistribution->numBoundary += NE1*(1 + periodic[0]);
791     }
792    
793     out->FaceElements->elementDistribution->numLocal = out->FaceElements->elementDistribution->numInternal + out->FaceElements->elementDistribution->numBoundary;
794    
795     /* face elements: */
796    
797     if (useElementsOnFace) {
798     NUMNODES=4;
799     } else {
800     NUMNODES=2;
801     }
802     totalNECount=numElementsLocal*NE1;
803     faceNECount=0;
804    
805     /* do the internal elements first */
806     if (!periodic[1])
807     {
808     /* elements on boundary 010 (x1=0): */
809     #pragma omp parallel for private(i0,k,node0)
810     for (i0=0;i0<numElementsLocal-numNodesExternal;i0++) {
811     k=i0+faceNECount;
812     node0=i0;
813    
814     out->FaceElements->Id[k]=i0+totalNECount;
815     out->FaceElements->Tag[k] = 10;
816     out->FaceElements->Color[k] = 0; // i0%2;
817    
818     if (useElementsOnFace) {
819     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
820     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
821     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal+1;
822     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal;
823     } else {
824     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
825     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
826     }
827     }
828     totalNECount += numElementsLocal-numNodesExternal;
829     faceNECount += numElementsLocal-numNodesExternal;
830    
831     /* ** elements on boundary 020 (x1=1): */
832    
833     #pragma omp parallel for private(i0,k,node0)
834     for (i0=0;i0<numElementsLocal-numNodesExternal;i0++) {
835     k=i0+faceNECount;
836     node0=i0+(NE1-1)*numNodesLocal;
837    
838     out->FaceElements->Id[k]=i0+totalNECount;
839     out->FaceElements->Tag[k] = 20;
840     out->FaceElements->Color[k] = 0; // i0%2+2;
841    
842     if (useElementsOnFace) {
843     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numNodesLocal+1;
844     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal;
845     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;
846     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
847     } else {
848     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numNodesLocal+1;
849     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal;
850     }
851     }
852     totalNECount += numElementsLocal-numNodesExternal;
853     faceNECount += numElementsLocal-numNodesExternal;
854     }
855     if ( !periodic[0])
856     {
857     /* ** elements on boundary 001 (x0=0): */
858     if( domLeft )
859     {
860     #pragma omp parallel for private(i1,k,node0)
861     for (i1=0;i1<NE1;i1++) {
862     k=i1+faceNECount;
863     node0=i1*numNodesLocal;
864    
865     out->FaceElements->Id[k]=i1+totalNECount;
866     out->FaceElements->Tag[k] = 1;
867     out->FaceElements->Color[k] = 0; //i1%2+4;
868    
869     if (useElementsOnFace) {
870     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numNodesLocal;
871     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
872     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
873     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal+1;
874     } else {
875     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numNodesLocal;
876     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
877    
878     }
879     }
880     totalNECount+=NE1;
881     faceNECount+=NE1;
882     }
883     /* ** elements on boundary 002 (x0=1): */
884     if( domRight )
885     {
886     #pragma omp parallel for private(i1,k,node0)
887     for (i1=0;i1<NE1;i1++) {
888     k=i1+faceNECount;
889     node0=(numNodesLocal-2)+i1*numNodesLocal;
890    
891     out->FaceElements->Id[k]=i1+totalNECount;
892     out->FaceElements->Tag[k] = 2;
893     out->FaceElements->Color[k] = 0; //i1%2+6;
894    
895     if (useElementsOnFace) {
896     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
897     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal+1;
898     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal;
899     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
900     } else {
901     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
902     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal+1;
903     }
904     }
905     totalNECount+=NE1;
906     faceNECount+=NE1;
907     }
908     }
909    
910     if( mpi_info->size>1 )
911     {
912     if( !periodic[1] && (domInternal || domRight) )
913     {
914     // bottom left
915     k = faceNECount;
916     node0 = numNodesLocal*N1;
917    
918     out->FaceElements->Id[k]=totalNECount;
919     out->FaceElements->Tag[k] = 10;
920     out->FaceElements->Color[k] = 0; //i1%2+6;
921    
922     if (useElementsOnFace) {
923     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
924     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=0;
925     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=numNodesLocal;
926     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
927     } else {
928     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
929     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=0;
930     }
931     totalNECount++;
932     faceNECount++;
933    
934     // top left
935     k = faceNECount;
936     node0 = (numNodesLocal+1)*N1 - 2;
937    
938     out->FaceElements->Id[k]=totalNECount;
939     out->FaceElements->Tag[k] = 20;
940     out->FaceElements->Color[k] = 0; //i1%2+6;
941    
942     if (useElementsOnFace) {
943     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=numNodesLocal*(N1-1);
944     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
945     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;
946     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=numNodesLocal*(N1-2);
947     } else {
948     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=numNodesLocal*(N1-1);
949     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
950     }
951     totalNECount++;
952     faceNECount++;
953     }
954     if( !periodic[1] && (domInternal || domLeft) )
955     {
956     // bottom right
957     k = faceNECount;
958     node0 = (numNodesLocal+nodesExternal[0])*N1;
959    
960     out->FaceElements->Id[k]=totalNECount;
961     out->FaceElements->Tag[k] = 10;
962     out->FaceElements->Color[k] = 0; //i1%2+6;
963    
964     if (useElementsOnFace) {
965     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=numNodesLocal-1;
966     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
967     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
968     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=2*numNodesLocal-1;
969     } else {
970     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=numNodesLocal-1;
971     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
972     }
973     totalNECount++;
974     faceNECount++;
975    
976     // top right
977     k = faceNECount;
978     node0 = (numNodesLocal+1+nodesExternal[0])*N1-2;
979    
980     out->FaceElements->Id[k]=totalNECount;
981     out->FaceElements->Tag[k] = 20;
982     out->FaceElements->Color[k] = 0; //i1%2+6;
983    
984     if (useElementsOnFace) {
985     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
986     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=numNodesLocal*N1-1;
987     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=numNodesLocal*(N1-1)-1;
988     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
989     } else {
990     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
991     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=numNodesLocal*N1-1;
992     }
993     totalNECount++;
994     faceNECount++;
995     }
996     }
997     out->FaceElements->minColor = 0;
998     out->FaceElements->maxColor = 0; //7;
999    
1000     if( domInternal )
1001     {
1002     if( !periodic[1] )
1003     {
1004     out->FaceElements->elementDistribution->numInternal = 2*(numElementsLocal-2);
1005     out->FaceElements->elementDistribution->numBoundary = 4;
1006     }
1007     }
1008     else if( mpi_info->size==1 )
1009     {
1010     if( !periodic[0] )
1011     out->FaceElements->elementDistribution->numInternal += NE1*2;
1012     if( !periodic[1] )
1013     out->FaceElements->elementDistribution->numInternal += numElementsLocal*2;
1014     }
1015     else
1016     {
1017     if( !periodic[0] )
1018     out->FaceElements->elementDistribution->numInternal += NE1;
1019     if( !periodic[1] )
1020     {
1021     out->FaceElements->elementDistribution->numInternal += 2*(numElementsLocal-1-periodic[0]);
1022     out->FaceElements->elementDistribution->numBoundary += 2*(1 + periodic[0]);
1023     }
1024     }
1025    
1026     out->FaceElements->elementDistribution->numLocal = out->FaceElements->elementDistribution->numInternal + out->FaceElements->elementDistribution->numBoundary;
1027    
1028     if( out->FaceElements->elementDistribution->numLocal!=faceNECount )
1029     printf( "ERROR : face element numbers generated %d + %d = %d != %d\n",out->FaceElements->elementDistribution->numInternal, out->FaceElements->elementDistribution->numBoundary, out->FaceElements->elementDistribution->numLocal, faceNECount );
1030    
1031    
1032    
1033     /* face elements done: */
1034    
1035     /* condense the nodes: */
1036     Finley_Mesh_resolveNodeIds( out );
1037    
1038     /* prepare mesh for further calculatuions:*/
1039    
1040     //Finley_Mesh_prepare(out) ;
1041    
1042     #ifdef Finley_TRACE
1043     printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
1044     #endif
1045    
1046     if ( !Finley_MPI_noError( mpi_info )) {
1047     Finley_Mesh_dealloc(out);
1048     return NULL;
1049     }
1050     Paso_MPIInfo_dealloc( mpi_info );
1051    
1052     return out;
1053     }
1054     #endif
1055    
1056    
1057    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26