# Diff of /trunk/finley/src/Mesh_line3.c

revision 781 by bcumming, Mon Jun 26 01:46:34 2006 UTC revision 782 by bcumming, Tue Jul 18 00:47:47 2006 UTC
# Line 26  Line 26
26
27  #include "RectangularMesh.h"  #include "RectangularMesh.h"
28

/**************************************************************
The code for Mesh_line[2/3].c is messy, it was used as a test
bed for getting the MPI distributed mesh figured out before
doing the rectangle/brick domains. The code is efficient, it just
uses too many if/else if/else statements. Unfortunately, calculating
a distributed mesh is an easy thing to visualise, but a devil to implement!

Sorry if I don't get around
to neatening this code up, but rest assured, the more important
(and complex) rectangle and brick codes are much better organised.
**************************************************************/

29  #ifdef PASO_MPI  #ifdef PASO_MPI
30  /* get the number of nodes/elements for domain with rank=rank, of size processors  /* get the number of nodes/elements for domain with rank=rank, of size processors
31     where n is the total number of nodes/elements in the global domain */     where n is the total number of nodes/elements in the global domain */
# Line 134  static void domain_calculateDimension( i Line 120  static void domain_calculateDimension( i
120    {    {
121      DOFBoundary[0] = DOFBoundary[1] = 0;      DOFBoundary[0] = DOFBoundary[1] = 0;
122    }    }
/* some debugging printf statements */
//printf( "rank/size = %d/%d\nNodes : %d Local, %d External[%d %d], First = %d\nElements : %d Local\nDOF : %d Local, External [%d %d], Boundary [%d %d]\nperiodicLocal [%d %d]\n\n", rank, size, *numNodesLocal, *numNodesExternal, nodesExternal[0], nodesExternal[1], *firstNode, *numElementsLocal, *numDOFLocal, DOFExternal[0], DOFExternal[1], DOFBoundary[0], DOFBoundary[1], periodicLocal[0], periodicLocal[1] );
123  }  }
124  #endif  #endif
125
126
127    #ifdef PASO_MPI
128    Finley_Mesh* Finley_RectangularMesh_Line3_singleCPU(dim_t* numElements,double* Length,bool_t* periodic,index_t order,bool_t useElementsOnFace,Paso_MPIInfo *mpi_info)
129    #else
130  Finley_Mesh* Finley_RectangularMesh_Line3(dim_t* numElements,double* Length,bool_t* periodic,index_t order,bool_t useElementsOnFace)  Finley_Mesh* Finley_RectangularMesh_Line3(dim_t* numElements,double* Length,bool_t* periodic,index_t order,bool_t useElementsOnFace)
131  #ifndef PASO_MPI  #endif
132  {  {
133    dim_t N0,NE0,i0,NDOF0,NFaceElements,NUMNODES;    dim_t N0,NE0,i0,NDOF0,NFaceElements,NUMNODES;
134    index_t k,node0;    index_t k,node0;
# Line 161  Finley_Mesh* Finley_RectangularMesh_Line Line 148  Finley_Mesh* Finley_RectangularMesh_Line
148    /*  allocate mesh: */    /*  allocate mesh: */
149
150    sprintf(name,"Rectangular mesh with %d nodes",N0);    sprintf(name,"Rectangular mesh with %d nodes",N0);
151    /* TEMPFIX */  #ifdef PASO_MPI
152      out=Finley_Mesh_alloc(name,1,order,mpi_info);
153      if (! Finley_noError()) return NULL;
154
155      out->Elements=Finley_ElementFile_alloc(Line3,out->order,mpi_info);
156      if (useElementsOnFace) {
157        out->FaceElements=Finley_ElementFile_alloc(Line3Face,out->order,mpi_info);
158        out->ContactElements=Finley_ElementFile_alloc(Line3Face_Contact,out->order,mpi_info);
159      } else {
160        out->FaceElements=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
161        out->ContactElements=Finley_ElementFile_alloc(Point1_Contact,out->order,mpi_info);
162      }
163      out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
164      if (! Finley_noError()) {
165           Finley_Mesh_dealloc(out);
166           return NULL;
167      }
168
169      /*  allocate tables: */
170
171      Finley_NodeFile_allocTable(out->Nodes,N0);
172      Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, N0, 0, 0);
173      Finley_ElementFile_allocTable(out->Elements,NE0);
174      Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
175    #else
176    out=Finley_Mesh_alloc(name,1,order);    out=Finley_Mesh_alloc(name,1,order);
177    if (! Finley_noError()) return NULL;    if (! Finley_noError()) return NULL;
178
# Line 186  Finley_Mesh* Finley_RectangularMesh_Line Line 196  Finley_Mesh* Finley_RectangularMesh_Line
196
197    Finley_ElementFile_allocTable(out->Elements,NE0);    Finley_ElementFile_allocTable(out->Elements,NE0);
198    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
199    #endif
200    if (! Finley_noError()) {    if (! Finley_noError()) {
201        Finley_Mesh_dealloc(out);        Finley_Mesh_dealloc(out);
202        return NULL;        return NULL;
# Line 200  Finley_Mesh* Finley_RectangularMesh_Line Line 211  Finley_Mesh* Finley_RectangularMesh_Line
211       out->Nodes->Id[k]=k;       out->Nodes->Id[k]=k;
212       out->Nodes->Tag[k]=0;       out->Nodes->Tag[k]=0;
213       out->Nodes->degreeOfFreedom[k]=(i0%NDOF0);       out->Nodes->degreeOfFreedom[k]=(i0%NDOF0);
214    #ifdef PASO_MPI
215             out->Nodes->Dom[k] = NODE_INTERNAL;
216    #endif
217    }    }
218    if (!periodic[0]) {    if (!periodic[0]) {
219       out->Nodes->Tag[0]=1;       out->Nodes->Tag[0]=1;
# Line 216  Finley_Mesh* Finley_RectangularMesh_Line Line 230  Finley_Mesh* Finley_RectangularMesh_Line
230      out->Elements->Id[k]=k;      out->Elements->Id[k]=k;
231      out->Elements->Tag[k]=0;      out->Elements->Tag[k]=0;
232      out->Elements->Color[k]=COLOR_MOD(i0);      out->Elements->Color[k]=COLOR_MOD(i0);
233    #ifdef PASO_MPI
234             out->Elements->Dom[k] = ELEMENT_INTERNAL;
235    #endif
236
237      out->Elements->Nodes[INDEX2(0,k,3)]=node0;      out->Elements->Nodes[INDEX2(0,k,3)]=node0;
238      out->Elements->Nodes[INDEX2(1,k,3)]=node0+2;      out->Elements->Nodes[INDEX2(1,k,3)]=node0+2;
# Line 235  Finley_Mesh* Finley_RectangularMesh_Line Line 252  Finley_Mesh* Finley_RectangularMesh_Line
252       out->FaceElements->Id[0]=NE0;       out->FaceElements->Id[0]=NE0;
253       out->FaceElements->Tag[0]=1;       out->FaceElements->Tag[0]=1;
254       out->FaceElements->Color[0]=0;       out->FaceElements->Color[0]=0;
255    #ifdef PASO_MPI
256             out->FaceElements->Dom[0] = ELEMENT_INTERNAL;
257    #endif
258       if (useElementsOnFace) {       if (useElementsOnFace) {
259         out->FaceElements->Nodes[INDEX2(0,0,NUMNODES)]=0;         out->FaceElements->Nodes[INDEX2(0,0,NUMNODES)]=0;
260         out->FaceElements->Nodes[INDEX2(1,0,NUMNODES)]=2;         out->FaceElements->Nodes[INDEX2(1,0,NUMNODES)]=2;
# Line 246  Finley_Mesh* Finley_RectangularMesh_Line Line 266  Finley_Mesh* Finley_RectangularMesh_Line
266       out->FaceElements->Id[1]=NE0+1;       out->FaceElements->Id[1]=NE0+1;
267       out->FaceElements->Tag[1]=2;       out->FaceElements->Tag[1]=2;
268       out->FaceElements->Color[1]=1;       out->FaceElements->Color[1]=1;
269    #ifdef PASO_MPI
270             out->FaceElements->Dom[1] = ELEMENT_INTERNAL;
271    #endif
272       if (useElementsOnFace) {       if (useElementsOnFace) {
273          out->FaceElements->Nodes[INDEX2(0,1,NUMNODES)]=N0-1;          out->FaceElements->Nodes[INDEX2(0,1,NUMNODES)]=N0-1;
274          out->FaceElements->Nodes[INDEX2(1,1,NUMNODES)]=N0-3;          out->FaceElements->Nodes[INDEX2(1,1,NUMNODES)]=N0-3;
# Line 257  Finley_Mesh* Finley_RectangularMesh_Line Line 280  Finley_Mesh* Finley_RectangularMesh_Line
280    out->FaceElements->minColor=0;    out->FaceElements->minColor=0;
281    out->FaceElements->maxColor=1;    out->FaceElements->maxColor=1;
282
283    /*  face elements done: */  #ifdef PASO_MPI
284          Finley_ElementFile_setDomainFlags( out->Elements );
285        Finley_ElementFile_setDomainFlags( out->FaceElements );
286        Finley_ElementFile_setDomainFlags( out->ContactElements );
287        Finley_ElementFile_setDomainFlags( out->Points );
288
289        /* reorder the degrees of freedom */
290        Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
291    #endif
292
293    /*   condense the nodes: */    /*   condense the nodes: */

294    Finley_Mesh_resolveNodeIds(out);    Finley_Mesh_resolveNodeIds(out);
295
296    /* prepare mesh for further calculations:*/    /* prepare mesh for further calculations:*/

297    Finley_Mesh_prepare(out) ;    Finley_Mesh_prepare(out) ;
298
299    #ifdef Finley_TRACE    #ifdef Finley_TRACE
# Line 277  Finley_Mesh* Finley_RectangularMesh_Line Line 306  Finley_Mesh* Finley_RectangularMesh_Line
306    }    }
307    return out;    return out;
308  }  }
#else
309  /* MPI version */  /* MPI version */
310    #ifdef PASO_MPI
311    Finley_Mesh* Finley_RectangularMesh_Line3(dim_t* numElements,double* Length,bool_t* periodic,index_t order,bool_t useElementsOnFace)
312  {  {
313    dim_t N0, NE0, NE0_local, i0, NDOF0, NFaceElements, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal, DOFBoundary[2];    dim_t N0, NE0, NE0_local, i0, NDOF0, NFaceElements, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal, DOFBoundary[2], N0t, NDOF0t;
314    index_t *numForward=NULL, *numBackward=NULL;    index_t *numForward=NULL, *numBackward=NULL;
315    index_t NUMNODES, node0,k, i,firstNode=0, DOFcount=0, forwardDOF[4], backwardDOF[4];    index_t NUMNODES, node0,k, i,firstNode=0, DOFcount=0, forwardDOF[4], backwardDOF[4], firstNodeConstruct;
316    index_t targetDomain=-1;    index_t targetDomain=-1;
317    bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE;    bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft, boundaryRight;
318    Finley_Mesh* out=NULL;    Finley_Mesh* out=NULL;
319    char name[50];    char name[50];
320    double time0=Finley_timer();    double time0=Finley_timer();
# Line 303  Finley_Mesh* Finley_RectangularMesh_Line Line 333  Finley_Mesh* Finley_RectangularMesh_Line
333          Finley_Mesh_dealloc(out);          Finley_Mesh_dealloc(out);
334          return NULL;          return NULL;
335    }    }
336
337        /* use the serial code to generate the mesh in the 1-CPU case */
338        if( mpi_info->size==1 ){
339            out = Finley_RectangularMesh_Line3_singleCPU(numElements,Length,periodic,order,useElementsOnFace,mpi_info);
340            return out;
341        }
342
343    if( mpi_info->rank==0 )    if( mpi_info->rank==0 )
344      domLeft = TRUE;      domLeft = TRUE;
# Line 324  Finley_Mesh* Finley_RectangularMesh_Line Line 360  Finley_Mesh* Finley_RectangularMesh_Line
360        NFaceElements++;        NFaceElements++;
361    }      }
362
363        boundaryLeft = !domLeft || periodicLocal[0];
364        boundaryRight = !domRight || periodicLocal[1];
365        N0t = numNodesLocal + boundaryRight + 2*boundaryLeft;
366        NDOF0t = numDOFLocal + boundaryRight + 2*boundaryLeft;
367        firstNodeConstruct = firstNode - 2*boundaryLeft;
368        firstNodeConstruct = firstNodeConstruct<0 ? N0-3 : firstNodeConstruct;
369
370    /*  allocate mesh: */    /*  allocate mesh: */
371    sprintf(name,"Rectangular mesh with %d nodes",N0);    sprintf(name,"Rectangular mesh with %d nodes",N0);
372    out=Finley_Mesh_alloc(name,1,order,mpi_info);    out=Finley_Mesh_alloc(name,1,order,mpi_info);
# Line 345  Finley_Mesh* Finley_RectangularMesh_Line Line 388  Finley_Mesh* Finley_RectangularMesh_Line
388
389    /*  allocate tables: */    /*  allocate tables: */
390    Finley_NodeFile_allocTable(out->Nodes,numNodesLocal+numNodesExternal);    Finley_NodeFile_allocTable(out->Nodes,numNodesLocal+numNodesExternal);
391    Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal, DOFExternal[0]+DOFExternal[1], 0 );    Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal, 3, 0 );
392    Finley_ElementFile_allocTable(out->Elements,numElementsLocal);    Finley_ElementFile_allocTable(out->Elements,numElementsLocal);
393    if( NFaceElements )    if( NFaceElements )
394      Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);      Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
# Line 356  Finley_Mesh* Finley_RectangularMesh_Line Line 399  Finley_Mesh* Finley_RectangularMesh_Line
399
400    /*  set nodes: */    /*  set nodes: */
401    #pragma omp parallel for private(i0,k)      #pragma omp parallel for private(i0,k)
402    /* local nodes */    for (i0=0;i0<N0t;i0++) {
for (i0=0;i0<numNodesLocal;i0++) {
403       k=i0;       k=i0;
404       out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(i0+firstNode)/DBLE(N0-1)*Length[0];       out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE((i0+firstNodeConstruct) % N0)/DBLE(N0-1)*Length[0];
405       out->Nodes->Id[k]=k;       out->Nodes->Id[k]=k;
406       out->Nodes->Tag[k]=0;       out->Nodes->Tag[k]=0;
407       out->Nodes->degreeOfFreedom[k]=i0%(numDOFLocal);       out->Nodes->degreeOfFreedom[k]=k;
408             out->Nodes->Dom[k]=NODE_INTERNAL;
409    }    }
410
411    /* external nodes - do left then right hand side */    /* setup boundary DOF data */
412    /* the following only applies if more than one domain */      if( boundaryLeft ){
413            out->Nodes->Dom[0]= out->Nodes->Dom[1]  = NODE_EXTERNAL;
414    /* this is messy, could be done cleaner */          out->Nodes->Dom[2] = NODE_BOUNDARY;
415    if( mpi_info->size>1 )      }
416    {      else
417      DOFcount = numNodesLocal;          out->Nodes->Tag[0] += 1;
418      k=numNodesLocal;      if( boundaryRight ){
419      if( mpi_info->rank!=0 || periodicLocal[0] )          out->Nodes->Dom[N0t-1] = NODE_EXTERNAL;
420      {          out->Nodes->Dom[N0t-2] = out->Nodes->Dom[N0t-3] = NODE_BOUNDARY;
421        /* left hand boundary is periodic -      }
422           add the nodes/DOF that define the element on the right hand boundary */      else
423        if( periodicLocal[0] )          out->Nodes->Tag[N0t-1] += 2;
424        {      if( periodicLocal[0] ){
425          k--;          out->Nodes->degreeOfFreedom[3] = out->Nodes->degreeOfFreedom[2];
426          out->Nodes->Coordinates[INDEX2(0,k,1)]=Length[0];          out->Nodes->Dom[3] = NODE_BOUNDARY;
427          out->Nodes->Id[k]=k;      }
428          out->Nodes->Tag[k]=0;
429          out->Nodes->degreeOfFreedom[k]=0;      if( !(mpi_info->size==2 && periodicLocal[0])){
430          DOFcount--;          if( boundaryLeft  ) {
431          k++;              targetDomain = mpi_info->rank-1 < 0 ? mpi_info->size-1 : mpi_info->rank-1;
432          out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(N0-2)/DBLE(N0-1)*Length[0];              forwardDOF[0] = out->Nodes->degreeOfFreedom[2];
433          out->Nodes->Id[k]=k;              backwardDOF[0] = out->Nodes->degreeOfFreedom[0];
434          out->Nodes->Tag[k]=0;              backwardDOF[1] = out->Nodes->degreeOfFreedom[1];
435          out->Nodes->degreeOfFreedom[k]=DOFcount++;              Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, forwardDOF );
436          k++;              Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, backwardDOF );
437          out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(N0-3)/DBLE(N0-1)*Length[0];          }
438          out->Nodes->Id[k]=k;          if( boundaryRight ) {
439          out->Nodes->Tag[k]=0;              targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;
440          out->Nodes->degreeOfFreedom[k]=DOFcount++;              forwardDOF[0] = out->Nodes->degreeOfFreedom[N0t-3];
441          k++;              forwardDOF[1] = out->Nodes->degreeOfFreedom[N0t-2];
442        }              backwardDOF[0] = out->Nodes->degreeOfFreedom[N0t-1];
443        /* left hand boundary with another subdomain, need to add the nodes/DOFs that              Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, forwardDOF );
444           defines the element that spans the boundary */              Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, backwardDOF );
445        else          }
446        {    } else{
447          out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(firstNode-1)/DBLE(N0-1)*Length[0];          /* periodic boundary conditions with 2 domains, need to change the order in which domain 0 passes boundary data */
448          out->Nodes->Id[k]=k;          targetDomain = 1;
449          out->Nodes->Tag[k]=0;
450          out->Nodes->degreeOfFreedom[k]=DOFcount++;          forwardDOF[0] = out->Nodes->degreeOfFreedom[N0t-3];
451          k++;          forwardDOF[1] = out->Nodes->degreeOfFreedom[N0t-2];
452          out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(firstNode-2)/DBLE(N0-1)*Length[0];          backwardDOF[0] = out->Nodes->degreeOfFreedom[N0t-1];
453          out->Nodes->Id[k]=k;          Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, forwardDOF );
454          out->Nodes->Tag[k]=0;          Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, backwardDOF );
455          out->Nodes->degreeOfFreedom[k]=DOFcount++;
456          k++;          forwardDOF[0] = out->Nodes->degreeOfFreedom[2];
457        }          backwardDOF[0] = out->Nodes->degreeOfFreedom[0];
458      }          backwardDOF[1] = out->Nodes->degreeOfFreedom[1];
459      if( mpi_info->rank!=(mpi_info->size-1) || periodicLocal[1] )          Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, forwardDOF );
460      {          Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, backwardDOF );
461        /* right hand boundary is periodic - add the external reference to the distribution */      }
if( periodicLocal[1] )
{
out->Nodes->Coordinates[INDEX2(0,k,1)]=Length[0];
out->Nodes->Id[k]=k;
out->Nodes->Tag[k]=0;
out->Nodes->degreeOfFreedom[k]=k;
k++;
}
/* right hand boundary with another subdomain, need to add the nodes/DOFs that
defines the element that spans the boundary */
else
{
out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(firstNode+numNodesLocal-periodicLocal[0])/DBLE(N0-1)*Length[0];
out->Nodes->Id[k]=k;
out->Nodes->Tag[k]=0;
out->Nodes->degreeOfFreedom[k]=DOFcount;
k++;
}
}
/* setup boundary DOF data */
if( domInternal )
{
targetDomain = mpi_info->rank-1;
forwardDOF[0] = 0;
backwardDOF[0] = numNodesLocal+1; backwardDOF[1] = numNodesLocal;
Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, forwardDOF );
Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, backwardDOF );

