CO2GridGenerator Class Reference

#include <co2gridgenerator.h>

List of all members.

Public Member Functions
	CO2GridGenerator ()
	~CO2GridGenerator ()
Static Public Member Functions
static void	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)
static void	subdivided_hyper_rectangle_with_hole (Triangulation< 3 > &tria, const std::vector< unsigned int > &repetitions, const Point3D &p1, const Point3D &p2, VecWellInfo &wells)

---

Detailed Description

Definition at line 18 of file co2gridgenerator.h.

---

Constructor & Destructor Documentation

CO2GridGenerator::CO2GridGenerator

(

)

CO2GridGenerator::~CO2GridGenerator

(

)

---

Member Function Documentation

void CO2GridGenerator::subdivided_hyper_rectangle_with_hole	(	Triangulation< 3 > &	tria,
		const std::vector< unsigned int > &	repetitions,
		const Point3D &	p1,
		const Point3D &	p2,
		VecWellInfo &	wells
	)			[static]

Definition at line 89 of file co2gridgenerator.cpp.

{
  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());  

}

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
	)			[static]

Construct a rectangular mesh and adds a hole in it. Important to note that meshes with holes are not structured meshes in the cense that the distribution of its vertices cannot be seen as a matrix of points. So, the creation a general mesh with holes demands data structures and complexities inherent to unstrucured meshes constructions wich is not our objective to develop. Actually if one day we will used complex geometries is better to use third parties meshes generation code lake gmsh. To avoid the complexities of unstructured meshes we restric the hole generation to just one array of contiguos elements along the y axes. In this way we can use the structured mesh generation code of the library deal.II and then add a hole to it deleting cells. Note that if we want to not have dangling vertices the hole must be of the size of one element such that we dont need to delete any vertices. Because to delete vertices we have to do a renumeration of the vertices indices and update all cells that refer to them. This demands much powerfull data structures tipical of unstructured meshes constructrions alghorithms. The hole is specified give the x,z point of the hole than informing the lenght of the well specifying an interval in the y axes direction [y1-y2]. The well is then specified as a line aligned with the y axes. All cells that contain this line are deleted from the mesh.

Parameters:

	tria	The triangulation to be done
	repetitions	The vector containing the number of discretizations for each dimension. (repetitions.size() == 3)
	p1	The lower point of the rectangular mesh.
	p2	The top point of the rectangular mesh.
	repetitions	Vector os size three containing the discretization in each axis direction
	wx	The X coordinate of the hole (well).
	wz	The z coordinate of the hole (well).
	wy1	The lower cordinate of well in Y.
	wy2	The upper coordinate of well in Y.

Definition at line 16 of file co2gridgenerator.cpp.

{
  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());  

}

---

The documentation for this class was generated from the following files:

• co2gridgenerator.h
• co2gridgenerator.cpp