cancer_exp_3d_nd.h

// #define __UREA__
// #define __SIMPATURS__
//   
// #include <Urea_3d4d.h>
// #include <MacroCell.h>

void ExampleFile()
{
  
  OutPut("Example: cancer_exp_3d.h" << endl); //OutPut("GRID_TYPE set to " << TDatabase::ParamDB->GRID_TYPE << endl);
}

// ========================================================================
// definitions for the temperature
// ========================================================================

void Exact(double x, double y, double z, double *values)
{

  values[0] = 0;
}

void V_Exact(double x, double y, double z, double *values)
{
 values[0] =  0;
}

void W_Exact(double x, double y, double z, double *values)
{
  values[0] =  0;
}

// =================initial conditon=============================================================================
void InitialCondition(double x, double y, double z, double *values)
{ 
//===============cube=============================================================================================
//        double temp; 
//          temp = ((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)+(z-0.5)*(z-0.5));
//   
//    if (temp >0.1)
//     { values[0] = 0.; }
//    else
//     { values[0] = exp(-temp/0.01);}
    
//===============cube-1D-Test Case==================================================================================
//          double temp; 
//          temp = ((x)*(x)+(y)*(y)+(z)*(z));
//   
//         if (temp >0.25)
//          { values[0] = 0.; }
//         else
//          { values[0] = exp(-temp/0.01);}
         
   
//==============Breast Geometry=====================================================================================
  //values[0] = exp(-((x-286.530470)*(x-286.530470)+(y-34.297560)*(y-34.297560)+(z-162.509553)*(z-162.509553))/20); 
    /* values[0] = exp(-((x-239)*(x-239)+(y-24)*(y-24)+(z-167)*(z-167))/20); */
    
     values[0] = exp(-((x-240)*(x-240)+(y-50)*(y-50)+(z-147)*(z-147))/20); 
    
//    values[0] = exp(-((x-225)*(x-225)+(y-5)*(y-5)+(z-150)*(z-150))/100);    //near boundary
// values[0] = exp(-((x-230)*(x-230)+(y+30)*(y+30)+(z-177)*(z-177))/100);    //on boundary
     
     
//==============Breast Geometry - breast2.mesh========================================================================
//   values[0] =  exp(-((x-264.792889)*(x-264.792889)+(y-144.038695)*(y-144.038695)+(z+10.472588)*(z+10.472588))/100); 
  
}

void V_InitialCondition(double x, double y, double z, double *values)
{ 
//===============cube==================================================================================================== 
//       double temp; 
//        temp = ((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)+(z-0.5)*(z-0.5));
//   
//    if (temp >0.1)
//     { values[0] = 1.; }
//    else
//     { values[0] = 1- 0.5 * exp(-temp/0.01);}
  
//===============cube-1D-Test Case=========================================================================================
//     double temp; 
//     temp = ((x)*(x)+(y)*(y)+(z)*(z));
//   
//     if (temp >0.25)
//      { values[0] = 1.; }
//     else
//      { values[0] = 1- 0.5 * exp(-temp/0.01);}
  
//==============Breast Geometry=============================================================================================
  // values[0] = 1 - exp(-((x-286.530470)*(x-286.530470)+(y-34.297560)*(y-34.297560)+(z-162.509553)*(z-162.509553))/20); 
  //    values[0] = 1- exp(-((x-239)*(x-239)+(y-24)*(y-24)+(z-167)*(z-167))/20); 
  
  
   values[0] = 1 - 0.5 * exp(-((x-240)*(x-240)+(y-50)*(y-50)+(z-147)*(z-147))/20); 
   
//==============Breast Geometry - breast2.mesh===============================================================================
//     values[0] = 1- 0.5* exp(-((x-264.792889)*(x-264.792889)+(y-144.038695)*(y-144.038695)+(z+10.472588)*(z+10.472588))/100);  
//          values[0] = 1- 0.5 * exp(-((x-268)*(x-268)+(y-165)*(y-165)+(z+10.)*(z+10.))/20); 
      
}

void W_InitialCondition(double x, double y, double z, double *values)
{
//===============cube========================================================================================================
//      double temp; 
//          temp = ((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)+(z-0.5)*(z-0.5));
//   
//   if (temp >0.1)
//     { values[0] = 0.; }
//   else
//     { values[0] = 0.5 * exp(-temp/0.01);}
  
//===============cube-1D-Test Case=============================================================================================
//     double temp; 
//     temp = ((x)*(x)+(y)*(y)+(z)*(z));
//  
//     if (temp >0.25)
//     { values[0] = 0.; }
//     else
//     { values[0] = 0.5 * exp(-temp/0.01);}
    
   
//==============Breast Geometry=================================================================================================
     //  values[0] = 5 * exp(-((x-286.530470)*(x-286.530470)+(y-34.297560)*(y-34.297560)+(z-162.509553)*(z-162.509553))/20); 
    //   values[0] = 0.5 * exp(-((x-239)*(x-239)+(y-24)*(y-24)+(z-167)*(z-167))/20);
  
  
   values[0] =   0.5 * exp(-((x-240)*(x-240)+(y-50)*(y-50)+(z-147)*(z-147))/20);   
   
//==============Breast Geometry - breast2.mesh==================================================================================
//      values[0] = 0.5 * exp(-((x-264.792889)*(x-264.792889)+(y-144.038695)*(y-144.038695)+(z+10.472588)*(z+10.472588))/100); 
     
       
}


void BoundCondition(double x, double y, double z, BoundCond &cond)
{
  cond = NEUMANN;
}

void V_BoundCondition(double x, double y, double z, BoundCond &cond)
{
  cond = NEUMANN;
}


void W_BoundCondition(double x, double y, double z, BoundCond &cond)
{
  cond = NEUMANN;
}

// value of boundary condition
void BoundValue(double x, double y, double z, double &value)
{
 value = 0; 
}

void V_BoundValue(double x, double y, double z, double &value)
{
 value = 0; 
}

void W_BoundValue(double x, double y, double z, double &value)
{
 value = 0; 
}
// ========================================================================
// BilinearCoeffs for Heat 
// ========================================================================
void BilinearCoeffs(int n_points, double *X, double *Y, double *Z,
               double **parameters, double **coeffs)
{
  double eps= 0.0/TDatabase::ParamDB->RE_NR;
  int i;
  double *coeff;                                  // *param;
  double x, y, z, c, a[3], b[3], s[3], h;
  double t = TDatabase::TimeDB->CURRENTTIME;
  
  b[0] = 0;
  b[1] = 0;
  b[2] = 0;  
  c = 0;  
  for(i=0;i<n_points;i++)
  {
    coeff = coeffs[i];
    // param = parameters[i];

//     x = X[i];
//     y = Y[i];
//     z = Z[i];

    // diffusion
    coeff[0] = eps;
    // convection in x direction
    coeff[1] = b[0];
    // convection in y direction
    coeff[2] = b[1];
    // convection in z direction
    coeff[3] = b[2];
    // reaction
    coeff[4] = c;
     // rhs
    coeff[5] = 0;  
    
   
    coeff[6] = 0;

  }
}

void V_BilinearCoeffs(int n_points, double *X, double *Y, double *Z,
               double **parameters, double **coeffs)
{
  double eps= 0.0/TDatabase::ParamDB->RE_NR;
  int i;
  double *coeff;                                  // *param;
  double x, y, z, c, a[3], b[3], s[3], h;
  double t = TDatabase::TimeDB->CURRENTTIME;
  b[0] = 0;
  b[1] = 0;
  b[2] = 0;  
  c = 0;
  for(i=0;i<n_points;i++)
  {
    coeff = coeffs[i];
    // param = parameters[i];

//     x = X[i];
//     y = Y[i];
//     z = Z[i];

    // diffusion
    coeff[0] = eps;
    // convection in x direction
    coeff[1] = b[0];
    // convection in y direction
    coeff[2] = b[1];
    // convection in z direction
    coeff[3] = b[2];
    // reaction
    coeff[4] = c;
     // rhs
    coeff[5] = 0;
    coeff[6] = 0;

  }
}

void W_BilinearCoeffs(int n_points, double *X, double *Y, double *Z,
               double **parameters, double **coeffs)
{
  double eps= 5./TDatabase::ParamDB->RE_NR;
  int i;
  double *coeff;                                  // *param;
  double x, y, z, c, a[3], b[3], s[3], h;
  double t = TDatabase::TimeDB->CURRENTTIME;
  
  b[0] = 0;
  b[1] = 0;
  b[2] = 0;  
  c = 0;

  for(i=0;i<n_points;i++)
  {
    coeff = coeffs[i];
    // param = parameters[i];

//     x = X[i];
//     y = Y[i];
//     z = Z[i];

    // diffusion
    coeff[0] = eps;
    // convection in x direction
    coeff[1] = b[0];
    // convection in y direction
    coeff[2] = b[1];
    // convection in z direction
    coeff[3] = b[2];
    // reaction
    coeff[4] = c;
     // rhs
    coeff[5] = 0 ;
    coeff[6] = 0;

  }
}
void ReadAdditionalNModes(char *NODE, tetgenio &InAddPts)
{
  int i, N_Vertices, index, N_Atribude, BDMarker;
  char line[100];  
 std::ifstream dat(NODE);

  if (!dat)
  {
    cerr << "cannot open '" << NODE << "' for input" << endl;
    exit(-1);
  }      
  
  dat.getline (line, 99); // first line is command line    
  dat >> N_Vertices >> N_Atribude >> BDMarker;  
  dat.getline (line, 99);  
 
  InAddPts.numberofpoints = N_Vertices;
  InAddPts.pointlist = new double[3*N_Vertices];  

  //read from file
   for(i=0;i<N_Vertices; i++)
    {
     dat.getline (line, 99);
     dat >> index >> InAddPts.pointlist[3*i] >> InAddPts.pointlist[3*i+1] >> InAddPts.pointlist[3*i+2]; 
     cout<< i << " vert X: " <<InAddPts.pointlist[3*i] << " vert Y: " <<InAddPts.pointlist[3*i+1] <<" vert Z: " <<InAddPts.pointlist[3*i+2] <<endl;
 
    }
 
   dat.close();

// cout << "N_Vertices " <<N_Vertices <<endl;
// 
// exit(0);

} // ReadAdditionalNModes