targetDomain = mpi_info->rank+1;
forwardDOF[0] = numNodesLocal-2; forwardDOF[1] = numNodesLocal-1;
backwardDOF[0] = numNodesLocal+2;
Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, forwardDOF );
Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, backwardDOF );
}
else if( mpi_info->size>2 || (mpi_info->size==2&&!periodic[0]) )
if( domLeft )
{
if( periodicLocal[0] )
{
targetDomain = mpi_info->size-1;
forwardDOF[0] = 0;
backwardDOF[0] = numNodesLocal; backwardDOF[1] = numNodesLocal-1;
Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, forwardDOF );
Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, backwardDOF );
}
targetDomain = mpi_info->rank+1;
forwardDOF[0] = numNodesLocal-2-periodicLocal[0];
forwardDOF[1] = numNodesLocal-1-periodicLocal[0];
backwardDOF[0] = numNodesLocal + periodicLocal[0];
Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, forwardDOF );
Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, backwardDOF );
}
else
{
targetDomain = mpi_info->rank-1;
forwardDOF[0] = 0;
backwardDOF[0] = numNodesLocal+1;
backwardDOF[1] = numNodesLocal;
Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, forwardDOF );
Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, backwardDOF );

if( periodicLocal[1] )
{
targetDomain = 0;
forwardDOF[0] = numNodesLocal-1; forwardDOF[1] = numNodesLocal-1;
backwardDOF[0] = numNodesLocal+1;
Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, forwardDOF );
Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, backwardDOF );
}
}
else if( mpi_info->size==2 && periodic[0])
if( domLeft )
{
targetDomain = mpi_info->size-1;
forwardDOF[0] = 0;
forwardDOF[1] = numDOFLocal-2;
forwardDOF[2] = numDOFLocal-1;
backwardDOF[0] = numDOFLocal+1;
backwardDOF[1] = numDOFLocal;
backwardDOF[2] = numDOFLocal+2;
Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 3, forwardDOF );
Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 3, backwardDOF );
}
else
{
targetDomain = 0;
forwardDOF[0] = numDOFLocal-2;
forwardDOF[1] = numDOFLocal-1;
forwardDOF[2] = 0;
backwardDOF[0] = numDOFLocal+2;
backwardDOF[1] = numDOFLocal+1;
backwardDOF[2] = numDOFLocal;
Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 3, forwardDOF );
Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 3, backwardDOF );
}

