co2gridgenerator.cpp

Go to the documentation of this file.

#include "co2gridgenerator.h"
#include <limits.h>


void CO2GridGenerator::subdivided_hyper_rectangle_with_hole (Triangulation<3> &tria, const std::vector<unsigned int> &  repetitions,const Point3D  &p1,const Point3D &p2, const double wx, double wz, double wy1, double wy2) 
{
  assert(repetitions.size() == 3);
  assert(p1[0] < p2[0] && p1[1] < p2[1] && p1[2] < p2[2]);
  assert(wy2 > wy1);

                                   // then check that all repetitions
                                   // are >= 1, and calculate deltas
                                   // convert repetitions from double
                                   // to int by taking the ceiling.
  double delta[3];
  for (unsigned int i=0; i<3; ++i)
    {
      assert (repetitions[i] >= 1);
      delta[i] = (p2[i]-p1[i])/repetitions[i];
    }
 
  // then generate the necessary
  // points
  std::vector<Point3D > points;
  for (unsigned int z=0; z<=repetitions[2]; ++z)
    for (unsigned int y=0; y<=repetitions[1]; ++y)
      for (unsigned int x=0; x<=repetitions[0]; ++x)
      {
        points.push_back (Point3D    (p1[0]+x*delta[0],
                                      p1[1]+y*delta[1],
                                      p1[2]+z*delta[2]));
      }
  // next create the cells
  // Prepare cell data
  std::vector<CellData<3> > cells;
  CellData<3> cell;
  const unsigned int n_x  = (repetitions[0]+1);
  const unsigned int n_xy = (repetitions[0]+1)*(repetitions[1]+1);
        
  cells.reserve (repetitions[2]*repetitions[1]*repetitions[0]);
  for (unsigned int z=0; z<repetitions[2]; ++z)
    for (unsigned int y=0; y<repetitions[1]; ++y)
      for (unsigned int x=0; x<repetitions[0]; ++x)
      {
        cell.vertices[0] = z*n_xy + y*n_x + x;
        cell.vertices[1] = z*n_xy + y*n_x + x+1;
        cell.vertices[2] = z*n_xy + (y+1)*n_x + x;
        cell.vertices[3] = z*n_xy + (y+1)*n_x + x+1;
        cell.vertices[4] = (z+1)*n_xy + y*n_x + x;
        cell.vertices[5] = (z+1)*n_xy + y*n_x + x+1;
        cell.vertices[6] = (z+1)*n_xy + (y+1)*n_x + x;
        cell.vertices[7] = (z+1)*n_xy + (y+1)*n_x + x+1;
        cell.material_id = 0;

        //If cell is in the hole, dont insert it in the array of cells.
        Point3D vertex_000 = points[cell.vertices[VERTEX_000]];
        Point3D vertex_111 = points[cell.vertices[VERTEX_111]];

        if ( wx > vertex_000[0] && wx < vertex_111[0] &&
             wz > vertex_000[2] && wz < vertex_111[2] &&
             !( vertex_000[1] >=  wy2 || vertex_111[1] <= wy1)) 
        {
          //Ok the cell contains at least on point in the line X=wx,y in [wy1,wy2], Z = wz
          //So dont add the cell to the array of cells
          continue;
        }
        else 
        {
          cells.push_back(cell);
        }
      }

  tria.create_triangulation (points, cells, SubCellData());  

}