void ReadMeditMeshWithCells(char *SMESH, tetgenio &In)
{  
  int i, j, k, dimension, N_FVert, N_Faces, N_Vertices, N_CVert, N_Cells;
  int BDComp_Min=10000000;
  int *plist;
  
  double x;
  
  char line[100];
  tetgenio::facet *F;
  tetgenio::polygon *P;    

 std::ifstream dat(SMESH);

  if (!dat)
  {
    cerr << "cannot open '" << SMESH << "' for input" << endl;
    exit(-1);
  }     
   
  // check the dimension
  while (!dat.eof())
  {
    dat >> line;
    
    if ( (!strcmp(line, "Dimension"))  ||  (!strcmp(line, "dimension")) ||  (!strcmp(line, "DIMENSION")))
    {
     dat.getline (line, 99);
     dat >> dimension;
     break;
    }    
    // read until end of line
    dat.getline (line, 99);   
  }
  
  if(dimension!=3)
   {
    cerr << "dimension: " << dimension << endl;
    exit(-1);
   }
  
  // find the N_Vertices
  while (!dat.eof())
  {
    dat >> line;
    
    if ( (!strcmp(line, "Vertices")) ||  (!strcmp(line, "vertices"))   ||  (!strcmp(line, "VERTICES"))   ) 
    {
     dat.getline (line, 99);
     dat >> N_Vertices;
     break;
    }    
    // read until end of line
    dat.getline (line, 99);   
  } 
  
  In.numberofpoints = N_Vertices;
  In.pointlist = new double[3*N_Vertices];  

  //read from file
   for(i=0;i<N_Vertices; i++)
    {
     dat.getline (line, 99);
     dat >> x >> In.pointlist[3*i+1] >> In.pointlist[3*i+2];
     
     In.pointlist[3*i] = x - 3.606596999999999830777142; //center point
     
//      if(i==1063 || i==1061 || i==1062 || i==162 || i==1175 || i==517 || i==1064 || i==518 ) // i==1063  as center
     if(i==1062 || i==916 || i==341 || i==914 || i==162 || i==1063 || i==1061) // i==1062  as center
      cout<< i << " vert X: " <<In.pointlist[3*i] << " vert Y: " <<In.pointlist[3*i+1] <<" vert Z: " <<In.pointlist[3*i+2] <<endl;   
    }
    
  // find the N_Triangles
  // find the N_Vertices
  while (!dat.eof())
  {
    dat >> line;
    
    if ( (!strcmp(line, "Triangles")) ||  (!strcmp(line, "triangles"))   ||  (!strcmp(line, "TRIANGLES"))   ) 
    {
     N_FVert = 3;
     dat.getline (line, 99);
     dat >> N_Faces;
     break;
    }    
    else if ( (!strcmp(line, "Quadrilaterals")) ||  (!strcmp(line, "quadrilaterals"))   ||  (!strcmp(line, "QUADRILATERALS"))   ) 
    {
     N_FVert = 4;
     dat.getline (line, 99);
     dat >> N_Faces;
     break;
    }    
   
    // read until end of line
    dat.getline (line, 99);   
  } 
 
    In.numberoffacets = N_Faces;
    In.facetlist = new tetgenio::facet[In.numberoffacets];
    In.facetmarkerlist = new int[In.numberoffacets];
    
    for(i=0;i<N_Faces; i++)
     {
      dat.getline (line, 99);       
      
      F = &In.facetlist[i];
      tetgenio::init(F);      
      F->numberofpolygons = 1;
      F->polygonlist = new tetgenio::polygon[F->numberofpolygons];
      F->numberofholes = 0;
      F->holelist = NULL;
      P = &F->polygonlist[0];
      tetgenio::init(P);
      P->numberofvertices = N_FVert;
      P->vertexlist = new int[P->numberofvertices];
      
      for(j=0;j<N_FVert;j++)
      {
        dat >> k;
        P->vertexlist[j] = k-1;  // c numbering 
      }
      
      dat >>  In.facetmarkerlist[i];        
      
      if(BDComp_Min > In.facetmarkerlist[i])
        BDComp_Min=In.facetmarkerlist[i];
      
      // cout << i<<  " P->vertexlist[j]:  " <<  In.facetmarkerlist[i]  << endl;  
     } //   for(i=0;i<N_Faces; i++)
   
    BDComp_Min--;
   
    for(i=0;i<N_Faces; i++)
      In.facetmarkerlist[i] -= BDComp_Min;
         
    
  while (!dat.eof())
  {
    dat >> line;
    
    if ( (!strcmp(line, "Tetrahedron")) ||  (!strcmp(line, "tetrahedron"))   ||  (!strcmp(line, "TETRAHEDRON"))   ) 
    {
     N_CVert = 4;
     dat.getline (line, 99);
     dat >> N_Cells;
     break;
    }    
    // read until end of line
    dat.getline (line, 99);   
  }     
 
     In.numberoftetrahedra = N_Cells;
     In.numberofcorners = N_CVert;
     In.numberoftetrahedronattributes = 1;
     In.tetrahedronlist = new int[N_Cells * 4];
     In.tetrahedronattributelist = new REAL[N_Cells];   
     In.tetrahedronvolumelist = new REAL[N_Cells];  
     
     for(i=0; i<N_Cells; i++) 
      {  
       dat.getline (line, 99);       
       plist = &(In.tetrahedronlist[i * 4]);

       for(j=0;j<N_CVert;j++)
        {
         dat >> k;
         plist[j] = k -1;       
        }//  for(j=0;j<N_CV
        
        
       dat >> In.tetrahedronattributelist[i];   
//        In.tetrahedronvolumelist [i] = 100.;
      } // for(i=0; i<N_
     
     
     
//      cout << i<<  "  N_Cells :  " <<  N_Cells << endl;     
//      exit(0);
     
       dat.close();
//    exit(0);    
} // ReadMeditMeshWithCells 


void ReadMeditMesh(char *SMESH, tetgenio &In)
{  
  int i, j, k, dimension, N_FVert, N_Faces, N_Vertices;
  int BDComp_Min=10000000;
  
  char line[100];
  tetgenio::facet *F;
  tetgenio::polygon *P;    

 std::ifstream dat(SMESH);

  if (!dat)
  {
    cerr << "cannot open '" << SMESH << "' for input" << endl;
    exit(-1);
  }     
   
  // check the dimension
  while (!dat.eof())
  {
    dat >> line;
    
    if ( (!strcmp(line, "Dimension"))  ||  (!strcmp(line, "dimension")) ||  (!strcmp(line, "DIMENSION")))
    {
     dat.getline (line, 99);
     dat >> dimension;
     break;
    }    
    // read until end of line
    dat.getline (line, 99);   
  }
  
  if(dimension!=3)
   {
    cerr << "dimension: " << dimension << endl;
    exit(-1);
   }
  
  // find the N_Vertices
  while (!dat.eof())
  {
    dat >> line;
    
    if ( (!strcmp(line, "Vertices")) ||  (!strcmp(line, "vertices"))   ||  (!strcmp(line, "VERTICES"))   ) 
    {
     dat.getline (line, 99);
     dat >> N_Vertices;
     break;
    }    
    // read until end of line
    dat.getline (line, 99);   
  } 
  
  In.numberofpoints = N_Vertices;
  In.pointlist = new double[3*N_Vertices];  

  //read from file
   for(i=0;i<N_Vertices; i++)
    {
     dat.getline (line, 99);
     dat >> In.pointlist[3*i] >> In.pointlist[3*i+1] >> In.pointlist[3*i+2];
     //cout<< i << " vert X: " <<In.pointlist[3*i] << " vert Y: " <<In.pointlist[3*i+1] <<endl;   
    }
    
  // find the N_Triangles
  // find the N_Vertices
  while (!dat.eof())
  {
    dat >> line;
    
    if ( (!strcmp(line, "Triangles")) ||  (!strcmp(line, "triangles"))   ||  (!strcmp(line, "TRIANGLES"))   ) 
    {
     N_FVert = 3;
     dat.getline (line, 99);
     dat >> N_Faces;
     break;
    }    
    else if ( (!strcmp(line, "Quadrilaterals")) ||  (!strcmp(line, "quadrilaterals"))   ||  (!strcmp(line, "QUADRILATERALS"))   ) 
    {
     N_FVert = 4;
     dat.getline (line, 99);
     dat >> N_Faces;
     break;
    }    
   
    // read until end of line
    dat.getline (line, 99);   
  } 
   
    In.numberoffacets = N_Faces;
    In.facetlist = new tetgenio::facet[In.numberoffacets];
    In.facetmarkerlist = new int[In.numberoffacets];
    
    for(i=0;i<N_Faces; i++)
     {
      dat.getline (line, 99);       
      
      F = &In.facetlist[i];
      tetgenio::init(F);      
      F->numberofpolygons = 1;
      F->polygonlist = new tetgenio::polygon[F->numberofpolygons];
      F->numberofholes = 0;
      F->holelist = NULL;
      P = &F->polygonlist[0];
      tetgenio::init(P);
      P->numberofvertices = N_FVert;
      P->vertexlist = new int[P->numberofvertices];
      
      for(j=0;j<N_FVert;j++)
      {
        dat >> k;
        P->vertexlist[j] = k-1;  // c numbering 
      }
      
      dat >>  In.facetmarkerlist[i];        
      
      if(BDComp_Min > In.facetmarkerlist[i])
        BDComp_Min=In.facetmarkerlist[i];
      
      // cout << i<<  " P->vertexlist[j]:  " <<  In.facetmarkerlist[i]  << endl;  
     } //   for(i=0;i<N_Faces; i++)
   
    BDComp_Min--;
   
    for(i=0;i<N_Faces; i++)
      In.facetmarkerlist[i] -= BDComp_Min;
         
//         cout << i<<  " P->vertexlist[j]:  " <<  BDComp_Min << endl;     
     
     
       dat.close();
//    exit(0);    
} // ReadMeditMesh


void DeleteDomain(TDomain *&Domain)
{
  int i, j, k, N_Cells, N_RootCells;
  int CurrVertex, N_Faces, N_Vertices, ID;
  
  TBaseCell **CellTree,  **SurfCellTree, *cell;
  TGridCell **DelCell;
  TVertex **VertexDel;  
  TCollection *coll;
  
  
    Domain->GetTreeInfo(CellTree,N_RootCells);
    coll = Domain->GetCollection(It_Finest, 0);
    N_Cells = coll->GetN_Cells();  

    
 //     cout<<"N_RootCells: "<<N_RootCells<<endl;
   // remove all existing vertices and joints
    VertexDel = new TVertex*[8*N_RootCells];
    CurrVertex = 0;   
    
   for(i=0;i<N_Cells;i++)
     {
       cell = coll->GetCell(i);
       N_Faces = cell->GetN_Faces();
       N_Vertices = cell->GetN_Vertices();
       for(j=0;j<N_Faces;j++)
         {
          if(CurrVertex==0)
           {
               VertexDel[CurrVertex++] = cell->GetVertex(j);
            }
           else
            {
             ID = 0;
             for(k=0;k<CurrVertex;k++)
             if(VertexDel[k]==cell->GetVertex(j))
              {
               ID = 1; break;
              }
             if(ID!=1)
              {
               VertexDel[CurrVertex++] = cell->GetVertex(j);
              }
          } // else if(CurrVertex==0)

           ID = 0;
           for(k=0;k<CurrVertex;k++)
           if(VertexDel[k]==cell->GetVertex((j+1)%N_Vertices))
            {
             ID = 1; break;
            }
             if(ID!=1)
            {
             VertexDel[CurrVertex++] = cell->GetVertex((j+1)%N_Vertices);
            }
 
           ID = 0;
           for(k=0;k<CurrVertex;k++)
           if(VertexDel[k]==cell->GetVertex((j+2)%N_Vertices))
            {
             ID = 1; break;
            }
             if(ID!=1)
            {
             VertexDel[CurrVertex++] = cell->GetVertex((j+2)%N_Vertices);
            }
         if(N_Faces==6) // If old cell is hexahedrol
          {  
           ID = 0;
           for(k=0;k<CurrVertex;k++)
           if(VertexDel[k]==cell->GetVertex((j+4)%N_Vertices))
            {
             ID = 1; break;
            }
             if(ID!=1)
            {
             VertexDel[CurrVertex++] = cell->GetVertex((j+4)%N_Vertices);
            }        
        }
      } // for j
    } // for i

    for(i=0;i<CurrVertex;i++)
     delete VertexDel[i];
     delete [] VertexDel;
     OutPut(CurrVertex<<" vertices were deleted"<<endl);

  // remove all existing cells and joints

   for(i=0;i<N_RootCells;i++)
    delete (TGridCell*)CellTree[i];
    delete [] CellTree;
    OutPut(N_RootCells<<" cells were deleted"<<endl);    
    
    
/*    
    cout << "DeleteDomain" <<endl;
    exit(0);
    */
}


