//---------------------------------------------------------------------------


#pragma hdrstop

#include "MC.h"
#include <math.h>
//---------------------------------------------------------------------------
void findFirstCube(TProcess *p);

TCubes *findCube(TProcess *p, int i, int j, int k);
TCubeEdge *findEdge(TProcess *p, TCubeVertex *v1, TCubeVertex *v2);

void subdivideCube(TProcess *p, TCubes *Cube);
TCubes *findCubeAt(TProcess *p, int i, int j, int k, int subdivision);
TCubes *createCubeAt(TProcess *p, int i, int j, int k, int subdivision);
TCubes *findNeighbour(TProcess *p, TCube *Cube, int dir);

//---------------------------------------------------------------------------



void CreateCubes(TProcess *p){
  TCubes *LastCube; // last cube of the list
  int Directions[6];

  p->Stats.Cubes.Start = clock();

  findFirstCube(p);

  p->ActiveCube = p->ActiveCubes;
// loop at active cubes
  while(p->ActiveCube != NULL){
    if(p->quit) break;
    TCubes *AC = p->ActiveCube;

    if(AC->Cube.SubdivisionLevel == 4){
      int i = 0;
    }

    if(AC->Cube.Evaluate){
//....determine configuration
      AC->Cube.Configuration = (AC->Cube.Vertices[0]->Value < 0.0) +
															((AC->Cube.Vertices[1]->Value < 0.0) << 1 ) +
															((AC->Cube.Vertices[2]->Value < 0.0) << 2 ) +
															((AC->Cube.Vertices[3]->Value < 0.0) << 3 ) +
															((AC->Cube.Vertices[4]->Value < 0.0) << 4 ) +
															((AC->Cube.Vertices[5]->Value < 0.0) << 5 ) +
															((AC->Cube.Vertices[6]->Value < 0.0) << 6 ) +
															((AC->Cube.Vertices[7]->Value < 0.0) << 7 );



      MCALLBACK(p->Callback, p, S_MC, A_MC_EVALUATE)

      if(AC->Cube.SubdivisionLevel == 0){
//......find neighbours for root cubes
        for(int i = 0; i < 6; i++){ // loop trough all traversal faces
          int dir = C_CONF_DIR[ CCALL[AC->Cube.Configuration][C_CON] ][i];
          TCubes *NewCube;
          
          if(dir == NO_DIR) //..........no more traversal faces, so have a
            break;

          dir = C_ROT_DIR[ CCALL[AC->Cube.Configuration][C_ROT] ][dir];
//........find/create cube in direction dir
          NewCube = findCube(p, AC->Cube.Index.I + (dir == 0 ? -1 : dir == 1 ? 1 : 0) * p->Config.RC_RATIO,
                                AC->Cube.Index.J + (dir == 2 ? 1 : dir == 3 ? -1 : 0) * p->Config.RC_RATIO,
                                AC->Cube.Index.K + (dir == 4 ? 1 : dir == 5 ? -1 : 0) * p->Config.RC_RATIO);

//........its new cube
          for(int vdir = 0; vdir < 8; vdir++){
            if(NewCube->Cube.Vertices[vdir] == NULL){
							NewCube->Cube.Vertices[vdir] = findCubeVertex(p,
                    NewCube->Cube.Index.I + IS_RIGHT(vdir) * p->Config.RC_RATIO,
                    NewCube->Cube.Index.J + IS_TOP(vdir) * p->Config.RC_RATIO,
                    NewCube->Cube.Index.K + IS_FRONT(vdir) * p->Config.RC_RATIO);
						}
          }
//........set neighbourhood
//	  			AC->Cube.Neighbour[dir] = &(NewCube->Cube);
//  				NewCube->Cube.Neighbour[BACK(dir)] = &(AC->Cube);
        }
      } // endif(AC->Cube.SubdivisionLevel == 0)
//......find edges
        for(int i = 0; i < 12; i++){
          int edge = C_TRAVERSAL_EDGES[CCALL[AC->Cube.Configuration][C_CON]][i]; //
//........no more traversal edges
          if(edge == NO_EDGE) break;
//........set left edge (left in CCW direction)
          int edge_left = CROT[ CCALL[AC->Cube.Configuration][C_ROT] ]
                              [ CCONF[CCALL[AC->Cube.Configuration][C_CON] ][edge][CLEFT] ];
          int edge_right = CROT[ CCALL[AC->Cube.Configuration][C_ROT] ]
                              [ CCONF[CCALL[AC->Cube.Configuration][C_CON] ][edge][CRIGHT] ];
//........rotate edge
          edge = CROT[CCALL[AC->Cube.Configuration][C_ROT]][edge];
//........find/create edge
          AC->Cube.Edges[edge] = findEdge(p,
                                AC->Cube.Vertices[C_EDGE_VERTICES[edge][0]],
                                AC->Cube.Vertices[C_EDGE_VERTICES[edge][1]]);

//........set adjecant cube for edge
          AC->Cube.Edges[edge]->Cubes[EDGE_NEIGHBOUR[edge]] = AC;
//........set neighbours for edges
//........all neighbours have same level
          AC->Cube.Edges[edge]->Neighbours[EDGE_NEIGHBOUR[edge]]
            = AC->Cube.Edges[edge_left];

          if(AC->Cube.Edges[edge_right] != NULL){
            AC->Cube.Edges[edge_right]->Neighbours[EDGE_NEIGHBOUR[edge_right]]
              = AC->Cube.Edges[edge];
          }

          /**      AC->Cube.Edges[edge]->Neighbours[ EDGE_NEIGHBOUR[ EDGE_BACK_NEIGHBOUR[edge][edge_right] ]]
            = AC->Cube.Edges[edge_right]; /**/

/**      AC->Cube.Edges[edge_left]->Neighbours[BACK_IN_PATCH(C..N1[edge_left])]
            = AC->Cube.Edges[edge]; /**/
        }


      if(p->Config.SUBDIVIDE){
//......check subdivision
        for(int i = 0; i < 12; i++){
          int edge = C_TRAVERSAL_EDGES[CCALL[AC->Cube.Configuration][C_CON]][i];
//........no more traversal edges
          if(edge == NO_EDGE) break;
//........rotate edge
          edge = CROT[CCALL[AC->Cube.Configuration][C_ROT]][edge];

//........check if all neighbours at same sbdLvl are set!!
/*          if((AC->Cube.Edges[edge]->Neighbours[EAST] == NULL)
            || (AC->Cube.Edges[edge]->Neighbours[NORTH] == NULL)
            || (AC->Cube.Edges[edge]->Neighbours[WEST] == NULL)
            || (AC->Cube.Edges[edge]->Neighbours[SOUTH] == NULL)
            || AC->Cube.Edges[edge]->Subdivide)
            continue;       */
          if((AC->Cube.Edges[edge]->Cubes[EAST] == NULL)
            || (AC->Cube.Edges[edge]->Cubes[NORTH] == NULL)
            || (AC->Cube.Edges[edge]->Cubes[WEST] == NULL)
            || (AC->Cube.Edges[edge]->Cubes[SOUTH] == NULL)
            || AC->Cube.Edges[edge]->Subdivide)
            continue;
//........if there was a cube vertex f(v)==0
//..........special method !!!!!!!!!!!!!!!!!!

//........check subdivision
          if( (fabs((AC->Cube.Edges[edge]->Neighbours[EAST]->Intersection - AC->Cube.Edges[edge]->Intersection)
              * (AC->Cube.Edges[edge]->Intersection - AC->Cube.Edges[edge]->Neighbours[WEST]->Intersection))
              / ( sqrt((AC->Cube.Edges[edge]->Neighbours[EAST]->Intersection - AC->Cube.Edges[edge]->Intersection)
                  * (AC->Cube.Edges[edge]->Neighbours[EAST]->Intersection - AC->Cube.Edges[edge]->Intersection) )
                 * sqrt((AC->Cube.Edges[edge]->Intersection - AC->Cube.Edges[edge]->Neighbours[WEST]->Intersection)
                  * (AC->Cube.Edges[edge]->Intersection - AC->Cube.Edges[edge]->Neighbours[WEST]->Intersection))
              )
            < pow(p->Config.ALPHA_SUBDIVISON, 1.0 / (AC->Cube.SubdivisionLevel + 1)) )
//              < p->Config.ALPHA_SUBDIVISON )
            || (fabs((AC->Cube.Edges[edge]->Neighbours[NORTH]->Intersection - AC->Cube.Edges[edge]->Intersection)
              * (AC->Cube.Edges[edge]->Intersection - AC->Cube.Edges[edge]->Neighbours[SOUTH]->Intersection))
              / ( sqrt((AC->Cube.Edges[edge]->Neighbours[NORTH]->Intersection - AC->Cube.Edges[edge]->Intersection)
                  * (AC->Cube.Edges[edge]->Neighbours[NORTH]->Intersection - AC->Cube.Edges[edge]->Intersection) )
                * sqrt((AC->Cube.Edges[edge]->Intersection - AC->Cube.Edges[edge]->Neighbours[SOUTH]->Intersection)
                  * (AC->Cube.Edges[edge]->Intersection - AC->Cube.Edges[edge]->Neighbours[SOUTH]->Intersection)) )
            < pow(p->Config.ALPHA_SUBDIVISON, 1.0 / (AC->Cube.SubdivisionLevel + 1)) )
//              < p->Config.ALPHA_SUBDIVISON)
              ){

            AC->Cube.Edges[edge]->Subdivide = 1;

//..........set all cubes sharing this edge to subdide
            for(int i = 0; i < 4; i++){
//............if edge wasn't eyt set (subd/evaulate) add to list
              if((AC->Cube.Edges[edge]->Cubes[i]->Cube.Subdivide == 0)
                  && (AC->Cube.Edges[edge]->Cubes[i]->Cube.Evaluate == 0)){
                AC->Cube.Edges[edge]->Cubes[i]->Next = NULL;
                p->ActiveCubes->Next = AC->Cube.Edges[edge]->Cubes[i];
                p->ActiveCubes = AC->Cube.Edges[edge]->Cubes[i];
              }
              AC->Cube.Edges[edge]->Cubes[i]->Cube.Subdivide = 1;
            }
          }
        } // end of subdivision check
      } // endif(p->Config.SUBDIVIDE)
      AC->Cube.Evaluate = 0;
    } // endif(AC->Cube.Evaluate)



    if(AC->Cube.Subdivide){
      if((p->Config.RC_RATIO >> AC->Cube.SubdivisionLevel) > 1){
//      if(AC->Cube.SubdivisionLevel < 4){
//......subdivision needed
        MCALLBACK(p->Callback, p, S_MC, A_MC_SUBDIVIDE)
//......subdivide (create only if not exists)
        subdivideCube(p, AC);

        for(int face = 0; face < 6; face++){
          Directions[face] = -1;
        }

//......check if new intersectoins were found (NOT IMPLEMENTED YET)
        for(int edge = 0; edge < 12; edge++){
//........check non-existing edges
          if(AC->Cube.Edges[edge] != NULL) continue; // with next edge

          if((AC->Cube.Children[EDGE_CHILDREN[edge][0]]->Cube.Vertices[C_EDGE_VERTICES[edge][0]]->Value < 0.0)
              == (AC->Cube.Children[EDGE_CHILDREN[edge][0]]->Cube.Vertices[C_EDGE_VERTICES[edge][1]]->Value < 0.0)) continue;
//........create cubes up to this subdivision level (if needed)
//........set neighbours for subdivision
          Directions[C_EDGE_DIRS[edge][0]] = 1;
          Directions[C_EDGE_DIRS[edge][1]] = 1;

          for(int dir = 0; dir < 3; dir++){

/**/        TCubes * FoundCube = createCubeAt(p,
                       AC->Cube.Index.I + C_EDGE_DELTA[edge][dir][0] * (p->Config.RC_RATIO >> AC->Cube.SubdivisionLevel),
                       AC->Cube.Index.J + C_EDGE_DELTA[edge][dir][1] * (p->Config.RC_RATIO >> AC->Cube.SubdivisionLevel),
                       AC->Cube.Index.K + C_EDGE_DELTA[edge][dir][2] * (p->Config.RC_RATIO >> AC->Cube.SubdivisionLevel),
                         AC->Cube.SubdivisionLevel);

            if(FoundCube->Cube.Subdivide != 1){
              FoundCube->Cube.Subdivide = 1;
              if(FoundCube->Cube.Evaluate != 1){
                FoundCube->Next = NULL;
                p->ActiveCubes->Next = FoundCube;
                p->ActiveCubes = FoundCube;
              }
            } /**/
          }

        }

/**/        for(int face = 0; face < 6; face++){
//........check faces
          int dir = C_CONF_IS_DIR[ CCALL[AC->Cube.Configuration][C_CON] ][face];
          TCubes *NewCube;
          if(dir != NO_DIR) //..........no more traversal faces, so continue with another
            continue;

          dir = C_ROT_DIR[ CCALL[AC->Cube.Configuration][C_ROT] ][face];
          if(Directions[dir] != NO_DIR) continue;
//........check if subdivision is needed
          if( (AC->Cube.Vertices[DIR_VERTICES[dir][0]]->Value < 0.0)
            == (AC->Cube.Children[DIR_VERTICES[dir][0]]->Cube.Vertices[DIR_VERTICES[dir][3]]->Value < 0.0) ) continue;
//........create cubes up to this subdivision level (if needed)
// face 0-left
//        1-right
//        2-up
//        3-down
//        4-front
//        5-back
          TCubes * FoundCube = createCubeAt(p,
                       AC->Cube.Index.I + (dir == 0 ? -1 : (dir == 1 ? 1 : 0)) * (p->Config.RC_RATIO >> AC->Cube.SubdivisionLevel),
                       AC->Cube.Index.J + (dir == 2 ?  1 : (dir == 3 ?-1 : 0)) * (p->Config.RC_RATIO >> AC->Cube.SubdivisionLevel),
                       AC->Cube.Index.K + (dir == 4 ?  1 : (dir == 5 ?-1 : 0)) * (p->Config.RC_RATIO >> AC->Cube.SubdivisionLevel),
                         AC->Cube.SubdivisionLevel);

          if(FoundCube->Cube.Subdivide != 1){
            FoundCube->Cube.Subdivide = 1;
            if(FoundCube->Cube.Evaluate != 1){
              p->ActiveCubes->Next = FoundCube;
              p->ActiveCubes = FoundCube;
            }
          }

        }    /**/
      }
      else {
        AC->Cube.Subdivide = 0;
//......edge detection
      }
    }

    p->ActiveCube = p->ActiveCube->Next;
  }

  p->Stats.Cubes.End = clock();
}