if (! Finley_MPI_noError(mpi_info)) {
Finley_Mesh_dealloc(out);
return NULL;
}
out->Nodes->degreeOfFreedomDistribution->numBoundary = DOFBoundary[0] + DOFBoundary[1];
out->Nodes->degreeOfFreedomDistribution->numInternal = numDOFLocal - out->Nodes->degreeOfFreedomDistribution->numBoundary;

/*printf( "\n============NODES=============\n" );
for( k=0; k<numNodesLocal; k++ )
printf( "\tI\tId %d\tDOF %d\tcoord [%g]\n", out->Nodes->Id[k], out->Nodes->degreeOfFreedom[k] , out->Nodes->Coordinates[INDEX2(0,k,1)] );
for( k=numNodesLocal; k<numNodesLocal+numNodesExternal; k++ )
printf( "\tE\tId %d\tDOF %d\tcoord [%g]\n", out->Nodes->Id[k], out->Nodes->degreeOfFreedom[k] , out->Nodes->Coordinates[INDEX2(0,k,1)] );

for( k=0; k<out->Nodes->degreeOfFreedomDistribution->numNeighbours; k++ )
{
if( out->Nodes->degreeOfFreedomDistribution->neighbours[k]>=0 )
{
printf( "\t%d boundary DOF { ", out->Nodes->degreeOfFreedomDistribution->edges[k]->numForward );
for( i0=0; i0<out->Nodes->degreeOfFreedomDistribution->edges[k]->numForward; i0++ )
printf( "%d ", out->Nodes->degreeOfFreedomDistribution->edges[k]->indexForward[i0] );
printf("} to %d\n", out->Nodes->degreeOfFreedomDistribution->neighbours[k] );
printf( "\t%d boundary DOF { ", out->Nodes->degreeOfFreedomDistribution->edges[k]->numBackward );
for( i0=0; i0<out->Nodes->degreeOfFreedomDistribution->edges[k]->numBackward; i0++ )
printf( "%d ", out->Nodes->degreeOfFreedomDistribution->edges[k]->indexBackward[i0] );
printf("} from %d\n", out->Nodes->degreeOfFreedomDistribution->neighbours[k] );
}
}*/
}
462
463    /*   set the elements: */    /*   set the elements: */
464    /*   form internal elements */    /*   form internal elements */
465    //#pragma omp parallel for private(i0,k)    #pragma omp parallel for private(i0,k)
466    for (i0=0;i0<numElementsInternal;i0++)    for (i0=0;i0<numElementsLocal;i0++)
467    {    {
468      k=i0;      k=i0;
469      node0=2*i0;          node0 = (periodicLocal[0] && !i0) ? 0 :  2*i0 + periodicLocal[0];
470      out->Elements->Id[k]=k;      out->Elements->Id[k]=k;
471      out->Elements->Tag[k]=0;      out->Elements->Tag[k]=0;
472      out->Elements->Color[k]=0;      out->Elements->Color[k]=COLOR_MOD(i0);
473            out->Elements->Dom[k]=ELEMENT_INTERNAL;
474
475      out->Elements->Nodes[INDEX2(0,k,3)]=node0;      out->Elements->Nodes[INDEX2(0,k,3)]=node0;
476      out->Elements->Nodes[INDEX2(1,k,3)]=node0+2;      out->Elements->Nodes[INDEX2(1,k,3)]=node0+2;
477      out->Elements->Nodes[INDEX2(2,k,3)]=node0+1;      out->Elements->Nodes[INDEX2(2,k,3)]=node0+1;
478    }    }
479
480    /* followed by boundary elements... */      if( boundaryLeft )
481    i0 = numElementsInternal;          out->Elements->Dom[0] = ELEMENT_BOUNDARY;
482    if( mpi_info->size>1 )      if( boundaryRight )
483    {          out->Elements->Dom[numElementsLocal-1] = ELEMENT_BOUNDARY;
/* left hand boundary */
if( mpi_info->rank>0 ) /* left hand boundary is an internal boundary */
{
k=i0;
node0=numNodesLocal;
out->Elements->Id[k]=k;
out->Elements->Tag[k]=0;
out->Elements->Color[k]=0;

out->Elements->Nodes[INDEX2(0,k,3)]=node0+1;
out->Elements->Nodes[INDEX2(1,k,3)]=0;
out->Elements->Nodes[INDEX2(2,k,3)]=node0;
i0++;
}
else if( periodicLocal[0] ) /* left hand boundary is a periodic boundary */
{
k=i0;
node0 = numNodesLocal;
out->Elements->Id[k]=k;
out->Elements->Tag[k]=0;
out->Elements->Color[k]=0;

out->Elements->Nodes[INDEX2(0,k,3)]=node0+1;
out->Elements->Nodes[INDEX2(1,k,3)]=node0-1;
out->Elements->Nodes[INDEX2(2,k,3)]=node0;
i0++;
}
484
/* right hand boundary */
if( mpi_info->rank<mpi_info->size-1 ) /* right hand boundary is an internal boundary */
{
k=i0;
out->Elements->Id[k]=k;
out->Elements->Tag[k]=0;
out->Elements->Color[k]=0;

out->Elements->Nodes[INDEX2(0,k,3)]=numNodesLocal-2-periodicLocal[0];
out->Elements->Nodes[INDEX2(1,k,3)]=nodesExternal[0]+numNodesLocal;
out->Elements->Nodes[INDEX2(2,k,3)]=numNodesLocal-1-periodicLocal[0];
}
else if( periodicLocal[1] ) /* right hand boundary is a periodic boundary */
{
/* no need to do anything */;
k=i0;
out->Elements->Id[k]=k;
out->Elements->Tag[k]=0;
out->Elements->Color[k]=0;

out->Elements->Nodes[INDEX2(0,k,3)]=numNodesLocal-2;
out->Elements->Nodes[INDEX2(1,k,3)]=numNodesLocal+2;
out->Elements->Nodes[INDEX2(2,k,3)]=numNodesLocal-1;
}
}
485    out->Elements->minColor=0;    out->Elements->minColor=0;
486    out->Elements->maxColor=0;    out->Elements->maxColor=COLOR_MOD(0);
487
488    out->Elements->elementDistribution->numLocal    = numElementsLocal;      Finley_ElementFile_setDomainFlags( out->Elements );
out->Elements->elementDistribution->numInternal = numElementsInternal;
out->Elements->elementDistribution->numBoundary = numElementsLocal - numElementsInternal;
489
490    /*   face elements: */    /*   face elements: */
491    if (useElementsOnFace) {    if (useElementsOnFace) {
# Line 628  Finley_Mesh* Finley_RectangularMesh_Line Line 493  Finley_Mesh* Finley_RectangularMesh_Line
493    } else {    } else {
494       NUMNODES=1;       NUMNODES=1;
495    }    }
496        i0 = numElementsLocal;
497    if ( domLeft && !periodicLocal[0] )    if ( domLeft && !periodicLocal[0] )
498    {    {
499      out->FaceElements->Id[0]=i0-1;      out->FaceElements->Id[0]=i0;
500      out->FaceElements->Tag[0]=1;      out->FaceElements->Tag[0]=1;
501      out->FaceElements->Color[0]=0;      out->FaceElements->Color[0]=0;
502            out->FaceElements->Dom[0]=ELEMENT_INTERNAL;
503      if (useElementsOnFace) {      if (useElementsOnFace) {
504         out->FaceElements->Nodes[INDEX2(0,0,NUMNODES)]=0;         out->FaceElements->Nodes[INDEX2(0,0,NUMNODES)]=0;
505         out->FaceElements->Nodes[INDEX2(1,0,NUMNODES)]=2;         out->FaceElements->Nodes[INDEX2(1,0,NUMNODES)]=2;
# Line 647  Finley_Mesh* Finley_RectangularMesh_Line Line 514  Finley_Mesh* Finley_RectangularMesh_Line
514      out->FaceElements->Id[domLeft]=i0;      out->FaceElements->Id[domLeft]=i0;
515      out->FaceElements->Tag[domLeft]=2;      out->FaceElements->Tag[domLeft]=2;
516      out->FaceElements->Color[domLeft]=1;      out->FaceElements->Color[domLeft]=1;
517            out->FaceElements->Dom[domLeft]=ELEMENT_INTERNAL;
518      if (useElementsOnFace) {      if (useElementsOnFace) {
519         out->FaceElements->Nodes[INDEX2(0,domLeft,NUMNODES)]=numNodesLocal-3;         out->FaceElements->Nodes[INDEX2(0,domLeft,NUMNODES)]=N0t-1;
520         out->FaceElements->Nodes[INDEX2(1,domLeft,NUMNODES)]=numNodesLocal-1;         out->FaceElements->Nodes[INDEX2(1,domLeft,NUMNODES)]=N0t-3;
521         out->FaceElements->Nodes[INDEX2(2,domLeft,NUMNODES)]=numNodesLocal-2;         out->FaceElements->Nodes[INDEX2(2,domLeft,NUMNODES)]=N0t-2;
522      } else {      } else {
523         out->FaceElements->Nodes[INDEX2(0,domLeft,NUMNODES)]=numNodesLocal-1;         out->FaceElements->Nodes[INDEX2(0,domLeft,NUMNODES)]=N0t-1;
524      }      }
525    }    }
526        out->FaceElements->numElements=NFaceElements;
527    out->FaceElements->maxColor=0;    out->FaceElements->maxColor=0;
528    out->FaceElements->minColor=0;    out->FaceElements->minColor=1;
529    out->FaceElements->elementDistribution->numBoundary = 0;
530    if( domLeft || domRight && !periodic[0] )      Finley_ElementFile_setDomainFlags( out->FaceElements );
531      out->FaceElements->elementDistribution->numLocal = out->FaceElements->elementDistribution->numInternal = domLeft + domRight;
else
out->FaceElements->elementDistribution->numLocal = out->FaceElements->elementDistribution->numInternal = 0;

532    /* setup distribution info for other elements */    /* setup distribution info for other elements */
533    out->ContactElements->elementDistribution->numLocal = out->ContactElements->elementDistribution->numInternal = out->ContactElements->elementDistribution->numInternal = 0;      Finley_ElementFile_setDomainFlags( out->ContactElements );
534    out->Points->elementDistribution->numLocal = out->Points->elementDistribution->numInternal = out->Points->elementDistribution->numInternal = 0;      Finley_ElementFile_setDomainFlags( out->Points );
535
536    /*printf( "\n============ELEMENTS (%d)=============\n", out->Elements->numElements );      /* reorder the degrees of freedom */
537    for( k=0; k<out->Elements->elementDistribution->numInternal; k++ )      Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
538    {
printf( "I\tId %d : nodes [%d %d %d]->DOF [%d %d %d]\n", out->Elements->Id[k],  out->Elements->Nodes[INDEX2(0,k,3)], out->Elements->Nodes[INDEX2(1,k,3)], out->Elements->Nodes[INDEX2(2,k,3)], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(0,k,3)]], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(1,k,3)]], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(2,k,3)]] );
}
for( k=out->Elements->elementDistribution->numInternal; k<out->Elements->elementDistribution->numLocal; k++ )
{
printf( "B\tId %d : nodes [%d %d %d]->DOF [%d %d %d]\n", out->Elements->Id[k],  out->Elements->Nodes[INDEX2(0,k,3)], out->Elements->Nodes[INDEX2(1,k,3)], out->Elements->Nodes[INDEX2(2,k,3)], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(0,k,3)]], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(1,k,3)]], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(2,k,3)]] );
}
printf( "\n" );
for( k=0; k<out->FaceElements->numElements; k++ )
{
if( NUMNODES==3 )
printf( "F\tId %d : nodes [%d %d %d]->DOF [%d %d %d]\n", out->FaceElements->Id[k],  out->FaceElements->Nodes[INDEX2(0,k,3)], out->FaceElements->Nodes[INDEX2(1,k,3)], out->FaceElements->Nodes[INDEX2(2,k,3)], out->Nodes->degreeOfFreedom[out->FaceElements->Nodes[INDEX2(0,k,3)]], out->Nodes->degreeOfFreedom[out->FaceElements->Nodes[INDEX2(1,k,3)]], out->Nodes->degreeOfFreedom[out->FaceElements->Nodes[INDEX2(2,k,3)]] );
else
printf( "F\tId %d : nodes [%d]->DOF [%d]\n", out->FaceElements->Id[k], out->FaceElements->Nodes[INDEX2(0,k,1)], out->Nodes->degreeOfFreedom[out->FaceElements->Nodes[INDEX2(0,k,1)]] );
}*/
/*  face elements done: */