void TetrameshCreate(TDomain *&Domain)
{
  int i, j, k, l, dimension, N_Vertices;
  int N_FVert, N_Faces, *Facemarkerlist, *Facelist, N_RootCells;
  int v1, v2, v3, v4, CellMarker, RefLevel=0, *Tetrahedrals, *PointNeighb, maxEpV, *Tetrahedrals_loc;
  int a, b, c, Neib[2], Neighb_tmp, CurrNeib, len1, len2, len3, CurrComp, N_Points ;
  int N_Cells, MaxLen, jj, N_SourcePts;  
  const int *EdgeVertex, *TmpFV, *TmpLen;
  
  double X, Y, Z, DispX, DispY, DispZ;
  double Xmin = 1e10, Xmax = -1e10, Ymin = 1e10, Ymax = -1e10;  
  double Zmin = 1e10, Zmax = -1e10, T[4]={0,0,0,0}, S[4]={0,0,0,0}; 
  double StartX, StartY, StartZ, BoundX, BoundY, BoundZ; 
  
  char *MESH, line[100];
  MESH = TDatabase::ParamDB->SMESHFILE;

  TVertex **NewVertices;
  TBaseCell **CellTree, *cell;  
  TJoint *Joint;  
  TBoundComp3D *bdcomp; 
  TShapeDesc *ShapeDesc;
  TCollection *coll;
  
#ifdef _MPI  
  int rank, size;  
  MPI_Comm Comm = TDatabase::ParamDB->Comm;
    
  MPI_Comm_rank(Comm, &rank);
  MPI_Comm_size(Comm, &size);
 
#endif     
  
  DeleteDomain(Domain);
    
  //   cout << " SMESHFILE is " << SMESH << endl;
  std::ifstream dat(MESH);

  if (!dat)
  {
    cerr << "cannot open '" << MESH << "' for input" << endl;
    exit(-1);
  }     
  
    
  // check the dimension
  while (!dat.eof())
  {
    dat >> line;
    
    if ( (!strcmp(line, "Dimension"))  ||  (!strcmp(line, "dimension")) ||  (!strcmp(line, "DIMENSION")))
    {
     dat.getline (line, 99);
     dat >> dimension;
     break;
    }    
    // read until end of line
    dat.getline (line, 99);   
  }
  
  if(dimension!=3)
   {
    cerr << "dimension: " << dimension << endl;
    exit(-1);
   }
  
  // find the N_Vertices
  while (!dat.eof())
  {
    dat >> line;
    
    if ( (!strcmp(line, "Vertices")) ||  (!strcmp(line, "vertices"))   ||  (!strcmp(line, "VERTICES"))   ) 
    {
     dat.getline (line, 99);
     dat >> N_Vertices;
     break;
    }    
    // read until end of line
    dat.getline (line, 99);   
  } 
  
   // generate new vertices
   NewVertices = new TVertex*[N_Vertices];
  
  //read from file
   for(i=0;i<N_Vertices; i++)
    {
     dat.getline (line, 99);
     dat >> X >> Y >> Z;
     
     if(TDatabase::ParamDB->P13==0)
      { //Lside
        //center point is (291.977  -1.8733333333333333 159.89 ) //old, wrong
                      // 2.922130000000000e+02 -1.8733333333333333e+00 1.583800000000000e+02 1        
       DispX = 0.;
       DispY = -1.8733333333333333;
       DispZ = 0.;
       
       
      }
     else // default
      {
       //Front side 
       //center point is  side  300.565, 69.66666666666666, 157.038  
       DispX = 0.;
       DispY = 69.66666666666666;
       DispZ = 0.;       
      }  
//        if( sqrt((X -2.922130000000000e+02)*(X -291.977) + (Y+1.8733333333333333)*(Y+-1.8733333333333333) + (Z-1.583800000000000e+02)*(Z-1.583800000000000e+02)  ) < 1.e-5)
//           cout<< i << " vert X: " <<X << " vert Y: " <<Y <<" vert Z: " <<Z <<endl;      
//   98 vert X: 292.213 vert Y: -1.87333 vert Z: 158.38
     
       
        
     NewVertices[i] = new TVertex( (X - DispX), (Y - DispY), (Z- DispZ));

//       // Pt source is at (291.977 -1.87333 159.89 )
//      if(i==98 || i==123 || i==122 )      
//       cout<< i << " vert X: " <<X << " vert Y: " <<Y <<" vert Z: " <<Z <<endl; 
//           
     
      // set bounding box
      if (X > Xmax) Xmax = X;
      if (X < Xmin) Xmin = X;
      if (Y > Ymax) Ymax = Y;
      if (Y < Ymin) Ymin = Y;
      if (Z > Zmax) Zmax = Z;
      if (Z < Zmin) Zmin = Z;
   
     //cout<< i << " vert X: " <<X << " vert Y: " <<Y <<" vert Z: " <<Z <<endl;   
    } // for(i=0;i<N_Vertices; i++)
    
//     exit(0);
    
    
   // set bounding box
    StartX = Xmin;
    StartY = Ymin;
    StartZ = Zmin;
    BoundX = Xmax - Xmin;
    BoundY = Ymax - Ymin;
    BoundZ = Zmax - Zmin;


   Domain->SetBoundBox(StartX, StartY, StartZ, BoundX, BoundY, BoundZ);   
//        cout<<Xmin <<"  "<<Ymin <<"  "<<Zmin<<endl;
//        cout<<Xmax <<"  "<<Ymax <<"  "<<Zmax<<endl;    
 
//   cout<<"  N_Vertices "<<N_Vertices <<endl;
 
 
    
 
   // find the N_RootCells
  while (!dat.eof())
  {
    dat >> line;
    
    if ( (!strcmp(line, "Tetrahedron")) ||  (!strcmp(line, "tetrahedron"))   ||  (!strcmp(line, "TETRAHEDRON"))   ) 
    {
     dat.getline (line, 99);
     dat >> N_RootCells;
     break;
    }    
    // read until end of line
    dat.getline (line, 99);   
  }   
 
  // generate new cells
   CellTree = new TBaseCell*[N_RootCells];
   Tetrahedrals = new int[4*N_RootCells];
   
   for (i=0;i<N_RootCells;i++)
   {
     dat.getline (line, 99);
     dat >> v1 >> v2 >> v3 >> v4 >> CellMarker;  
     Tetrahedrals[4*i    ] = v1 -1;
     Tetrahedrals[4*i + 1] = v2 -1;
     Tetrahedrals[4*i + 2] = v3 -1;
     Tetrahedrals[4*i + 3] = v4 -1;
     
     
     CellTree[i] = new TMacroCell(TDatabase::RefDescDB[Tetrahedron], RefLevel);
     CellTree[i]->SetRegionID(CellMarker);
     CellTree[i]->SetAsLayerCell(1);
      
     CellTree[i]->SetVertex(0, NewVertices[Tetrahedrals[4*i    ]]);
     CellTree[i]->SetVertex(1, NewVertices[Tetrahedrals[4*i + 1]]);
     CellTree[i]->SetVertex(2, NewVertices[Tetrahedrals[4*i + 2]]);
     CellTree[i]->SetVertex(3, NewVertices[Tetrahedrals[4*i + 3]]);

     CellTree[i]->SetClipBoard(i);
     ((TMacroCell *) CellTree[i])->SetSubGridID(0);
   }   

   dat.close();
  
  
   Domain->SetTreeInfo(CellTree, N_RootCells);

   // initialize iterators
   TDatabase::IteratorDB[It_EQ]->SetParam(Domain);
   TDatabase::IteratorDB[It_LE]->SetParam(Domain);
   TDatabase::IteratorDB[It_Finest]->SetParam(Domain);
   TDatabase::IteratorDB[It_Between]->SetParam(Domain);
   TDatabase::IteratorDB[It_OCAF]->SetParam(Domain);

   // search neighbours
   PointNeighb = new int[N_Vertices];
#ifdef _MPI    
     if(rank==0)   
#endif       
   cout<<"Numberofpoints "<<N_Vertices<<endl;
   memset(PointNeighb, 0, N_Vertices*SizeOfInt);     
  
     for (i=0;i<4*N_RootCells;i++)
     PointNeighb[Tetrahedrals[i]]++;

   maxEpV = 0;
   for (i=0;i<N_Vertices;i++)
     if (PointNeighb[i] > maxEpV) maxEpV = PointNeighb[i];
   delete [] PointNeighb;
 
#ifdef _MPI    
     if(rank==0)   
#endif   
   cout<<"maxEpV "<< maxEpV<<endl;   
   
   PointNeighb = new int[++maxEpV * N_Vertices];

   memset(PointNeighb, 0, maxEpV*N_Vertices*SizeOfInt);

   // every vertex contains "maxEpV" columns
   // for every vertex at first colomn contains the number of cells containing this vertex
   // at further columns we set the index of corresponding cells containing this vertex
   // cout<<"maxEpV*N_Vertices "<<maxEpV*N_Vertices<<endl;
   for(i=0;i<4*N_RootCells;i++)
    {
     j = Tetrahedrals[i]*maxEpV;
     PointNeighb[j]++;
     //cout<<"j + PointNeighb[j] " << j <<endl;
     PointNeighb[j + PointNeighb[j]] = i / 4;
    }
 

   coll=Domain->GetCollection(It_Finest, 0);
   N_Cells = coll->GetN_Cells();       
     
    for(i=0;i<N_Cells;i++)
     {
      cell = coll->GetCell(i);
      ShapeDesc= cell->GetShapeDesc();   
      ShapeDesc->GetFaceVertex(TmpFV, TmpLen, MaxLen);      
      N_Faces = cell->GetN_Faces();
      Tetrahedrals_loc = Tetrahedrals+4*i;
      
      
       for(jj=0;jj<N_Faces;jj++)
         if(cell->GetJoint(jj) == NULL)
          {
           N_Points = TmpLen[jj];
   
           if(N_Points!=3)
            {     
             cerr << "Only Tria faces are allowed!!! N_FVert: " << N_Points << endl;
             exit(-1);     
            }
            
           //printf(" TetrameshGen  Null joint  \n" );             
           a = Tetrahedrals_loc[TmpFV[jj*MaxLen]];
           b = Tetrahedrals_loc[TmpFV[jj*MaxLen + 1]];
           c = Tetrahedrals_loc[TmpFV[jj*MaxLen + 2]];

           //cout<<"  "<< a<<"  "<< b<<"  "<< c<<endl;

           Neib[0] = -1;
           Neib[1] = -1;
           CurrNeib = 0;

           len1 = PointNeighb[a*maxEpV];
           len2 = PointNeighb[b*maxEpV];
           len3 = PointNeighb[c*maxEpV]; 

           // find the index of the cells containing current face with point indices a,b,c
           for (j=1;j<=len1;j++)
            {
             Neighb_tmp = PointNeighb[a*maxEpV + j];
             for (k=1;k<=len2;k++)
              {
               if(Neighb_tmp == PointNeighb[b*maxEpV + k])
                {
                 for(l=1;l<=len3;l++)
                  if(Neighb_tmp == PointNeighb[c*maxEpV + l])
                   {
                    Neib[CurrNeib++] = Neighb_tmp;
                    break;
                   }
                 }
               }
             if (CurrNeib == 2) break;
            } // for (j=1;j<=len1;j++)
     
        if(CurrNeib==1)
         {
 
           bdcomp = Domain->GetBdPart(0)->GetBdComp(CurrComp);
       

           if(bdcomp->GetTSofXYZ(NewVertices[a]->GetX(), NewVertices[a]->GetY(),
                             NewVertices[a]->GetY(), T[1], S[1]) ||
              bdcomp->GetTSofXYZ(NewVertices[b]->GetX(), NewVertices[b]->GetY(),
                             NewVertices[b]->GetY(), T[2], S[2]) ||
              bdcomp->GetTSofXYZ(NewVertices[c]->GetX(), NewVertices[c]->GetY(),
                             NewVertices[c]->GetY(), T[3], S[3])    )
             {
              cerr<<"Error: could not set parameter values"<<endl;
              OutPut(NewVertices[a]<<endl);
              OutPut(NewVertices[b]<<endl);
              OutPut(NewVertices[c]<<endl);
              exit(0);
             }
 
//        if (CurrNeib == 2)
//         {
//          if(bdcomp->IsFreeBoundary())
//            Joint = new TIsoInterfaceJoint3D(bdcomp, T, S,
//                        CellTree[Neib[0]], CellTree[Neib[1]] );
//           else
//            Joint = new TInterfaceJoint3D(bdcomp, T, S, 
//                        CellTree[Neib[0]], CellTree[Neib[1]] );
//         }
//        else
//         {
//          if(bdcomp->IsFreeBoundary())
//            Joint = new TIsoBoundFace(bdcomp, T, S);
//          else
           Joint = new TBoundFace(bdcomp, T, S);
//         }
         }           
        else
        {
         // cout<<"Inner face"<<endl;
         if (CurrNeib != 2)
         cerr << "Error !!!!!!!! not enough neighbours!" << endl;

         Joint = new TJointEqN(CellTree[Neib[0]], CellTree[Neib[1]]);
        }
 
        // First element containing the current face
        // find the local index for the point 'a' on the cell
        for (j=0;j<4;j++)
         if (Tetrahedrals[4*Neib[0]+j] == a) break;

        // find the local index for the point 'b' on the cell
        for (k=0;k<4;k++)
         if (Tetrahedrals[4*Neib[0]+k] == b) break;

        // find the local index for the point 'c' on the cell
        for (l=0;l<4;l++)
         if (Tetrahedrals[4*Neib[0]+l] == c) break;   

        l = l*100 + k*10 + j;  

        //cout<<" l "<< l <<endl;           

     switch (l) // j will contain the local index for the current face
      {
        case 210: case 21: case 102:
        case 120: case 12: case 201:
          j = 0;
          break;  
        case 310: case 31: case 103:
        case 130: case 13: case 301:
          j = 1;
          break;  
        case 321: case 132: case 213:
        case 231: case 123: case 312:
          j = 2;
          break;  
        case 230: case 23: case 302:
        case 320: case 32: case 203:
          j = 3;
          break; 

       default:
        Error("Unable to set the face !!!!!!!!!!!!" << endl);
        exit(0);
      } // switch (l)
     
     CellTree[Neib[0]]->SetJoint(j, Joint);

     // second cell containing the current face 
    if (Neib[1] != -1) // second element containing the current face
     {     
     // find the local index for the point 'a' on the cell
     for (j=0;j<4;j++)
      if (Tetrahedrals[4*Neib[1]+j] == a) break;

     // find the local index for the point 'b' on the cell
     for (k=0;k<4;k++)
      if (Tetrahedrals[4*Neib[1]+k] == b) break;

      // find the local index for the point 'c' on the cell
     for (l=0;l<4;l++)
      if (Tetrahedrals[4*Neib[1]+l] == c) break;   

     l = l*100 + k*10 + j;  

     //cout<<" l "<< l <<endl;

     switch (l) // j will contain the local index for the current face
      {
        case 210: case 21: case 102:
        case 120: case 12: case 201:
          j = 0;
          break;  
        case 310: case 31: case 103:
        case 130: case 13: case 301:
          j = 1;
          break;  
        case 321: case 132: case 213:
        case 231: case 123: case 312:
          j = 2;
          break;  
        case 230: case 23: case 302:
        case 320: case 32: case 203:
          j = 3;
          break; 

      default:
       Error("Unable to set the face !!!!!!!!!!!!" << endl);
       exit(0);
       }
       CellTree[Neib[1]]->SetJoint(j, Joint);   
     } // if (Neib[1] != -1) 
     
     
      if (Joint->GetType() == InterfaceJoint3D ||
       Joint->GetType() == IsoInterfaceJoint3D)
        {
         ((TInterfaceJoint3D*)Joint)->SetMapType();
         ((TInterfaceJoint3D*)(Joint))->CheckOrientation();
        }
       else 
        if (Joint->GetType() == JointEqN)
           Joint->SetMapType();      
       
       } //  if(cell->GetJoint(j)
      }  // for(i=0;i<N_Cells;i++)
     
delete [] PointNeighb;

// #ifdef _MPI    
//      if(rank==1)
//     printf(" TetrameshGen  complete \n" );  
//     MPI_Finalize(); 
// #endif    
//   exit(0);
// //       
  
}