//---------------------------------------------------------------------------

void findFirstCube(TProcess *p){
  TCubes *Cubes = findCube(p,
          (int)(floor(p->Config.STARTING_POSITION.X / p->Config.BIG_STEP)) * p->Config.RC_RATIO,
          (int)(floor(p->Config.STARTING_POSITION.Y / p->Config.BIG_STEP)) * p->Config.RC_RATIO,
          (int)(floor(p->Config.STARTING_POSITION.Z / p->Config.BIG_STEP)) * p->Config.RC_RATIO);


//Set first 8 Cube-Vertices
  for(int dir = 0; dir < 8; dir++){
    Cubes->Cube.Vertices[dir] = findCubeVertex(p,
                Cubes->Cube.Index.I + IS_RIGHT(dir) * p->Config.RC_RATIO,
                Cubes->Cube.Index.J + IS_TOP(dir) * p->Config.RC_RATIO,
                Cubes->Cube.Index.K + IS_FRONT(dir) * p->Config.RC_RATIO);

//    CV->Cubes[BACK(dir)] = &(p.ActiveCubes->Cube);
  }

  MCALLBACK(p->Callback, p, S_MC, A_MC_FIRST_CUBE)
}

//---------------------------------------------------------------------------

TCubeVertex* findCubeVertex(TProcess *p, int i, int j, int k)
{
  long int Index = HASH(i / (p->Config.RC_RATIO / 4),
                        j / (p->Config.RC_RATIO / 4),
                        k / (p->Config.RC_RATIO / 4));
  TCubeVertices *NewCV;
  TCubeVertices *CV = p->FoundVertices[Index];

  for(;CV != NULL;CV = CV->Next)
  {
    if((CV->Vertex.Index.I == i)
        && (CV->Vertex.Index.J == j)
        && (CV->Vertex.Index.K == k)){
      return &(CV->Vertex);
    }
  }
  NewCV = (TCubeVertices *) mycalloc(1, sizeof(TCubeVertices));
  NewCV->Vertex.Index.I = i;
  NewCV->Vertex.Index.J = j;
  NewCV->Vertex.Index.K = k;
  NewCV->Vertex.Value = p->Function(p->Config.STEP * i,
                                    p->Config.STEP * j,
                                    p->Config.STEP * k);

  if(NewCV->Vertex.Value == 0.0){
    NewCV->Vertex.Value = p->Config.MIN_VALUE;
  }

  NewCV->Next = p->FoundVertices[Index];
  p->FoundVertices[Index] = NewCV;

  return &(NewCV->Vertex);

}


