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

revision 750 by bcumming, Mon May 15 04:03:49 2006 UTC revision 751 by bcumming, Mon Jun 26 01:46:34 2006 UTC
# Line 26  Line 26
26
27  #include "RectangularMesh.h"  #include "RectangularMesh.h"
28
/**************************************************************/
29
30  Finley_Mesh* Finley_RectangularMesh_Line3(dim_t* numElements,double* Length,bool_t* periodic,index_t order,bool_t useElementsOnFace) {
31    /**************************************************************
32      The code for Mesh_line[2/3].c is messy, it was used as a test
33      bed for getting the MPI distributed mesh figured out before
34      doing the rectangle/brick domains. The code is efficient, it just
35      uses too many if/else if/else statements. Unfortunately, calculating
36      a distributed mesh is an easy thing to visualise, but a devil to implement!
37
38      Sorry if I don't get around
39      to neatening this code up, but rest assured, the more important
40      (and complex) rectangle and brick codes are much better organised.
41    **************************************************************/
42
43    #ifdef PASO_MPI
44    /* get the number of nodes/elements for domain with rank=rank, of size processors
45       where n is the total number of nodes/elements in the global domain */
46    static index_t domain_MODdim( index_t rank, index_t size, index_t n )
47    {
48      rank = size-rank-1;
49
50      if( rank < n%size )
51        return (index_t)floor(n/size)+1;
52      return (index_t)floor(n/size);
53    }
54
55
56    /* Determines the number of nodes/elements etc along an axis which is numElementsGlobal long for domain rank */
57    /* A bit messy, but it only has to be done once... */
58    static void domain_calculateDimension( index_t rank, dim_t size, dim_t numElementsGlobal, bool_t periodic, dim_t *numNodesLocal, dim_t *numDOFLocal, dim_t *numElementsLocal, dim_t *numElementsInternal, dim_t *firstNode, dim_t *nodesExternal, dim_t *DOFExternal, dim_t *numNodesExternal, bool_t *periodicLocal, dim_t *DOFBoundary )
59    {
60      index_t i0;
61      dim_t numNodesGlobal = numElementsGlobal+1;
62
63      (*numNodesLocal) = 2*domain_MODdim( rank, size, numNodesGlobal ) - 1;
64      if( rank<size-1 ) // add on node for right hand boundary
65        (*numNodesLocal) += 1;
66
67      numElementsLocal[0] = domain_MODdim( rank, size, numNodesGlobal )+1;
68      periodicLocal[0] = periodicLocal[1] = FALSE;
69      nodesExternal[0] = nodesExternal[1] = 1;
70      if( periodic )
71      {
72        if( size==1 )
73        {
74          numElementsLocal[0] = numElementsGlobal;
75          nodesExternal[0] = nodesExternal[1] = 0;
76          periodicLocal[0] = periodicLocal[1] = TRUE;
77        }
78        else
79        {
80          if( rank==0 )
81          {
82            periodicLocal[0] = TRUE;
83            numNodesLocal[0]++;
84          }
85          if( rank==(size-1) )
86          {
87            periodicLocal[1] = TRUE;
88            numNodesLocal[0]--;
89            numElementsLocal[0]--;
90          }
91        }
92      }
93      else if( !periodic )
94      {
95        if( rank==0 ){
96          nodesExternal[0]--;
97          numElementsLocal[0]--;
98        }
99        if( rank==(size-1) )
100        {
101          nodesExternal[1]--;
102          numElementsLocal[0]--;
103        }
104      }
105      nodesExternal[0]*=2;
106      numNodesExternal[0] = nodesExternal[0]+nodesExternal[1];
107
108      numElementsInternal[0] = numElementsLocal[0];
109      if( (rank==0) && (rank==size-1) );
110
111      else if( !periodic && ( (rank==0) ^ (rank==size-1) ) )
112          numElementsInternal[0] -= 1;
113      else
114        numElementsInternal[0] -= 2;
115
116      firstNode[0] = 0;
117      for( i0=0; i0<rank; i0++ )
118        firstNode[0] += domain_MODdim( i0, size, numNodesGlobal );
119      firstNode[0] *= 2;
120
121      numDOFLocal[0] = numNodesLocal[0];
122      if( periodicLocal[0] )
123        numDOFLocal[0]--;
124
125      DOFExternal[0] = nodesExternal[0];
126      DOFExternal[1] = nodesExternal[1];
127
128      if( size>1 )
129      {
130        DOFBoundary[0] = periodicLocal[0]*1 + (rank>0)*1;
131        DOFBoundary[1] = periodicLocal[1]*2 + (rank<(size-1))*2;
132      }
133      else
134      {
135        DOFBoundary[0] = DOFBoundary[1] = 0;
136      }
137      /* some debugging printf statements */
138      //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] );
139    }
140    #endif
141
142
143    Finley_Mesh* Finley_RectangularMesh_Line3(dim_t* numElements,double* Length,bool_t* periodic,index_t order,bool_t useElementsOnFace)
144    #ifndef PASO_MPI
145    {
146    dim_t N0,NE0,i0,NDOF0,NFaceElements,NUMNODES;    dim_t N0,NE0,i0,NDOF0,NFaceElements,NUMNODES;
147    index_t k,node0;    index_t k,node0;
148    Finley_Mesh* out;    Finley_Mesh* out;
# Line 48  Finley_Mesh* Finley_RectangularMesh_Line Line 162  Finley_Mesh* Finley_RectangularMesh_Line
162
163    sprintf(name,"Rectangular mesh with %d nodes",N0);    sprintf(name,"Rectangular mesh with %d nodes",N0);
164    /* TEMPFIX */    /* TEMPFIX */
165  #ifndef PASO_MPI
166    out=Finley_Mesh_alloc(name,1,order);    out=Finley_Mesh_alloc(name,1,order);
167    if (! Finley_noError()) return NULL;    if (! Finley_noError()) return NULL;
168
# Line 61  Finley_Mesh* Finley_RectangularMesh_Line Line 175  Finley_Mesh* Finley_RectangularMesh_Line
175      out->ContactElements=Finley_ElementFile_alloc(Point1_Contact,out->order);      out->ContactElements=Finley_ElementFile_alloc(Point1_Contact,out->order);
176    }    }
177    out->Points=Finley_ElementFile_alloc(Point1,out->order);    out->Points=Finley_ElementFile_alloc(Point1,out->order);
#else
/* TODO */
PASO_MPI_TODO;
out = NULL;
#endif
178    if (! Finley_noError()) {    if (! Finley_noError()) {
179         Finley_Mesh_dealloc(out);         Finley_Mesh_dealloc(out);
180         return NULL;         return NULL;
181    }    }
182
183    /*  allocate tables: */    /*  allocate tables: */
184   #ifndef PASO_MPI
185    Finley_NodeFile_allocTable(out->Nodes,N0);    Finley_NodeFile_allocTable(out->Nodes,N0);
186  #else
/* TODO */
#endif
187    Finley_ElementFile_allocTable(out->Elements,NE0);    Finley_ElementFile_allocTable(out->Elements,NE0);
188    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
189    if (! Finley_noError()) {    if (! Finley_noError()) {
# Line 170  Finley_Mesh* Finley_RectangularMesh_Line Line 277  Finley_Mesh* Finley_RectangularMesh_Line
277    }    }
278    return out;    return out;
279  }  }
280    #else
281    /* MPI version */
282    {
283      dim_t N0, NE0, NE0_local, i0, NDOF0, NFaceElements, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal, DOFBoundary[2];
284      index_t *numForward=NULL, *numBackward=NULL;
285      index_t NUMNODES, node0,k, i,firstNode=0, DOFcount=0, forwardDOF[4], backwardDOF[4];
286      index_t targetDomain=-1;
287      bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE;
288      Finley_Mesh* out=NULL;
289      char name[50];
290      double time0=Finley_timer();
291      Paso_MPIInfo *mpi_info = NULL;
292
293      /* dimensions the global domain */
294      NE0=MAX(1,numElements[0]);
295      N0=2*NE0+1;
296      NDOF0=N0;
297      if( periodic[0] )
298        NDOF0--;
299
300      /* get MPI information */
301      mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
302      if (! Finley_noError()) {
303            Finley_Mesh_dealloc(out);
304            return NULL;
305      }
306
307      if( mpi_info->rank==0 )
308        domLeft = TRUE;
309      if( mpi_info->rank==mpi_info->size-1 )
310        domRight = TRUE;
311      if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
312        domInternal = TRUE;
313
314      /* dimensions of the local subdomain */
315      domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal, DOFBoundary );
316
317      /* form face element if the local domain contains a periodic boundary */
318      NFaceElements = 0;
319      if( !periodic[0] )
320      {
321        if(domLeft)
322          NFaceElements++;
323        if(domRight)
324          NFaceElements++;
325      }
326
327      /*  allocate mesh: */
328      sprintf(name,"Rectangular mesh with %d nodes",N0);
329      out=Finley_Mesh_alloc(name,1,order,mpi_info);
330      if (! Finley_noError()) return NULL;
331
332      out->Elements=Finley_ElementFile_alloc(Line3,out->order,mpi_info);
333      if (useElementsOnFace) {
334         out->FaceElements=Finley_ElementFile_alloc(Line3Face,out->order,mpi_info);
335         out->ContactElements=Finley_ElementFile_alloc(Line3Face_Contact,out->order,mpi_info);
336      } else {
337         out->FaceElements=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
338         out->ContactElements=Finley_ElementFile_alloc(Point1_Contact,out->order,mpi_info);
339      }
340      out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
341      if (! Finley_noError()) {
342            Finley_Mesh_dealloc(out);
343            return NULL;
344      }
345
346      /*  allocate tables: */
347      Finley_NodeFile_allocTable(out->Nodes,numNodesLocal+numNodesExternal);
348      Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal, DOFExternal[0]+DOFExternal[1], 0 );
349      Finley_ElementFile_allocTable(out->Elements,numElementsLocal);
350      if( NFaceElements )
351        Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
352      if (! Finley_noError()) {
353          Finley_Mesh_dealloc(out);
354          return NULL;
355      }
356
357      /*  set nodes: */
358      #pragma omp parallel for private(i0,k)
359      /* local nodes */
360      for (i0=0;i0<numNodesLocal;i0++) {
361         k=i0;
362         out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(i0+firstNode)/DBLE(N0-1)*Length[0];
363         out->Nodes->Id[k]=k;
364         out->Nodes->Tag[k]=0;
365         out->Nodes->degreeOfFreedom[k]=i0%(numDOFLocal);
366      }
367
368      /* external nodes - do left then right hand side */
369      /* the following only applies if more than one domain */
370
371      /* this is messy, could be done cleaner */
372      if( mpi_info->size>1 )
373      {
374        DOFcount = numNodesLocal;
375        k=numNodesLocal;
376        if( mpi_info->rank!=0 || periodicLocal[0] )
377        {
378          /* left hand boundary is periodic -
379             add the nodes/DOF that define the element on the right hand boundary */
380          if( periodicLocal[0] )
381          {
382            k--;
383            out->Nodes->Coordinates[INDEX2(0,k,1)]=Length[0];
384            out->Nodes->Id[k]=k;
385            out->Nodes->Tag[k]=0;
386            out->Nodes->degreeOfFreedom[k]=0;
387            DOFcount--;
388            k++;
389            out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(N0-2)/DBLE(N0-1)*Length[0];
390            out->Nodes->Id[k]=k;
391            out->Nodes->Tag[k]=0;
392            out->Nodes->degreeOfFreedom[k]=DOFcount++;
393            k++;
394            out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(N0-3)/DBLE(N0-1)*Length[0];
395            out->Nodes->Id[k]=k;
396            out->Nodes->Tag[k]=0;
397            out->Nodes->degreeOfFreedom[k]=DOFcount++;
398            k++;
399          }
400          /* left hand boundary with another subdomain, need to add the nodes/DOFs that
401             defines the element that spans the boundary */
402          else
403          {
404            out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(firstNode-1)/DBLE(N0-1)*Length[0];
405            out->Nodes->Id[k]=k;
406            out->Nodes->Tag[k]=0;
407            out->Nodes->degreeOfFreedom[k]=DOFcount++;
408            k++;
409            out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(firstNode-2)/DBLE(N0-1)*Length[0];
410            out->Nodes->Id[k]=k;
411            out->Nodes->Tag[k]=0;
412            out->Nodes->degreeOfFreedom[k]=DOFcount++;
413            k++;
414          }
415        }
416        if( mpi_info->rank!=(mpi_info->size-1) || periodicLocal[1] )
417        {
418          /* right hand boundary is periodic - add the external reference to the distribution */
419          if( periodicLocal[1] )
420          {
421            out->Nodes->Coordinates[INDEX2(0,k,1)]=Length[0];
422            out->Nodes->Id[k]=k;
423            out->Nodes->Tag[k]=0;
424            out->Nodes->degreeOfFreedom[k]=k;
425            k++;
426          }
427          /* right hand boundary with another subdomain, need to add the nodes/DOFs that
428             defines the element that spans the boundary */
429          else
430          {
431            out->Nodes->Coordinates[INDEX2(0,k,1)]=DBLE(firstNode+numNodesLocal-periodicLocal[0])/DBLE(N0-1)*Length[0];
432            out->Nodes->Id[k]=k;
433            out->Nodes->Tag[k]=0;
434            out->Nodes->degreeOfFreedom[k]=DOFcount;
435            k++;
436          }
437        }
438        /* setup boundary DOF data */
439        if( domInternal )
440        {
441          targetDomain = mpi_info->rank-1;
442          forwardDOF[0] = 0;
443          backwardDOF[0] = numNodesLocal+1; backwardDOF[1] = numNodesLocal;
444          Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, forwardDOF );
445          Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, backwardDOF );
446
447          targetDomain = mpi_info->rank+1;
448          forwardDOF[0] = numNodesLocal-2; forwardDOF[1] = numNodesLocal-1;
449          backwardDOF[0] = numNodesLocal+2;
450          Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, forwardDOF );
451          Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, backwardDOF );
452        }
453        else if( mpi_info->size>2 || (mpi_info->size==2&&!periodic[0]) )
454          if( domLeft )
455          {
456            if( periodicLocal[0] )
457            {
458              targetDomain = mpi_info->size-1;
459              forwardDOF[0] = 0;
460              backwardDOF[0] = numNodesLocal; backwardDOF[1] = numNodesLocal-1;
461              Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, forwardDOF );
462              Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, backwardDOF );
463            }
464            targetDomain = mpi_info->rank+1;
465            forwardDOF[0] = numNodesLocal-2-periodicLocal[0];
466            forwardDOF[1] = numNodesLocal-1-periodicLocal[0];
467            backwardDOF[0] = numNodesLocal + periodicLocal[0];
468            Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, forwardDOF );
469            Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, backwardDOF );
470          }
471          else
472          {
473            targetDomain = mpi_info->rank-1;
474            forwardDOF[0] = 0;
475            backwardDOF[0] = numNodesLocal+1;
476            backwardDOF[1] = numNodesLocal;
477            Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, forwardDOF );
478            Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, backwardDOF );
479
480            if( periodicLocal[1] )
481            {
482              targetDomain = 0;
483              forwardDOF[0] = numNodesLocal-1; forwardDOF[1] = numNodesLocal-1;
484              backwardDOF[0] = numNodesLocal+1;
485              Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 2, forwardDOF );
486              Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 1, backwardDOF );
487            }
488          }
489        else if( mpi_info->size==2 && periodic[0])
490          if( domLeft )
491          {
492              targetDomain = mpi_info->size-1;
493              forwardDOF[0] = 0;
494              forwardDOF[1] = numDOFLocal-2;
495              forwardDOF[2] = numDOFLocal-1;
496              backwardDOF[0] = numDOFLocal+1;
497              backwardDOF[1] = numDOFLocal;
498              backwardDOF[2] = numDOFLocal+2;
499              Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 3, forwardDOF );
500              Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 3, backwardDOF );
501          }
502          else
503          {
504              targetDomain = 0;
505              forwardDOF[0] = numDOFLocal-2;
506              forwardDOF[1] = numDOFLocal-1;
507              forwardDOF[2] = 0;
508              backwardDOF[0] = numDOFLocal+2;
509              backwardDOF[1] = numDOFLocal+1;
510              backwardDOF[2] = numDOFLocal;
511              Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 3, forwardDOF );
512              Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, 3, backwardDOF );
513          }
514
515
516        if (! Finley_MPI_noError(mpi_info)) {
517          Finley_Mesh_dealloc(out);
518          return NULL;
519        }
520        out->Nodes->degreeOfFreedomDistribution->numBoundary = DOFBoundary[0] + DOFBoundary[1];
521        out->Nodes->degreeOfFreedomDistribution->numInternal = numDOFLocal - out->Nodes->degreeOfFreedomDistribution->numBoundary;
522
523        /*printf( "\n============NODES=============\n" );
524        for( k=0; k<numNodesLocal; k++ )
525          printf( "\tI\tId %d\tDOF %d\tcoord [%g]\n", out->Nodes->Id[k], out->Nodes->degreeOfFreedom[k] , out->Nodes->Coordinates[INDEX2(0,k,1)] );
526        for( k=numNodesLocal; k<numNodesLocal+numNodesExternal; k++ )
527          printf( "\tE\tId %d\tDOF %d\tcoord [%g]\n", out->Nodes->Id[k], out->Nodes->degreeOfFreedom[k] , out->Nodes->Coordinates[INDEX2(0,k,1)] );
528
529        for( k=0; k<out->Nodes->degreeOfFreedomDistribution->numNeighbours; k++ )
530        {
531          if( out->Nodes->degreeOfFreedomDistribution->neighbours[k]>=0 )
532          {
533            printf( "\t%d boundary DOF { ", out->Nodes->degreeOfFreedomDistribution->edges[k]->numForward );
534            for( i0=0; i0<out->Nodes->degreeOfFreedomDistribution->edges[k]->numForward; i0++ )
535              printf( "%d ", out->Nodes->degreeOfFreedomDistribution->edges[k]->indexForward[i0] );
536            printf("} to %d\n", out->Nodes->degreeOfFreedomDistribution->neighbours[k] );
537            printf( "\t%d boundary DOF { ", out->Nodes->degreeOfFreedomDistribution->edges[k]->numBackward );
538            for( i0=0; i0<out->Nodes->degreeOfFreedomDistribution->edges[k]->numBackward; i0++ )
539              printf( "%d ", out->Nodes->degreeOfFreedomDistribution->edges[k]->indexBackward[i0] );
540            printf("} from %d\n", out->Nodes->degreeOfFreedomDistribution->neighbours[k] );
541          }
542        }*/
543      }
544
545      /*   set the elements: */
546      /*   form internal elements */
547      //#pragma omp parallel for private(i0,k)
548      for (i0=0;i0<numElementsInternal;i0++)
549      {
550        k=i0;
551        node0=2*i0;
552        out->Elements->Id[k]=k;
553        out->Elements->Tag[k]=0;
554        out->Elements->Color[k]=0;
555
556        out->Elements->Nodes[INDEX2(0,k,3)]=node0;
557        out->Elements->Nodes[INDEX2(1,k,3)]=node0+2;
558        out->Elements->Nodes[INDEX2(2,k,3)]=node0+1;
559      }
560
561      /* followed by boundary elements... */
562      i0 = numElementsInternal;
563      if( mpi_info->size>1 )
564      {
565        /* left hand boundary */
566        if( mpi_info->rank>0 ) /* left hand boundary is an internal boundary */
567        {
568          k=i0;
569          node0=numNodesLocal;
570          out->Elements->Id[k]=k;
571          out->Elements->Tag[k]=0;
572          out->Elements->Color[k]=0;
573
574          out->Elements->Nodes[INDEX2(0,k,3)]=node0+1;
575          out->Elements->Nodes[INDEX2(1,k,3)]=0;
576          out->Elements->Nodes[INDEX2(2,k,3)]=node0;
577          i0++;
578        }
579        else if( periodicLocal[0] ) /* left hand boundary is a periodic boundary */
580        {
581          k=i0;
582          node0 = numNodesLocal;
583          out->Elements->Id[k]=k;
584          out->Elements->Tag[k]=0;
585          out->Elements->Color[k]=0;
586
587          out->Elements->Nodes[INDEX2(0,k,3)]=node0+1;
588          out->Elements->Nodes[INDEX2(1,k,3)]=node0-1;
589          out->Elements->Nodes[INDEX2(2,k,3)]=node0;
590          i0++;
591        }
592
593        /* right hand boundary */
594        if( mpi_info->rank<mpi_info->size-1 ) /* right hand boundary is an internal boundary */
595        {
596          k=i0;
597          out->Elements->Id[k]=k;
598          out->Elements->Tag[k]=0;
599          out->Elements->Color[k]=0;
600
601          out->Elements->Nodes[INDEX2(0,k,3)]=numNodesLocal-2-periodicLocal[0];
602          out->Elements->Nodes[INDEX2(1,k,3)]=nodesExternal[0]+numNodesLocal;
603          out->Elements->Nodes[INDEX2(2,k,3)]=numNodesLocal-1-periodicLocal[0];
604        }
605        else if( periodicLocal[1] ) /* right hand boundary is a periodic boundary */
606        {
607          /* no need to do anything */;
608          k=i0;
609          out->Elements->Id[k]=k;
610          out->Elements->Tag[k]=0;
611          out->Elements->Color[k]=0;
612
613          out->Elements->Nodes[INDEX2(0,k,3)]=numNodesLocal-2;
614          out->Elements->Nodes[INDEX2(1,k,3)]=numNodesLocal+2;
615          out->Elements->Nodes[INDEX2(2,k,3)]=numNodesLocal-1;
616        }
617      }
618      out->Elements->minColor=0;
619      out->Elements->maxColor=0;
620
621      out->Elements->elementDistribution->numLocal    = numElementsLocal;
622      out->Elements->elementDistribution->numInternal = numElementsInternal;
623      out->Elements->elementDistribution->numBoundary = numElementsLocal - numElementsInternal;
624
625      /*   face elements: */
626      if (useElementsOnFace) {
627         NUMNODES=3;
628      } else {
629         NUMNODES=1;
630      }
631      if ( domLeft && !periodicLocal[0] )
632      {
633        out->FaceElements->Id[0]=i0-1;
634        out->FaceElements->Tag[0]=1;
635        out->FaceElements->Color[0]=0;
636        if (useElementsOnFace) {
637           out->FaceElements->Nodes[INDEX2(0,0,NUMNODES)]=0;
638           out->FaceElements->Nodes[INDEX2(1,0,NUMNODES)]=2;
639           out->FaceElements->Nodes[INDEX2(2,0,NUMNODES)]=1;
640        } else {
641           out->FaceElements->Nodes[INDEX2(0,0,NUMNODES)]=0;
642        }
643        i0++;
644      }
645      if( domRight && !periodicLocal[1])
646      {
647        out->FaceElements->Id[domLeft]=i0;
648        out->FaceElements->Tag[domLeft]=2;
649        out->FaceElements->Color[domLeft]=1;
650        if (useElementsOnFace) {
651           out->FaceElements->Nodes[INDEX2(0,domLeft,NUMNODES)]=numNodesLocal-3;
652           out->FaceElements->Nodes[INDEX2(1,domLeft,NUMNODES)]=numNodesLocal-1;
653           out->FaceElements->Nodes[INDEX2(2,domLeft,NUMNODES)]=numNodesLocal-2;
654        } else {
655           out->FaceElements->Nodes[INDEX2(0,domLeft,NUMNODES)]=numNodesLocal-1;
656        }
657      }
658      out->FaceElements->maxColor=0;
659      out->FaceElements->minColor=0;
660      out->FaceElements->elementDistribution->numBoundary = 0;
661      if( domLeft || domRight && !periodic[0] )
662        out->FaceElements->elementDistribution->numLocal = out->FaceElements->elementDistribution->numInternal = domLeft + domRight;
663      else
664        out->FaceElements->elementDistribution->numLocal = out->FaceElements->elementDistribution->numInternal = 0;
665
666      /* setup distribution info for other elements */
667      out->ContactElements->elementDistribution->numLocal = out->ContactElements->elementDistribution->numInternal = out->ContactElements->elementDistribution->numInternal = 0;
668      out->Points->elementDistribution->numLocal = out->Points->elementDistribution->numInternal = out->Points->elementDistribution->numInternal = 0;
669
670      /*printf( "\n============ELEMENTS (%d)=============\n", out->Elements->numElements );
671      for( k=0; k<out->Elements->elementDistribution->numInternal; k++ )
672      {
673        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)]] );
674      }
675      for( k=out->Elements->elementDistribution->numInternal; k<out->Elements->elementDistribution->numLocal; k++ )
676      {
677        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)]] );
678      }
679      printf( "\n" );
680      for( k=0; k<out->FaceElements->numElements; k++ )
681      {
682        if( NUMNODES==3 )
683                    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)]] );
684        else
685          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)]] );
686      }*/
687      /*  face elements done: */
688
689      /*   condense the nodes: */
690
691      Finley_Mesh_resolveNodeIds(out);
692
693      /* setup the CommBuffer */
694      Finley_NodeDistribution_formCommBuffer( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
695      if ( !Finley_MPI_noError( mpi_info )) {
696        if( Finley_noError() )
697          Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );
698        Paso_MPIInfo_dealloc( mpi_info );
699        Finley_Mesh_dealloc(out);
700        return NULL;
701      }
702
703      Finley_NodeDistribution_calculateIndexExternal( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
704
705      /* prepare mesh for further calculatuions:*/
706      Finley_Mesh_prepare(out);
707
708      if( !Finley_MPI_noError(mpi_info) )
709      {
710        if( Finley_noError() )
711          Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );
712        Paso_MPIInfo_dealloc( mpi_info );
713        Finley_Mesh_dealloc(out);
714        return NULL;
715      }
716
717      /* free up memory */
718      Paso_MPIInfo_dealloc( mpi_info );
719
720      #ifdef Finley_TRACE
721      printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
722      #endif
723
724      return out;
725    }
726    #endif
727
728  /*  /*
729  * Revision 1.3  2005/09/01 03:31:36  jgs  * Revision 1.3  2005/09/01 03:31:36  jgs

Legend:
 Removed from v.750 changed lines Added in v.751