void TetrameshGen(TDomain *&Domain)
{
 //======================================================================
 // Tetgen for grid generation begin
 //======================================================================
  int i, j, k, l, N_Coord, *N_FVerts, N_Faces, *Facets;
  int N, N_RootCells, N_Cells, CurrVertex, N_Vertices, ID, N_G, RefLevel=0;
  int CurrNeib, len1, len2, len3, maxEpV = 0, a, b, c, Neib[2], Neighb_tmp, CurrComp;
  int *Tetrahedrals, *PointNeighb, dimension;
  int *Facelist, *Facemarkerlist;
 
  double *Coordinates, N_x, N_y, N_z;   
  double *Vertices;
  double Xmin = 1e10, Xmax = -1e10, Ymin = 1e10, Ymax = -1e10;  
  double Zmin = 1e10, Zmax = -1e10, T[4]={0,0,0,0}, S[4]={0,0,0,0}; 
  double StartX, StartY, StartZ, BoundX, BoundY, BoundZ; 
  
  tetgenio In, Out;
   
  TBaseCell **CellTree,  **SurfCellTree;
  TGridCell **DelCell;
  TVertex **VertexDel, **NewVertices, **NewSurfVertices; 
  TBoundPart *BoundPart;
  TBdPlane **UpdateFaceParams;
  TJoint *Joint;
  TBoundComp3D *bdcomp; 
  TCollection *coll, *SurfColl;
  TBaseCell *cell;
  TBdSphere *UpdateParam;  

  char *SMESH, line[100];
  SMESH = TDatabase::ParamDB->SMESHFILE;

#ifdef _MPI  
  int rank, size;  
  MPI_Comm Comm = TDatabase::ParamDB->Comm;
    
  MPI_Comm_rank(Comm, &rank);
  MPI_Comm_size(Comm, &size);
 
#endif    

  

#ifdef _MPI    
  if(rank==0)
#endif  
   {
   
//     In.initialize();
//     Out.initialize();
   
    std::ostringstream opts;
    opts << " "; 
  
    opts.seekp(std::ios::beg);
    opts<<'p'; // Tetrahedralize the PLC. Switches are chosen to read a PLC (p)
//     opts<<'r'; // -r  Reconstructs a previously generated mesh     
//     opts<<"q"<<1.49; // quality mesh generation(q) with a specified quality bound
//     opts<<"a"<<0.1; // maximum volume constraint
//     opts<<'i'; // Inserts a list of additional points into mesh.
    opts<<'z'; // numbers all output items starting from zero
//     opts<<'d'; // Detect intersections of PLC facets.
    opts<<'f'; // Outputs all  faces (including non-boundary) 
    opts<<'e'; // Outputs a list of edges of the triangulation
//     opts<<'I'; // Suppresses mesh iteration numbers.
    opts<<'C'; // Checks the consistency of the final mesh.
//     opts<<'Q'; // Quiet: No terminal output except errors.
//     opts<<'g'; // Outputs mesh to .mesh file for viewing by Medit
    opts<<'Y'; // Suppresses boundary facets/segments splitting
//     opts<<'V';  //verbose mode
    opts<<ends;
    
  
//   cout << " SMESHFILE is " << SMESH << endl;
  ReadMeditMesh(SMESH, In);
  
  
//     for(i=0;i<In.numberofpoints;i++)
//       OutPut(i<<" (x, y, z) =  "<<
//        In.pointlist[3*i]<<' '<<In.pointlist[3*i+1]<<' '<<In.pointlist[3*i+2]<<endl);

 // Calling  tetrahedralize function of 3dtetgen mesh generator
    tetrahedralize((char*)opts.str().c_str(), &In, &Out);
    
   } // if(rank==0)
     
    exit(0);
 //   output: coordinates of all vertices
//  for(i=0;i<Out.numberofpoints;i++)
//   OutPut(" (x, y, z) =  "<< Out.pointlist[3*i]<<' '<<Out.pointlist[3*i+1]<<' '<<Out.pointlist[3*i+2]<<endl);

    Domain->GetTreeInfo(CellTree,N_RootCells);
    coll = Domain->GetCollection(It_Finest, 0);
    N_Cells = coll->GetN_Cells();

//     cout<<"N_RootCells: "<<N_RootCells<<endl;
   // remove all existing vertices and joints
    VertexDel = new TVertex*[8*N_RootCells];
    CurrVertex = 0;

   for(i=0;i<N_Cells;i++)
     {
       cell = coll->GetCell(i);
       N_Faces = cell->GetN_Faces();
       N_Vertices = cell->GetN_Vertices();
       for(j=0;j<N_Faces;j++)
         {
          if(CurrVertex==0)
           {
               VertexDel[CurrVertex++] = cell->GetVertex(j);
            }
           else
            {
             ID = 0;
             for(k=0;k<CurrVertex;k++)
             if(VertexDel[k]==cell->GetVertex(j))
              {
               ID = 1; break;
              }
             if(ID!=1)
              {
               VertexDel[CurrVertex++] = cell->GetVertex(j);
              }
          } // else if(CurrVertex==0)

           ID = 0;
           for(k=0;k<CurrVertex;k++)
           if(VertexDel[k]==cell->GetVertex((j+1)%N_Vertices))
            {
             ID = 1; break;
            }
             if(ID!=1)
            {
             VertexDel[CurrVertex++] = cell->GetVertex((j+1)%N_Vertices);
            }
 
           ID = 0;
           for(k=0;k<CurrVertex;k++)
           if(VertexDel[k]==cell->GetVertex((j+2)%N_Vertices))
            {
             ID = 1; break;
            }
             if(ID!=1)
            {
             VertexDel[CurrVertex++] = cell->GetVertex((j+2)%N_Vertices);
            }
         if(N_Faces==6) // If old cell is hexahedrol
          {  
           ID = 0;
           for(k=0;k<CurrVertex;k++)
           if(VertexDel[k]==cell->GetVertex((j+4)%N_Vertices))
            {
             ID = 1; break;
            }
             if(ID!=1)
            {
             VertexDel[CurrVertex++] = cell->GetVertex((j+4)%N_Vertices);
            }        
        }
      } // for j
    } // for i

    for(i=0;i<CurrVertex;i++)
     delete VertexDel[i];
     delete [] VertexDel;
     OutPut(CurrVertex<<" vertices were deleted"<<endl);

  // remove all existing cells and joints

   for(i=0;i<N_RootCells;i++)
    delete (TGridCell*)CellTree[i];
    delete [] CellTree;
    OutPut(N_RootCells<<" cells were deleted"<<endl);

   
#ifdef _MPI     
     if(rank==0) 
#endif  
     {
       N_RootCells = Out.numberoftetrahedra;
       N_G = Out.numberofpoints;   

       // allocate auxillary fields
       Coordinates = Out.pointlist;  
       Tetrahedrals = Out.tetrahedronlist;       
     }

     
     
#ifdef _MPI
   MPI_Bcast(&N_RootCells, 1, MPI_INT, 0, Comm);
   MPI_Bcast(&N_G, 1, MPI_INT, 0, Comm);  
   
   
  if(rank!=0)
   {
    Coordinates = new double [3*N_G];
    Tetrahedrals = new int[4*N_RootCells];
   }
   
   MPI_Bcast(Coordinates, 3*N_G, MPI_DOUBLE, 0, Comm);  
   MPI_Bcast(Tetrahedrals, 4*N_RootCells, MPI_INT, 0, Comm);     
#endif      
   
  // generate new vertices
   NewVertices = new TVertex*[N_G];

   for (i=0;i<N_G;i++)
    {
      NewVertices[i] = new TVertex(Coordinates[3*i], Coordinates[3*i+1], Coordinates[3*i+2]);

      // set bounding box
      if (Coordinates[3*i] > Xmax) Xmax = Coordinates[3*i];
      if (Coordinates[3*i] < Xmin) Xmin = Coordinates[3*i];
      if (Coordinates[3*i+1] > Ymax) Ymax = Coordinates[3*i+1];
      if (Coordinates[3*i+1] < Ymin) Ymin = Coordinates[3*i+1];
      if (Coordinates[3*i+2] > Zmax) Zmax = Coordinates[3*i+2];
      if (Coordinates[3*i+2] < Zmin) Zmin = Coordinates[3*i+2];
   }

   // set bounding box
    StartX = Xmin;
    StartY = Ymin;
    StartZ = Zmin;
    BoundX = Xmax - Xmin;
    BoundY = Ymax - Ymin;
    BoundZ = Zmax - Zmin;


   Domain->SetBoundBox(StartX, StartY, StartZ, BoundX, BoundY, BoundZ);   
//        cout<<Xmin <<"  "<<Ymin <<"  "<<Zmin<<endl;
//        cout<<Xmax <<"  "<<Ymax <<"  "<<Zmax<<endl;

   CellTree = new TBaseCell*[N_RootCells];
 //  output of each tetraheron vertex indices (four vertices for each)
//   for (i=0;i<N_RootCells;i++)
//        cout<< Tetrahedrals[4*i]<<"  "<<Tetrahedrals[4*i + 1]<<"  "
//          <<Tetrahedrals[4*i + 2]<<"  "<<Tetrahedrals[4*i + 3]<<endl;


   for (i=0;i<N_RootCells;i++)
   {
     CellTree[i] = new TMacroCell(TDatabase::RefDescDB[Tetrahedron],
                                    RefLevel);

     CellTree[i]->SetVertex(0, NewVertices[Tetrahedrals[4*i    ]]);
     CellTree[i]->SetVertex(1, NewVertices[Tetrahedrals[4*i + 1]]);
     CellTree[i]->SetVertex(2, NewVertices[Tetrahedrals[4*i + 2]]);
     CellTree[i]->SetVertex(3, NewVertices[Tetrahedrals[4*i + 3]]);

    CellTree[i]->SetClipBoard(i);
     ((TMacroCell *) CellTree[i])->SetSubGridID(0);
   }

   Domain->SetTreeInfo(CellTree, N_RootCells);


   // initialize iterators
   TDatabase::IteratorDB[It_EQ]->SetParam(Domain);
   TDatabase::IteratorDB[It_LE]->SetParam(Domain);
   TDatabase::IteratorDB[It_Finest]->SetParam(Domain);
   TDatabase::IteratorDB[It_Between]->SetParam(Domain);
   TDatabase::IteratorDB[It_OCAF]->SetParam(Domain);
   
  
   // search neighbours
   PointNeighb = new int[N_G];
#ifdef _MPI    
     if(rank==0)   
#endif       
   cout<<"numberofpoints "<<N_G<<endl;
   memset(PointNeighb, 0, N_G *SizeOfInt);

     for (i=0;i<4*N_RootCells;i++)
     PointNeighb[Tetrahedrals[i]]++;

   for (i=0;i<N_G;i++)
     if (PointNeighb[i] > maxEpV) maxEpV = PointNeighb[i];
   delete [] PointNeighb;

#ifdef _MPI    
     if(rank==0)   
#endif   
   cout<<"maxEpV "<< maxEpV<<endl;

   PointNeighb = new int[++maxEpV * N_G];

   memset(PointNeighb, 0, maxEpV*N_G*SizeOfInt);

    // every vertex contains "maxEpV" columns
    // for every vertex at first colomn contains the number of cells containing this vertex
    // at further columns we set the index of corresponding cells containing this vertex
    // cout<<"maxEpV*N_G "<<maxEpV*N_G<<endl;

   for(i=0;i<4*N_RootCells;i++)
    {
     j = Tetrahedrals[i]*maxEpV;
     PointNeighb[j]++;
     //cout<<"j + PointNeighb[j] " << j <<endl;
     PointNeighb[j + PointNeighb[j]] = i / 4;
    }
 //  output of PointNeighb columns for each point
//   for (i=0;i<N_G;i++)
//    {
//     for (j=0;j<maxEpV;j++)
//     cout<<"  "<< PointNeighb[i*maxEpV+j];
//     cout<<endl;
//    }



   // generate new faces 
  
#ifdef _MPI     
     if(rank==0) 
#endif  
     {    
      N_G =  Out.numberoftrifaces;
      Facelist = Out.trifacelist;
      Facemarkerlist = Out.trifacemarkerlist;
     }
      
#ifdef _MPI
   MPI_Bcast(&N_G, 1, MPI_INT, 0, Comm);  

  if(rank!=0)
   {
    Facelist = new int [3*N_G];
    Facemarkerlist = new int[N_G];
   }   
   
   
   MPI_Bcast(Facelist, 3*N_G, MPI_INT, 0, Comm); 
   MPI_Bcast(Facemarkerlist, N_G, MPI_INT, 0, Comm); 
  
  if(rank==0) 
#endif      
   cout<<"numberoftrifaces "<<N_G<<endl;
    
   for (i=0;i<N_G;i++)
   {
     a = Facelist[3*i];
     b = Facelist[3*i+1];
     c = Facelist[3*i+2];

//      cout<<"  "<< a<<"  "<< b<<"  "<< c<<endl;

     Neib[0] = -1;
     Neib[1] = -1;
     CurrNeib = 0;

     len1 = PointNeighb[a*maxEpV];
     len2 = PointNeighb[b*maxEpV];
     len3 = PointNeighb[c*maxEpV];

   // find the index of the cells containing current face with point indices a,b,c
    for (j=1;j<=len1;j++)
     {
       Neighb_tmp = PointNeighb[a*maxEpV + j];
        for (k=1;k<=len2;k++)
         {
          if (Neighb_tmp == PointNeighb[b*maxEpV + k])
           {
            for (l=1;l<=len3;l++)
             if (Neighb_tmp == PointNeighb[c*maxEpV + l])
             {
              Neib[CurrNeib++] = Neighb_tmp;
              break;
             }
           }
          }
       if (CurrNeib == 2) break;
     }
 //   cout<<"CurrNeib " << CurrNeib <<endl;
// cout<<"Facemarkerlist[i] : "<<Facemarkerlist[i]<<endl;
     if (Facemarkerlist[i]) // 0 for inner edges and Boundcomp+1 for Boundedge respect
      {

       CurrComp = Facemarkerlist[i] - 1;
//        cout<<"Boundary face CurrComp: "<<CurrComp<<endl;
// exit(0);
       CurrComp = 0; // not yet implemented fully
       
       
       
       bdcomp = Domain->GetBdPart(0)->GetBdComp(CurrComp);
       

       if(bdcomp->GetTSofXYZ(NewVertices[a]->GetX(), NewVertices[a]->GetY(),
                             NewVertices[a]->GetY(), T[1], S[1]) ||
          bdcomp->GetTSofXYZ(NewVertices[b]->GetX(), NewVertices[b]->GetY(),
                             NewVertices[b]->GetY(), T[2], S[2]) ||
          bdcomp->GetTSofXYZ(NewVertices[c]->GetX(), NewVertices[c]->GetY(),
                             NewVertices[c]->GetY(), T[3], S[3])    )
         {
          cerr<<"Error: could not set parameter values"<<endl;
          OutPut(NewVertices[a]<<endl);
          OutPut(NewVertices[b]<<endl);
          OutPut(NewVertices[c]<<endl);
          exit(0);
         }
 
       if (CurrNeib == 2)
        {
         if(bdcomp->IsFreeBoundary())
           Joint = new TIsoInterfaceJoint3D(bdcomp, T, S,
                       CellTree[Neib[0]], CellTree[Neib[1]] );
          else
           Joint = new TInterfaceJoint3D(bdcomp, T, S, 
                       CellTree[Neib[0]], CellTree[Neib[1]] );
        }
       else
        {
         if(bdcomp->IsFreeBoundary())
           Joint = new TIsoBoundFace(bdcomp, T, S);
         else
           Joint = new TBoundFace(bdcomp, T, S);
        }
      }
     else
      {
// //       cout<<"Inner face"<<endl;
       if (CurrNeib != 2)
       cerr << "Error !!!!!!!! not enough neighbours!" << endl;

       Joint = new TJointEqN(CellTree[Neib[0]], CellTree[Neib[1]]);

      }

    // First element containing the current face
    // find the local index for the point 'a' on the cell
    for (j=0;j<4;j++)
      if (Tetrahedrals[4*Neib[0]+j] == a) break;

    // find the local index for the point 'b' on the cell
    for (k=0;k<4;k++)
      if (Tetrahedrals[4*Neib[0]+k] == b) break;

       // find the local index for the point 'c' on the cell
    for (l=0;l<4;l++)
      if (Tetrahedrals[4*Neib[0]+l] == c) break;   

     l = l*100 + k*10 + j;  

//      cout<<""<< l <<endl;

     switch (l) // j will contain the local index for the current face
      {
        case 210: case 21: case 102:
        case 120: case 12: case 201:
          j = 0;
          break;  
        case 310: case 31: case 103:
        case 130: case 13: case 301:
          j = 1;
          break;  
        case 321: case 132: case 213:
        case 231: case 123: case 312:
          j = 2;
          break;  
        case 230: case 23: case 302:
        case 320: case 32: case 203:
          j = 3;
          break; 

      default:
       Error("Unable to set the face !!!!!!!!!!!!" << endl);
       exit(0);
     }
      CellTree[Neib[0]]->SetJoint(j, Joint);

   if (Neib[1] != -1) // second element containing the current face
    {
          // find the local index for the point 'a' on the cell
    for (j=0;j<4;j++)
      if (Tetrahedrals[4*Neib[1]+j] == a) break;

    // find the local index for the point 'b' on the cell
    for (k=0;k<4;k++)
      if (Tetrahedrals[4*Neib[1]+k] == b) break;

       // find the local index for the point 'c' on the cell
    for (l=0;l<4;l++)
      if (Tetrahedrals[4*Neib[1]+l] == c) break;   

     l = l*100 + k*10 + j;  

//      cout<<""<< l <<endl;

     switch (l) // j will contain the local index for the current face
      {
        case 210: case 21: case 102:
        case 120: case 12: case 201:
          j = 0;
          break;  
        case 310: case 31: case 103:
        case 130: case 13: case 301:
          j = 1;
          break;  
        case 321: case 132: case 213:
        case 231: case 123: case 312:
          j = 2;
          break;  
        case 230: case 23: case 302:
        case 320: case 32: case 203:
          j = 3;
          break; 

      default:
       Error("Unable to set the face !!!!!!!!!!!!" << endl);
       exit(0);
      }
      CellTree[Neib[1]]->SetJoint(j, Joint);
     }

  if (Joint->GetType() == InterfaceJoint3D ||
      Joint->GetType() == IsoInterfaceJoint3D)
      {
        ((TInterfaceJoint3D*)Joint)->SetMapType();
        ((TInterfaceJoint3D*)(Joint))->CheckOrientation();
      }
      else 
       if (Joint->GetType() == JointEqN)
           Joint->SetMapType();

  }

#ifdef _MPI    
  if(rank==0)    
#endif
  {
//    In.deinitialize();
//    Out.deinitialize();
  }
  

 
 
#ifdef _MPI    
  if(rank!=0) 
   {
    delete [] Tetrahedrals ;  
    delete [] Coordinates;      
    delete [] Facelist;
    delete [] Facemarkerlist;
   }
#endif  
  
  delete [] PointNeighb;

// #ifdef _MPI    
//      if(rank==3)
//     printf(" TetrameshGen  %d \n",   N_G );  
//     MPI_Finalize(); 
// #endif    
//   exit(0);
//      

}