//---------------------------------------------------------------------------

TCubes * findCube(TProcess *p, int i, int j, int k){
  int Index = HASH(i / p->Config.RC_RATIO,
                   j / p->Config.RC_RATIO,
                   k / p->Config.RC_RATIO);
  TCubesHash *NewCube;
  TCubesHash *Cubes = p->FoundCubes[Index];

  for(;Cubes != NULL;Cubes = Cubes->Next)
  {
    if((Cubes->Cube.Cube.Index.I == i)
        && (Cubes->Cube.Cube.Index.J == j)
        && (Cubes->Cube.Cube.Index.K == k)){
      return &(Cubes->Cube);
    }
  }
  NewCube = (TCubesHash *) mycalloc(1, sizeof(TCubesHash));
  NewCube->Cube.Cube.Index.I = i;
  NewCube->Cube.Cube.Index.J = j;
  NewCube->Cube.Cube.Index.K = k;
  NewCube->Cube.Cube.Evaluate = 1;

  NewCube->Next = p->FoundCubes[Index];
  p->FoundCubes[Index] = NewCube;

//change to stack  !!!!!
  if(p->ActiveCubes != NULL)
    p->ActiveCubes->Next = &(NewCube->Cube);
  p->ActiveCubes = &(NewCube->Cube);

  p->Stats.Cubes.Count++;
  return &(NewCube->Cube);

}