539    /*   condense the nodes: */    /*   condense the nodes: */

540    Finley_Mesh_resolveNodeIds(out);    Finley_Mesh_resolveNodeIds(out);
541      if( !Finley_MPI_noError(mpi_info) )
542      {
543        Paso_MPIInfo_dealloc( mpi_info );
544        Finley_Mesh_dealloc(out);
545        return NULL;
546      }
547
548    /* setup the CommBuffer */    /* setup the CommBuffer */
549    Finley_NodeDistribution_formCommBuffer( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );    Finley_NodeDistribution_formCommBuffer( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
550    if ( !Finley_MPI_noError( mpi_info )) {    if ( !Finley_MPI_noError( mpi_info )) {
if( Finley_noError() )
Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );
551      Paso_MPIInfo_dealloc( mpi_info );      Paso_MPIInfo_dealloc( mpi_info );
552      Finley_Mesh_dealloc(out);      Finley_Mesh_dealloc(out);
553      return NULL;      return NULL;
554    }    }
555
556    Finley_NodeDistribution_calculateIndexExternal( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );    Finley_NodeDistribution_calculateIndexExternal( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
557      if( !Finley_MPI_noError(mpi_info) )
558      {
559        Paso_MPIInfo_dealloc( mpi_info );
560        Finley_Mesh_dealloc(out);
561        return NULL;
562      }
563
564    /* prepare mesh for further calculatuions:*/    /* prepare mesh for further calculatuions:*/
565    Finley_Mesh_prepare(out);    Finley_Mesh_prepare(out);

566    if( !Finley_MPI_noError(mpi_info) )    if( !Finley_MPI_noError(mpi_info) )
567    {    {
if( Finley_noError() )
Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );
568      Paso_MPIInfo_dealloc( mpi_info );      Paso_MPIInfo_dealloc( mpi_info );
569      Finley_Mesh_dealloc(out);      Finley_Mesh_dealloc(out);
570      return NULL;      return NULL;
571    }    }
572
573    /* free up memory */    /* free up memory */
574    Paso_MPIInfo_dealloc( mpi_info );    Paso_MPIInfo_dealloc( mpi_info );
575
576        //print_mesh_statistics( out, TRUE );
577
578    #ifdef Finley_TRACE    #ifdef Finley_TRACE
579    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
580    #endif    #endif

Legend:
 Removed from v.781 changed lines Added in v.782