void TetraGen(TDomain *&Domain)
{
 //======================================================================
 // Tetgen for grid generation begin
 //======================================================================
  int i, j, k, l, N_Coord, *N_FVerts, N_Faces, *Facets;
  int N, N_RootCells, N_Cells, CurrVertex, N_Vertices, ID, N_G, RefLevel=0;
  int CurrNeib, len1, len2, len3, maxEpV = 0, a, b, c, Neib[2], Neighb_tmp, CurrComp;
  int *Tetrahedrals, *PointNeighb, dimension;
  int jj, N_Points, MaxLen, *Tetrahedrals_loc;
  const int *EdgeVertex, *TmpFV, *TmpLen;

  double *Coordinates, N_x, N_y, N_z;   
  double *Vertices;
  double Xmin = 1e10, Xmax = -1e10, Ymin = 1e10, Ymax = -1e10;  
  double Zmin = 1e10, Zmax = -1e10, T[4]={0,0,0,0}, S[4]={0,0,0,0}; 
  double StartX, StartY, StartZ, BoundX, BoundY, BoundZ; 
  
  tetgenio In, Out;
   
  TBaseCell **CellTree,  **SurfCellTree;
  TGridCell **DelCell;
  TVertex **VertexDel, **NewVertices, **NewSurfVertices; 
  TBoundPart *BoundPart;
  TBdPlane **UpdateFaceParams;
  TJoint *Joint;
  TBoundComp3D *bdcomp; 
  TCollection *coll, *SurfColl;
  TBaseCell *cell;
  TBdSphere *UpdateParam;  
  TShapeDesc *ShapeDesc;
   
  char *SMESH, line[100];
  SMESH = TDatabase::ParamDB->SMESHFILE;

#ifdef _MPI  
  int rank, size;  
  MPI_Comm Comm = TDatabase::ParamDB->Comm;
    
  MPI_Comm_rank(Comm, &rank);
  MPI_Comm_size(Comm, &size);
 
#endif    

#ifdef _MPI    
  if(rank==0)
#endif  
   {
   
//     In.initialize();
//     Out.initialize();
   
    std::ostringstream opts;
    opts << " "; 
  
    opts.seekp(std::ios::beg);
    opts<<'p'; // Tetrahedralize the PLC. Switches are chosen to read a PLC (p)
//     opts<<'r'; // -r  Reconstructs a previously generated mesh     
    opts<<"q"<<1.20; // quality mesh generation(q) with a specified quality bound
//     opts<<"a"<<0.1; // maximum volume constraint
//     opts<<'i'; // Inserts a list of additional points into mesh.
    opts<<'z'; // numbers all output items starting from zero
//     opts<<'d'; // Detect intersections of PLC facets.
    opts<<'f'; // Outputs all  faces (including non-boundary) 
    opts<<'e'; // Outputs a list of edges of the triangulation
//     opts<<'I'; // Suppresses mesh iteration numbers.
    opts<<'C'; // Checks the consistency of the final mesh.
//     opts<<'Q'; // Quiet: No terminal output except errors.
//     opts<<'g'; // Outputs mesh to .mesh file for viewing by Medit
    opts<<'Y'; // Suppresses boundary facets/segments splitting
//     opts<<'V';  //verbose mode
    opts<<ends;
    
  
//   cout << " SMESHFILE is " << SMESH << endl;
  ReadMeditMesh(SMESH, In);
  
  
//     for(i=0;i<In.numberofpoints;i++)
//       OutPut(i<<" (x, y, z) =  "<<
//        In.pointlist[3*i]<<' '<<In.pointlist[3*i+1]<<' '<<In.pointlist[3*i+2]<<endl);

 // Calling  tetrahedralize function of 3dtetgen mesh generator
    tetrahedralize((char*)opts.str().c_str(), &In, &Out);

 //   output: coordinates of all vertices
//  for(i=0;i<Out.numberofpoints;i++)
//   OutPut(" (x, y, z) =  "<< Out.pointlist[3*i]<<' '<<Out.pointlist[3*i+1]<<' '<<Out.pointlist[3*i+2]<<endl);    
    
   } // if(rank==0)
     
  DeleteDomain(Domain);

   
#ifdef _MPI     
     if(rank==0) 
#endif  
     {
       N_RootCells = Out.numberoftetrahedra;
       N_G = Out.numberofpoints;   

       // allocate auxillary fields
       Coordinates = Out.pointlist;  
       Tetrahedrals = Out.tetrahedronlist;       
     }
 
#ifdef _MPI
   MPI_Bcast(&N_RootCells, 1, MPI_INT, 0, Comm);
   MPI_Bcast(&N_G, 1, MPI_INT, 0, Comm);  
   
   
  if(rank!=0)
   {
    Coordinates = new double [3*N_G];
    Tetrahedrals = new int[4*N_RootCells];
   }
   
   MPI_Bcast(Coordinates, 3*N_G, MPI_DOUBLE, 0, Comm);  
   MPI_Bcast(Tetrahedrals, 4*N_RootCells, MPI_INT, 0, Comm);     
#endif      
   
  // generate new vertices
   NewVertices = new TVertex*[N_G];

   for (i=0;i<N_G;i++)
    {
      NewVertices[i] = new TVertex(Coordinates[3*i], Coordinates[3*i+1], Coordinates[3*i+2]);

      // set bounding box
      if (Coordinates[3*i] > Xmax) Xmax = Coordinates[3*i];
      if (Coordinates[3*i] < Xmin) Xmin = Coordinates[3*i];
      if (Coordinates[3*i+1] > Ymax) Ymax = Coordinates[3*i+1];
      if (Coordinates[3*i+1] < Ymin) Ymin = Coordinates[3*i+1];
      if (Coordinates[3*i+2] > Zmax) Zmax = Coordinates[3*i+2];
      if (Coordinates[3*i+2] < Zmin) Zmin = Coordinates[3*i+2];
    }

   // set bounding box
    StartX = Xmin;
    StartY = Ymin;
    StartZ = Zmin;
    BoundX = Xmax - Xmin;
    BoundY = Ymax - Ymin;
    BoundZ = Zmax - Zmin;


   Domain->SetBoundBox(StartX, StartY, StartZ, BoundX, BoundY, BoundZ);   
//        cout<<Xmin <<"  "<<Ymin <<"  "<<Zmin<<endl;
//        cout<<Xmax <<"  "<<Ymax <<"  "<<Zmax<<endl;

   CellTree = new TBaseCell*[N_RootCells];
 //  output of each tetraheron vertex indices (four vertices for each)
//   for (i=0;i<N_RootCells;i++)
//        cout<< Tetrahedrals[4*i]<<"  "<<Tetrahedrals[4*i + 1]<<"  "
//          <<Tetrahedrals[4*i + 2]<<"  "<<Tetrahedrals[4*i + 3]<<endl;


   for (i=0;i<N_RootCells;i++)
    {
     CellTree[i] = new TMacroCell(TDatabase::RefDescDB[Tetrahedron],
                                    RefLevel);

     CellTree[i]->SetVertex(0, NewVertices[Tetrahedrals[4*i    ]]);
     CellTree[i]->SetVertex(1, NewVertices[Tetrahedrals[4*i + 1]]);
     CellTree[i]->SetVertex(2, NewVertices[Tetrahedrals[4*i + 2]]);
     CellTree[i]->SetVertex(3, NewVertices[Tetrahedrals[4*i + 3]]);

     CellTree[i]->SetClipBoard(i);
      ((TMacroCell *) CellTree[i])->SetSubGridID(0);
    }

   Domain->SetTreeInfo(CellTree, N_RootCells);


   // initialize iterators
   TDatabase::IteratorDB[It_EQ]->SetParam(Domain);
   TDatabase::IteratorDB[It_LE]->SetParam(Domain);
   TDatabase::IteratorDB[It_Finest]->SetParam(Domain);
   TDatabase::IteratorDB[It_Between]->SetParam(Domain);
   TDatabase::IteratorDB[It_OCAF]->SetParam(Domain);
   
  
   // search neighbours
   PointNeighb = new int[N_G];
#ifdef _MPI    
     if(rank==0)   
#endif       
   cout<<"numberofpoints "<<N_G<<endl;
   memset(PointNeighb, 0, N_G *SizeOfInt);

     for (i=0;i<4*N_RootCells;i++)
     PointNeighb[Tetrahedrals[i]]++;

   for (i=0;i<N_G;i++)
     if (PointNeighb[i] > maxEpV) maxEpV = PointNeighb[i];
   delete [] PointNeighb;

#ifdef _MPI    
     if(rank==0)   
#endif   
   cout<<"maxEpV "<< maxEpV<<endl;

   PointNeighb = new int[++maxEpV * N_G];

   memset(PointNeighb, 0, maxEpV*N_G*SizeOfInt);

    // every vertex contains "maxEpV" columns
    // for every vertex at first colomn contains the number of cells containing this vertex
    // at further columns we set the index of corresponding cells containing this vertex
    // cout<<"maxEpV*N_G "<<maxEpV*N_G<<endl;

   for(i=0;i<4*N_RootCells;i++)
    {
     j = Tetrahedrals[i]*maxEpV;
     PointNeighb[j]++;
     //cout<<"j + PointNeighb[j] " << j <<endl;
     PointNeighb[j + PointNeighb[j]] = i / 4;
    }


   // generate new faces    
   coll=Domain->GetCollection(It_Finest, 0);
   N_Cells = coll->GetN_Cells();       
   N_G = 0;
   for (i=0;i<N_Cells;i++)
   {
    cell = coll->GetCell(i);
    ShapeDesc= cell->GetShapeDesc();   
    ShapeDesc->GetFaceVertex(TmpFV, TmpLen, MaxLen);      
    N_Faces = cell->GetN_Faces();
    Tetrahedrals_loc = Tetrahedrals+4*i;
           
    for(jj=0;jj<N_Faces;jj++)
     if(cell->GetJoint(jj) == NULL)
      {
       N_G++;    
       N_Points = TmpLen[jj];
   
       if(N_Points!=3)
        {     
         cerr << "Only Tria faces are allowed!!! N_FVert: " << N_Points << endl;
         exit(-1);     
        }     
     
     
       //printf(" TetrameshGen  Null joint  \n" );             
       a = Tetrahedrals_loc[TmpFV[jj*MaxLen]];
       b = Tetrahedrals_loc[TmpFV[jj*MaxLen + 1]];
       c = Tetrahedrals_loc[TmpFV[jj*MaxLen + 2]];

//      cout<<"  "<< a<<"  "<< b<<"  "<< c<<endl;

       Neib[0] = -1;
       Neib[1] = -1;
       CurrNeib = 0;

       len1 = PointNeighb[a*maxEpV];
       len2 = PointNeighb[b*maxEpV];
       len3 = PointNeighb[c*maxEpV];

       // find the index of the cells containing current face with point indices a,b,c
       for (j=1;j<=len1;j++)
       {
        Neighb_tmp = PointNeighb[a*maxEpV + j];
        for (k=1;k<=len2;k++)
         {
          if (Neighb_tmp == PointNeighb[b*maxEpV + k])
           {
            for (l=1;l<=len3;l++)
             if (Neighb_tmp == PointNeighb[c*maxEpV + l])
             {
              Neib[CurrNeib++] = Neighb_tmp;
              break;
             }
            }
           }
         if (CurrNeib == 2) break;
        }// for (j=1;j<=len1;j++)
        
      //   cout<<"CurrNeib " << CurrNeib <<endl;
      // cout<<"Facemarkerlist[i] : "<<Facemarkerlist[i]<<endl;
      if (CurrNeib == 1) // 0 for inner edges and Boundcomp+1 for Boundedge respect
       {
        CurrComp = 0; // not yet implemented fully
       
       
       
        bdcomp = Domain->GetBdPart(0)->GetBdComp(CurrComp);
       
 
        if(bdcomp->GetTSofXYZ(NewVertices[a]->GetX(), NewVertices[a]->GetY(),
                             NewVertices[a]->GetY(), T[1], S[1]) ||
          bdcomp->GetTSofXYZ(NewVertices[b]->GetX(), NewVertices[b]->GetY(),
                             NewVertices[b]->GetY(), T[2], S[2]) ||
          bdcomp->GetTSofXYZ(NewVertices[c]->GetX(), NewVertices[c]->GetY(),
                             NewVertices[c]->GetY(), T[3], S[3])    )
         {
          cerr<<"Error: could not set parameter values"<<endl;
          OutPut(NewVertices[a]<<endl);
          OutPut(NewVertices[b]<<endl);
          OutPut(NewVertices[c]<<endl);
          exit(0);
         }
 
       if (CurrNeib == 2)
        {
         if(bdcomp->IsFreeBoundary())
           Joint = new TIsoInterfaceJoint3D(bdcomp, T, S,
                       CellTree[Neib[0]], CellTree[Neib[1]] );
          else
           Joint = new TInterfaceJoint3D(bdcomp, T, S, 
                       CellTree[Neib[0]], CellTree[Neib[1]] );
        }
       else
        {
         if(bdcomp->IsFreeBoundary())
           Joint = new TIsoBoundFace(bdcomp, T, S);
         else
           Joint = new TBoundFace(bdcomp, T, S);
        }
      }
     else
      {
       // cout<<"Inner face"<<endl;
       if (CurrNeib != 2)
       cerr << "Error !!!!!!!! not enough neighbours!" << endl;

       Joint = new TJointEqN(CellTree[Neib[0]], CellTree[Neib[1]]);
      }

     // First element containing the current face
     // find the local index for the point 'a' on the cell
     for (j=0;j<4;j++)
      if (Tetrahedrals[4*Neib[0]+j] == a) break;

     // find the local index for the point 'b' on the cell
     for (k=0;k<4;k++)
      if (Tetrahedrals[4*Neib[0]+k] == b) break;

      // find the local index for the point 'c' on the cell
     for (l=0;l<4;l++)
      if (Tetrahedrals[4*Neib[0]+l] == c) break;   

      l = l*100 + k*10 + j;  

//      cout<<""<< l <<endl;

      switch (l) // j will contain the local index for the current face
       {
        case 210: case 21: case 102:
        case 120: case 12: case 201:
          j = 0;
          break;  
        case 310: case 31: case 103:
        case 130: case 13: case 301:
          j = 1;
          break;  
        case 321: case 132: case 213:
        case 231: case 123: case 312:
          j = 2;
          break;  
        case 230: case 23: case 302:
        case 320: case 32: case 203:
          j = 3;
          break; 

       default:
        Error("Unable to set the face !!!!!!!!!!!!" << endl);
        exit(0);
      }
      CellTree[Neib[0]]->SetJoint(j, Joint);

    if (Neib[1] != -1) // second element containing the current face
     {
          // find the local index for the point 'a' on the cell
    for (j=0;j<4;j++)
      if (Tetrahedrals[4*Neib[1]+j] == a) break;

    // find the local index for the point 'b' on the cell
    for (k=0;k<4;k++)
      if (Tetrahedrals[4*Neib[1]+k] == b) break;

       // find the local index for the point 'c' on the cell
    for (l=0;l<4;l++)
      if (Tetrahedrals[4*Neib[1]+l] == c) break;   

     l = l*100 + k*10 + j;  

//      cout<<""<< l <<endl;

     switch (l) // j will contain the local index for the current face
      {
        case 210: case 21: case 102:
        case 120: case 12: case 201:
          j = 0;
          break;  
        case 310: case 31: case 103:
        case 130: case 13: case 301:
          j = 1;
          break;  
        case 321: case 132: case 213:
        case 231: case 123: case 312:
          j = 2;
          break;  
        case 230: case 23: case 302:
        case 320: case 32: case 203:
          j = 3;
          break; 

      default:
       Error("Unable to set the face !!!!!!!!!!!!" << endl);
       exit(0);
      }
      CellTree[Neib[1]]->SetJoint(j, Joint);
     } // if (Neib[1] != -1)

  if (Joint->GetType() == InterfaceJoint3D ||
      Joint->GetType() == IsoInterfaceJoint3D)
      {
        ((TInterfaceJoint3D*)Joint)->SetMapType();
        ((TInterfaceJoint3D*)(Joint))->CheckOrientation();
      }
      else 
       if (Joint->GetType() == JointEqN)
           Joint->SetMapType();
   } // if(cell->GetJoint(jj) == NULL)
  } //    for (i=0;i<N_Cells;i++)

  
#ifdef _MPI    
    if(rank==0) 
#endif      
   cout<<"numberoftrifaces after "<<N_G<<endl;
  
#ifdef _MPI    
  if(rank==0)    
#endif
  {
//    In.deinitialize();
//    Out.deinitialize();
  }
  

 
 
#ifdef _MPI    
  if(rank!=0) 
   {
    delete [] Tetrahedrals ;  
    delete [] Coordinates;      
//     delete [] Facelist;
//     delete [] Facemarkerlist;
   }
#endif  
  
  delete [] PointNeighb;

// #ifdef _MPI    
//      if(rank==3)
//     printf(" TetrameshGen  %d \n",   N_G );  
//     MPI_Finalize(); 
// #endif    
//   exit(0);
//      

}