//---------------------------------------------------------------------------

TCubeEdge * findEdge(TProcess *p, TCubeVertex *v1, TCubeVertex *v2){
  int Index = HASH((v1->Index.I + v2->Index.I) / (p->Config.RC_RATIO / 4),
                   (v1->Index.J + v2->Index.J) / (p->Config.RC_RATIO / 4),
                   (v1->Index.K + v2->Index.K) / (p->Config.RC_RATIO / 4));
  double t;


  TCubeEdges *NewEdge;
  TCubeEdges *Edges = p->FoundEdges[Index];

  for(;Edges != NULL;Edges = Edges->Next)
  {
    if(((Edges->Edge.Vertices[0] == v1) && (Edges->Edge.Vertices[1] == v2))
      || ((Edges->Edge.Vertices[1] == v1) && (Edges->Edge.Vertices[0] == v2))){
      return &(Edges->Edge);
    }
  }
  NewEdge = (TCubeEdges *) mycalloc(1, sizeof(TCubeEdges));
  NewEdge->Edge.Vertices[0] = v1;
  NewEdge->Edge.Vertices[1] = v2;

//find intersection
  if((v1->Value < 0.0) != (v2->Value < 0.0)){
    t = v1->Value / (v1->Value - v2->Value) ;
    NewEdge->Edge.Intersection.X = ((t * (double) v2->Index.I )
                                    + ((1 - t) * (double) v1->Index.I)) * p->Config.STEP;

    NewEdge->Edge.Intersection.Y = ((t * (double) v2->Index.J )
                                    + ((1 - t) * (double) v1->Index.J)) * p->Config.STEP;

    NewEdge->Edge.Intersection.Z = ((t * (double) v2->Index.K )
                                    + ((1 - t) * (double) v1->Index.K)) * p->Config.STEP;

//    ((v1->Value * p->Config.STEP * (double) v2->Index.K )
//                                  - (v2->Value * p->Config.STEP * (double) v1->Index.K))
//                                        / (v1->Value - v2->Value);



  }
/*
  if( (NewEdge->Edge.Intersection.X < 0.1) && //(NewEdge->Edge.Intersection.X > -0.1) &&
    (NewEdge->Edge.Intersection.Y < 0.1) && // (NewEdge->Edge.Intersection.Y > -0.1) &&
    (NewEdge->Edge.Intersection.Z < 0.1) //&& (NewEdge->Edge.Intersection.Z > -0.1)){
    ){
    int i;
    i =0;
  }
*/

  NewEdge->Next = p->FoundEdges[Index];
  p->FoundEdges[Index] = NewEdge;

  return &(NewEdge->Edge);

}