void CO2GridGenerator::subdivided_hyper_rectangle_with_hole (Triangulation<3> &tria, const std::vector<unsigned int> &  repetitions,const Point3D  &p1,const Point3D &p2, VecWellInfo &wells) 
{
  assert(repetitions.size() == 3);
  assert(p1[0] < p2[0] && p1[1] < p2[1] && p1[2] < p2[2]);

                                   // then check that all repetitions
                                   // are >= 1, and calculate deltas
                                   // convert repetitions from double
                                   // to int by taking the ceiling.
  

  unsigned  nVertices=(repetitions[0]+1)*(repetitions[1]+1)*(repetitions[2]+1);
         
  double delta[3];
  unsigned INVALID_INDEX = UINT_MAX;
  
  for (unsigned int i=0; i<3; ++i)
    {
      assert (repetitions[i] >= 1);
      delta[i] = (p2[i]-p1[i])/repetitions[i];
    }

  
  // then generate the necessary
  // points
  std::vector<Point3D > points;
  points.reserve(nVertices);
  std::vector<unsigned> pointsIndex(nVertices);
  for (unsigned int z=0; z<=repetitions[2]; ++z)
    for (unsigned int y=0; y<=repetitions[1]; ++y)
      for (unsigned int x=0; x<=repetitions[0]; ++x)
      {
        points.push_back (Point3D    (p1[0]+x*delta[0],
                                      p1[1]+y*delta[1],
                                      p1[2]+z*delta[2]));
      }

  //Verify which vertices are valid (not inside wells)
  //and enumerate all vertices not counting the invalid
  //ones. Put the indices in the vector pointsIndex.
  //pointsIndex[i] == indice of ith vertice in the vector "points"
  //pointsIndex[i] == INVALID_INDEX if the vertice is invalid.
  unsigned index=0;
  Point3D DX_4(delta[0]/4.0,delta[1]/4.0,delta[2]/4.0);
  bool bIsValid;
  for(unsigned i=0;i<nVertices;i++)
  {
    bIsValid=true;
    for(unsigned j=0;j<wells.size();j++)
    {
      WellInfo well = wells[j];
      well.P+=DX_4;
      well.Q-=DX_4;
      if (well.isPointInWell(points[i]))
      {
        bIsValid=false;
        break;
      }
    }
    if (bIsValid)
      pointsIndex[i]=index++;
    else
    {
      pointsIndex[i]=INVALID_INDEX;
      //printf("vertice ripped out %d\n",i);
    }
  }
  

  // next create the cells
  // Prepare cell data
  std::vector<CellData<3> > cells;
  CellData<3> cell;
  const unsigned int n_x  = (repetitions[0]+1);
  const unsigned int n_xy = (repetitions[0]+1)*(repetitions[1]+1);
        
  cells.reserve (repetitions[2]*repetitions[1]*repetitions[0]);
  for (unsigned int z=0; z<repetitions[2]; ++z)
    for (unsigned int y=0; y<repetitions[1]; ++y)
      for (unsigned int x=0; x<repetitions[0]; ++x)
      {
        cell.vertices[0] = z*n_xy + y*n_x + x;
        cell.vertices[1] = z*n_xy + y*n_x + x+1;
        cell.vertices[2] = z*n_xy + (y+1)*n_x + x;
        cell.vertices[3] = z*n_xy + (y+1)*n_x + x+1;
        cell.vertices[4] = (z+1)*n_xy + y*n_x + x;
        cell.vertices[5] = (z+1)*n_xy + y*n_x + x+1;
        cell.vertices[6] = (z+1)*n_xy + (y+1)*n_x + x;
        cell.vertices[7] = (z+1)*n_xy + (y+1)*n_x + x+1;
        cell.material_id = 0;

        //If cell is in the hole, dont insert it in the array of cells.
        //First find the cell's barycenter
        Point3D BC  = points[cell.vertices[VERTEX_000]];
        BC += points[cell.vertices[VERTEX_111]];
        BC/=2.0;

        //Find out if the cell is inside one of the wells
        bool isCellInWell=false;
        for (unsigned w=0;w<wells.size();w++)
        {
          if (wells[w].isPointInWell(BC))
          {
            isCellInWell=true;
            break;
          }
        }
        if (!isCellInWell)
        {
          //Insert it in the vector cells
          //but first map the vertices to the right indices
          for (unsigned k=0;k<8;k++)
          {
            cell.vertices[k]=pointsIndex[cell.vertices[k]];
            assert(cell.vertices[k] != INVALID_INDEX);
          }
          cells.push_back(cell);
        }
      }
  //Now adjust the points vector containning the vertices
  std::vector<Point3D>::iterator it = points.begin();
  for (unsigned i=0;i<nVertices;i++)
  {
    if (pointsIndex[i] != INVALID_INDEX)
      *(it++) = points[i];
  }
  if(it != points.end())
    points.erase(it,points.end());

  /*

  for (unsigned i=0;i<points.size();i++)
  {
    printf("%u) %g %g %g\n" ,i,points[i](0),points[i](1),points[i](2));
    
  }
  */
  
  tria.create_triangulation (points, cells, SubCellData());  

}

---