void TetraGenWithInputCells(TDomain *&Domain)
{
 //======================================================================
 // Tetgen for grid generation begin
 //======================================================================
  int i, j, k, l, N_Coord, *N_FVerts, N_Faces, *Facets;
  int N, N_RootCells, N_Cells, CurrVertex, N_Vertices, ID, N_G, RefLevel=0;
  int CurrNeib, len1, len2, len3, maxEpV = 0, a, b, c, Neib[2], Neighb_tmp, CurrComp;
  int *Tetrahedrals, *PointNeighb, dimension;
  int jj, N_Points, MaxLen, *Tetrahedrals_loc;
  const int *EdgeVertex, *TmpFV, *TmpLen;
  double *TetRegionlist;

  double *Coordinates, N_x, N_y, N_z;   
  double *Vertices;
  double Xmin = 1e10, Xmax = -1e10, Ymin = 1e10, Ymax = -1e10;  
  double Zmin = 1e10, Zmax = -1e10, T[4]={0,0,0,0}, S[4]={0,0,0,0}; 
  double StartX, StartY, StartZ, BoundX, BoundY, BoundZ; 
  
  tetgenio In, Out, InAddPts;
   
  TBaseCell **CellTree,  **SurfCellTree;
  TGridCell **DelCell;
  TVertex **VertexDel, **NewVertices, **NewSurfVertices; 
  TBoundPart *BoundPart;
  TBdPlane **UpdateFaceParams;
  TJoint *Joint;
  TBoundComp3D *bdcomp; 
  TCollection *coll, *SurfColl;
  TBaseCell *cell;
  TBdSphere *UpdateParam;  
  TShapeDesc *ShapeDesc;
   
  char *SMESH, line[100];
  SMESH = TDatabase::ParamDB->SMESHFILE;
//   NODE  = TDatabase::ParamDB->GRAPEBASENAME;
  
#ifdef _MPI  
  int rank, size;  
  MPI_Comm Comm = TDatabase::ParamDB->Comm;
    
  MPI_Comm_rank(Comm, &rank);
  MPI_Comm_size(Comm, &size);
 
#endif    

#ifdef _MPI    
  if(rank==0)
#endif  
   {
   
//     In.initialize();
//     Out.initialize();
   
    std::ostringstream opts;
    opts << " "; 
  
    opts.seekp(std::ios::beg);
//     opts<<'p'; // Tetrahedralize the PLC. Switches are chosen to read a PLC (p)
    opts<<'r'; // -r  Reconstructs a previously generated mesh     
//     opts<<'q'<<1.2; // quality mesh generation(q) with a specified quality bound
//     opts<<'a'<<100; // maximum volume constraint
//     opts<<'a'; // maximum volume constraint
//     opts<<'i'; // Inserts a list of additional points into mesh.
    opts<<'z'; // numbers all output items starting from zero
//     opts<<'d'; // Detect intersections of PLC facets.
    opts<<'f'; // Outputs all  faces (including non-boundary) 
    opts<<'e'; // Outputs a list of edges of the triangulation
//     opts<<'I'; // Suppresses mesh iteration numbers.
    opts<<'C'; // Checks the consistency of the final mesh.
//     opts<<'Q'; // Quiet: No terminal output except errors.
//     opts<<'g'; // Outputs mesh to .mesh file for viewing by Medit
    opts<<'Y'; // Suppresses boundary facets/segments splitting
//     opts<<'A'; // output region attribute list
//     opts<<'V';  //verbose mode
//     opts<<'V';
//     opts<<'V';
//     opts<<'V';
//     opts<<'V';    
    opts<<ends;
    
     In.initialize();
     Out.initialize();   
//      InAddPts.initialize();     
//   cout << " SMESHFILE is " << SMESH << endl;
  ReadMeditMeshWithCells(SMESH, In);
//   ReadAdditionalNModes(NODE, InAddPts);   
//      In.load_medit(SMESH, 1);
     

//      InAddPts.load_node(NODE);
//      cout<< " InAddPts.numberofpoints " << InAddPts.numberofpoints <<endl; 
//      exit(0);
//     for(i=0;i<In.numberofpoints;i++)
//       OutPut(i<<" (x, y, z) =  "<<
//        In.pointlist[3*i]<<' '<<In.pointlist[3*i+1]<<' '<<In.pointlist[3*i+2]<<endl);
//      exit(0);
 // Calling  tetrahedralize function of 3dtetgen mesh generator 


        tetrahedralize((char*)opts.str().c_str(), &In, &Out);
//     tetrahedralize((char*)opts.str().c_str(), &In, &Out, &InAddPts);
//     exit(0);
 //   output: coordinates of all vertices
//  for(i=0;i<Out.numberofpoints;i++)
//   OutPut(" (x, y, z) =  "<< Out.pointlist[3*i]<<' '<<Out.pointlist[3*i+1]<<' '<<Out.pointlist[3*i+2]<<endl);    
    
   } // if(rank==0)
     
  DeleteDomain(Domain);

   
#ifdef _MPI     
     if(rank==0) 
#endif  
     {
       N_RootCells = Out.numberoftetrahedra;
       N_G = Out.numberofpoints;   

       // allocate auxillary fields
       Coordinates = Out.pointlist;  
       Tetrahedrals = Out.tetrahedronlist;   
       
       TetRegionlist = Out.tetrahedronattributelist;
     }

#ifdef _MPI
   MPI_Bcast(&N_RootCells, 1, MPI_INT, 0, Comm);
   MPI_Bcast(&N_G, 1, MPI_INT, 0, Comm);  
   
   
  if(rank!=0)
   {
    Coordinates = new double [3*N_G];
    Tetrahedrals = new int[4*N_RootCells];
    TetRegionlist = new double [N_RootCells];    
   }
   
   MPI_Bcast(Coordinates, 3*N_G, MPI_DOUBLE, 0, Comm);  
   MPI_Bcast(Tetrahedrals, 4*N_RootCells, MPI_INT, 0, Comm);
   MPI_Bcast(TetRegionlist, N_RootCells, MPI_DOUBLE, 0, Comm);    
#endif      
   
  // generate new vertices
   NewVertices = new TVertex*[N_G];

   for (i=0;i<N_G;i++)
    {
      NewVertices[i] = new TVertex(Coordinates[3*i], Coordinates[3*i+1], Coordinates[3*i+2]);

      // set bounding box
      if (Coordinates[3*i] > Xmax) Xmax = Coordinates[3*i];
      if (Coordinates[3*i] < Xmin) Xmin = Coordinates[3*i];
      if (Coordinates[3*i+1] > Ymax) Ymax = Coordinates[3*i+1];
      if (Coordinates[3*i+1] < Ymin) Ymin = Coordinates[3*i+1];
      if (Coordinates[3*i+2] > Zmax) Zmax = Coordinates[3*i+2];
      if (Coordinates[3*i+2] < Zmin) Zmin = Coordinates[3*i+2];
    }

   // set bounding box
    StartX = Xmin;
    StartY = Ymin;
    StartZ = Zmin;
    BoundX = Xmax - Xmin;
    BoundY = Ymax - Ymin;
    BoundZ = Zmax - Zmin;


   Domain->SetBoundBox(StartX, StartY, StartZ, BoundX, BoundY, BoundZ);   
//        cout<<Xmin <<"  "<<Ymin <<"  "<<Zmin<<endl;
//        cout<<Xmax <<"  "<<Ymax <<"  "<<Zmax<<endl;

   CellTree = new TBaseCell*[N_RootCells];
 //  output of each tetraheron vertex indices (four vertices for each)
//   for (i=0;i<N_RootCells;i++)
//        cout<< Tetrahedrals[4*i]<<"  "<<Tetrahedrals[4*i + 1]<<"  "
//          <<Tetrahedrals[4*i + 2]<<"  "<<Tetrahedrals[4*i + 3]<<endl;


   for (i=0;i<N_RootCells;i++)
    {
     CellTree[i] = new TMacroCell(TDatabase::RefDescDB[Tetrahedron],
                                    RefLevel);

     CellTree[i]->SetVertex(0, NewVertices[Tetrahedrals[4*i    ]]);
     CellTree[i]->SetVertex(1, NewVertices[Tetrahedrals[4*i + 1]]);
     CellTree[i]->SetVertex(2, NewVertices[Tetrahedrals[4*i + 2]]);
     CellTree[i]->SetVertex(3, NewVertices[Tetrahedrals[4*i + 3]]);

     CellTree[i]->SetClipBoard(i);
      ((TMacroCell *) CellTree[i])->SetSubGridID(0);
      
//        cout << i<<  "  TetRegionlist :  " <<  TetRegionlist[i] << endl;     
       CellTree[i]->SetRegionID((int)TetRegionlist[i]);
    }
//      exit(0);
   Domain->SetTreeInfo(CellTree, N_RootCells);


   // initialize iterators
   TDatabase::IteratorDB[It_EQ]->SetParam(Domain);
   TDatabase::IteratorDB[It_LE]->SetParam(Domain);
   TDatabase::IteratorDB[It_Finest]->SetParam(Domain);
   TDatabase::IteratorDB[It_Between]->SetParam(Domain);
   TDatabase::IteratorDB[It_OCAF]->SetParam(Domain);
   
  
   // search neighbours
   PointNeighb = new int[N_G];
#ifdef _MPI    
     if(rank==0)   
#endif       
   cout<<"numberofpoints "<<N_G<<endl;
   memset(PointNeighb, 0, N_G *SizeOfInt);

     for (i=0;i<4*N_RootCells;i++)
     PointNeighb[Tetrahedrals[i]]++;

   for (i=0;i<N_G;i++)
     if (PointNeighb[i] > maxEpV) maxEpV = PointNeighb[i];
   delete [] PointNeighb;

#ifdef _MPI    
     if(rank==0)   
#endif   
   cout<<"maxEpV "<< maxEpV<<endl;

   PointNeighb = new int[++maxEpV * N_G];

   memset(PointNeighb, 0, maxEpV*N_G*SizeOfInt);

    // every vertex contains "maxEpV" columns
    // for every vertex at first colomn contains the number of cells containing this vertex
    // at further columns we set the index of corresponding cells containing this vertex
    // cout<<"maxEpV*N_G "<<maxEpV*N_G<<endl;

   for(i=0;i<4*N_RootCells;i++)
    {
     j = Tetrahedrals[i]*maxEpV;
     PointNeighb[j]++;
     //cout<<"j + PointNeighb[j] " << j <<endl;
     PointNeighb[j + PointNeighb[j]] = i / 4;
    }


   // generate new faces    
   coll=Domain->GetCollection(It_Finest, 0);
   N_Cells = coll->GetN_Cells();       
   N_G = 0;
   for (i=0;i<N_Cells;i++)
   {
    cell = coll->GetCell(i);
    ShapeDesc= cell->GetShapeDesc();   
    ShapeDesc->GetFaceVertex(TmpFV, TmpLen, MaxLen);      
    N_Faces = cell->GetN_Faces();
    Tetrahedrals_loc = Tetrahedrals+4*i;
           
    for(jj=0;jj<N_Faces;jj++)
     if(cell->GetJoint(jj) == NULL)
      {
       N_G++;    
       N_Points = TmpLen[jj];
   
       if(N_Points!=3)
        {     
         cerr << "Only Tria faces are allowed!!! N_FVert: " << N_Points << endl;
         exit(-1);     
        }     
     
     
       //printf(" TetrameshGen  Null joint  \n" );             
       a = Tetrahedrals_loc[TmpFV[jj*MaxLen]];
       b = Tetrahedrals_loc[TmpFV[jj*MaxLen + 1]];
       c = Tetrahedrals_loc[TmpFV[jj*MaxLen + 2]];

//      cout<<"  "<< a<<"  "<< b<<"  "<< c<<endl;

       Neib[0] = -1;
       Neib[1] = -1;
       CurrNeib = 0;

       len1 = PointNeighb[a*maxEpV];
       len2 = PointNeighb[b*maxEpV];
       len3 = PointNeighb[c*maxEpV];

       // find the index of the cells containing current face with point indices a,b,c
       for (j=1;j<=len1;j++)
       {
        Neighb_tmp = PointNeighb[a*maxEpV + j];
        for (k=1;k<=len2;k++)
         {
          if (Neighb_tmp == PointNeighb[b*maxEpV + k])
           {
            for (l=1;l<=len3;l++)
             if (Neighb_tmp == PointNeighb[c*maxEpV + l])
             {
              Neib[CurrNeib++] = Neighb_tmp;
              break;
             }
            }
           }
         if (CurrNeib == 2) break;
        }// for (j=1;j<=len1;j++)
        
      //   cout<<"CurrNeib " << CurrNeib <<endl;
      // cout<<"Facemarkerlist[i] : "<<Facemarkerlist[i]<<endl;
      if (CurrNeib == 1) // 0 for inner edges and Boundcomp+1 for Boundedge respect
       {
        CurrComp = 0; // not yet implemented fully
       
       
       
        bdcomp = Domain->GetBdPart(0)->GetBdComp(CurrComp);
       
 
        if(bdcomp->GetTSofXYZ(NewVertices[a]->GetX(), NewVertices[a]->GetY(),
                             NewVertices[a]->GetY(), T[1], S[1]) ||
          bdcomp->GetTSofXYZ(NewVertices[b]->GetX(), NewVertices[b]->GetY(),
                             NewVertices[b]->GetY(), T[2], S[2]) ||
          bdcomp->GetTSofXYZ(NewVertices[c]->GetX(), NewVertices[c]->GetY(),
                             NewVertices[c]->GetY(), T[3], S[3])    )
         {
          cerr<<"Error: could not set parameter values"<<endl;
          OutPut(NewVertices[a]<<endl);
          OutPut(NewVertices[b]<<endl);
          OutPut(NewVertices[c]<<endl);
          exit(0);
         }
 
       if (CurrNeib == 2)
        {
         if(bdcomp->IsFreeBoundary())
           Joint = new TIsoInterfaceJoint3D(bdcomp, T, S,
                       CellTree[Neib[0]], CellTree[Neib[1]] );
          else
           Joint = new TInterfaceJoint3D(bdcomp, T, S, 
                       CellTree[Neib[0]], CellTree[Neib[1]] );
        }
       else
        {
         if(bdcomp->IsFreeBoundary())
           Joint = new TIsoBoundFace(bdcomp, T, S);
         else
           Joint = new TBoundFace(bdcomp, T, S);
        }
      }
     else
      {
       // cout<<"Inner face"<<endl;
       if (CurrNeib != 2)
       cerr << "Error !!!!!!!! not enough neighbours!" << endl;

       Joint = new TJointEqN(CellTree[Neib[0]], CellTree[Neib[1]]);
      }

     // First element containing the current face
     // find the local index for the point 'a' on the cell
     for (j=0;j<4;j++)
      if (Tetrahedrals[4*Neib[0]+j] == a) break;

     // find the local index for the point 'b' on the cell
     for (k=0;k<4;k++)
      if (Tetrahedrals[4*Neib[0]+k] == b) break;

      // find the local index for the point 'c' on the cell
     for (l=0;l<4;l++)
      if (Tetrahedrals[4*Neib[0]+l] == c) break;   

      l = l*100 + k*10 + j;  

//      cout<<""<< l <<endl;

      switch (l) // j will contain the local index for the current face
       {
        case 210: case 21: case 102:
        case 120: case 12: case 201:
          j = 0;
          break;  
        case 310: case 31: case 103:
        case 130: case 13: case 301:
          j = 1;
          break;  
        case 321: case 132: case 213:
        case 231: case 123: case 312:
          j = 2;
          break;  
        case 230: case 23: case 302:
        case 320: case 32: case 203:
          j = 3;
          break; 

       default:
        Error("Unable to set the face !!!!!!!!!!!!" << endl);
        exit(0);
      }
      CellTree[Neib[0]]->SetJoint(j, Joint);

    if (Neib[1] != -1) // second element containing the current face
     {
          // find the local index for the point 'a' on the cell
    for (j=0;j<4;j++)
      if (Tetrahedrals[4*Neib[1]+j] == a) break;

    // find the local index for the point 'b' on the cell
    for (k=0;k<4;k++)
      if (Tetrahedrals[4*Neib[1]+k] == b) break;

       // find the local index for the point 'c' on the cell
    for (l=0;l<4;l++)
      if (Tetrahedrals[4*Neib[1]+l] == c) break;   

     l = l*100 + k*10 + j;  

//      cout<<""<< l <<endl;

     switch (l) // j will contain the local index for the current face
      {
        case 210: case 21: case 102:
        case 120: case 12: case 201:
          j = 0;
          break;  
        case 310: case 31: case 103:
        case 130: case 13: case 301:
          j = 1;
          break;  
        case 321: case 132: case 213:
        case 231: case 123: case 312:
          j = 2;
          break;  
        case 230: case 23: case 302:
        case 320: case 32: case 203:
          j = 3;
          break; 

      default:
       Error("Unable to set the face !!!!!!!!!!!!" << endl);
       exit(0);
      }
      CellTree[Neib[1]]->SetJoint(j, Joint);
     } // if (Neib[1] != -1)

  if (Joint->GetType() == InterfaceJoint3D ||
      Joint->GetType() == IsoInterfaceJoint3D)
      {
        ((TInterfaceJoint3D*)Joint)->SetMapType();
        ((TInterfaceJoint3D*)(Joint))->CheckOrientation();
      }
      else 
       if (Joint->GetType() == JointEqN)
           Joint->SetMapType();
   } // if(cell->GetJoint(jj) == NULL)
  } //    for (i=0;i<N_Cells;i++)

  
#ifdef _MPI    
    if(rank==0) 
#endif      
   cout<<"numberoftrifaces after "<<N_G<<endl;
  
#ifdef _MPI    
  if(rank==0)    
#endif
  {
//    In.deinitialize();
//    Out.deinitialize();
  }
  

 
 
#ifdef _MPI    
  if(rank!=0) 
   {
    delete [] Tetrahedrals ;  
    delete [] Coordinates;      
    delete [] TetRegionlist;
//     delete [] Facemarkerlist;
   }
#endif  
  
  delete [] PointNeighb;

// #ifdef _MPI    
//      if(rank==3)
//     printf(" TetrameshGen  %d \n",   N_G );  
//     MPI_Finalize(); 
// #endif    
//   exit(0);
//      

}