//---------------------------------------------------------------------------
void subdivideCube(TProcess *p, TCubes *Cube){
  for(int i = 0; i < 8; i++){
    if(Cube->Cube.Children[i] == NULL){
      TCubes *NewCube = (TCubes *) mycalloc(1, sizeof(TCubesHash));
      NewCube->Cube.Parent = Cube;
      NewCube->Cube.SubdivisionLevel = Cube->Cube.SubdivisionLevel + 1;
      NewCube->Cube.Index.I = Cube->Cube.Index.I + IS_RIGHT(i) * ( p->Config.RC_RATIO >> NewCube->Cube.SubdivisionLevel);
      NewCube->Cube.Index.J = Cube->Cube.Index.J + IS_TOP(i) * ( p->Config.RC_RATIO >> NewCube->Cube.SubdivisionLevel);
      NewCube->Cube.Index.K = Cube->Cube.Index.K + IS_FRONT(i) * ( p->Config.RC_RATIO >> NewCube->Cube.SubdivisionLevel);

      NewCube->Cube.Evaluate = 1;
      NewCube->Cube.Subdivide = 0;

      p->Stats.Cubes.Count++;

      p->ActiveCubes->Next = NewCube;
      p->ActiveCubes = NewCube;

      for(int vdir = 0; vdir < 8; vdir++){
        if(NewCube->Cube.Vertices[vdir] == NULL){
				 	NewCube->Cube.Vertices[vdir] = findCubeVertex(p,
             NewCube->Cube.Index.I + IS_RIGHT(vdir) * ( p->Config.RC_RATIO >> NewCube->Cube.SubdivisionLevel),
             NewCube->Cube.Index.J + IS_TOP(vdir) * ( p->Config.RC_RATIO >> NewCube->Cube.SubdivisionLevel),
             NewCube->Cube.Index.K + IS_FRONT(vdir) * ( p->Config.RC_RATIO >> NewCube->Cube.SubdivisionLevel));
        }
      }

      Cube->Cube.Children[i] = NewCube;
    }
  }
}


//---------------------------------------------------------------------------
// return NULL if not exists
TCubes *findCubeAt(TProcess *p, int i, int j, int k, int subdivision){

}

//---------------------------------------------------------------------------
TCubes *createCubeAt(TProcess *p, int i, int j, int k, int subdivision){
  int ratio = p->Config.RC_RATIO;
/**/  TCubes *Parent = findCube(p, (i / ratio) * ratio - ((i % ratio) < 0 ? ratio : 0),
                               (j / ratio) * ratio - ((j % ratio) < 0 ? ratio : 0),
                               (k / ratio) * ratio - ((k % ratio) < 0 ? ratio : 0));/**/
/**  TCubes *Parent = findCube(p, (i / ratio) * ratio,
                               (j / ratio) * ratio,
                               (k / ratio) * ratio);    /**/
  for(int vdir = 0; vdir < 8; vdir++){
    if(Parent->Cube.Vertices[vdir] == NULL){
  		Parent->Cube.Vertices[vdir] = findCubeVertex(p,
      Parent->Cube.Index.I + IS_RIGHT(vdir) * p->Config.RC_RATIO,
      Parent->Cube.Index.J + IS_TOP(vdir) * p->Config.RC_RATIO,
      Parent->Cube.Index.K + IS_FRONT(vdir) * p->Config.RC_RATIO);
		}
  }

  while(Parent->Cube.SubdivisionLevel != subdivision){

//!!!!!!!!!!!!!!
    int idx = (((i - ((i / ratio) * ratio - ((i % ratio) < 0 ? ratio : 0))) < (ratio >> 1))
              ? 0 : 1)
              + ((((j - ((j / ratio) * ratio - ((j % ratio) < 0 ? ratio : 0))) < (ratio >> 1))
              ? 0 : 1) << 1)
              + ((((k - ((k / ratio) * ratio - ((k % ratio) < 0 ? ratio : 0))) < (ratio >> 1))
              ? 0 : 1) << 2);
    ratio = ratio >> 1;


    if(Parent->Cube.Children[idx] == NULL){
//....create child
      TCubes *NewCube = (TCubes *) mycalloc(1, sizeof(TCubesHash));
      NewCube->Cube.Parent = Parent;
      NewCube->Cube.SubdivisionLevel = Parent->Cube.SubdivisionLevel + 1;
      NewCube->Cube.Index.I = (i / ratio) * ratio - ((i % ratio) < 0 ? ratio : 0);
      NewCube->Cube.Index.J = (j / ratio) * ratio - ((j % ratio) < 0 ? ratio : 0);
      NewCube->Cube.Index.K = (k / ratio) * ratio - ((k % ratio) < 0 ? ratio : 0);

      NewCube->Cube.Evaluate = 1;
      NewCube->Cube.Subdivide = 0;

      p->Stats.Cubes.Count++;

      p->ActiveCubes->Next = NewCube;
      p->ActiveCubes = NewCube;

      for(int vdir = 0; vdir < 8; vdir++){
        if(NewCube->Cube.Vertices[vdir] == NULL){
				 	NewCube->Cube.Vertices[vdir] = findCubeVertex(p,
             NewCube->Cube.Index.I + IS_RIGHT(vdir) * ( p->Config.RC_RATIO >> NewCube->Cube.SubdivisionLevel),
             NewCube->Cube.Index.J + IS_TOP(vdir) * ( p->Config.RC_RATIO >> NewCube->Cube.SubdivisionLevel),
             NewCube->Cube.Index.K + IS_FRONT(vdir) * ( p->Config.RC_RATIO >> NewCube->Cube.SubdivisionLevel));
        }
      }

      Parent->Cube.Children[idx] = NewCube;
//....if subdivision not set, then set subdivision and add to list
      if(Parent->Cube.Subdivide != 1){
        Parent->Cube.Subdivide = 1;
        if(Parent->Cube.Evaluate != 1){
          Parent->Next = NULL;
          p->ActiveCubes->Next = Parent;
          p->ActiveCubes = Parent;
        }
      }
    }
    Parent = Parent->Cube.Children[idx];

  }

  return Parent;

}

//---------------------------------------------------------------------------
// return the smallest posible cube in direction
TCubes *findNeighbour(TProcess *p, TCube *Cube, int dir){
}
//---------------------------------------------------------------------------


#pragma package(smart_init)


