Loading...
Searching...
No Matches
GridAdapter.hpp
/*
Copyright 2010 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_GRIDADAPTER_HEADER_INCLUDED
#define OPM_GRIDADAPTER_HEADER_INCLUDED
#include <opm/grid/UnstructuredGrid.h>
#include <opm/grid/CpGrid.hpp>
#include <stdexcept>
#include <vector>
class GridAdapter
{
public:
template <class Grid>
void init(const Grid& grid)
{
buildTopology(grid);
buildGeometry(grid);
buildGlobalCell(grid);
copyCartDims(grid);
}
UnstructuredGrid* c_grid()
{
return &g_;
}
const UnstructuredGrid* c_grid() const
{
return &g_;
}
// ------ Forwarding the same interface that init() expects ------
//
// This is only done in order to verify that init() works correctly.
//
enum { dimension = 3 }; // This is actually a hack used for testing (dim is a runtime parameter).
struct Vector
{
explicit Vector(const double* source)
{
for (int i = 0; i < dimension; ++i) {
data[i] = source[i];
}
}
double& operator[] (const int ix)
{
return data[ix];
}
double operator[] (const int ix) const
{
return data[ix];
}
double data[dimension];
};
// Topology
int numCells() const
{
return g_.number_of_cells;
}
int numFaces() const
{
return g_.number_of_faces;
}
int numVertices() const
{
return g_.number_of_nodes;
}
int numCellFaces(int cell) const
{
return cell_facepos_[cell + 1] - cell_facepos_[cell];
}
int cellFace(int cell, int local_index) const
{
return cell_faces_[cell_facepos_[cell] + local_index];
}
int faceCell(int face, int local_index) const
{
return face_cells_[2*face + local_index];
}
int numFaceVertices(int face) const
{
return face_nodepos_[face + 1] - face_nodepos_[face];
}
int faceVertex(int face, int local_index) const
{
return face_nodes_[face_nodepos_[face] + local_index];
}
// Geometry
Vector vertexPosition(int vertex) const
{
return Vector(&node_coordinates_[g_.dimensions*vertex]);
}
double faceArea(int face) const
{
return face_areas_[face];
}
Vector faceCentroid(int face) const
{
return Vector(&face_centroids_[g_.dimensions*face]);
}
Vector faceNormal(int face) const
{
Vector fn(&face_normals_[g_.dimensions*face]);
// We must renormalize since the stored normals are
// 'unit normal * face area'.
double invfa = 1.0 / faceArea(face);
for (int i = 0; i < dimension; ++i) {
fn[i] *= invfa;
}
return fn;
}
double cellVolume(int cell) const
{
return cell_volumes_[cell];
}
Vector cellCentroid(int cell) const
{
return Vector(&cell_centroids_[g_.dimensions*cell]);
}
bool operator==(const GridAdapter& other)
{
return face_nodes_ == other.face_nodes_
&& face_nodepos_ == other.face_nodepos_
&& face_cells_ == other.face_cells_
&& cell_faces_ == other.cell_faces_
&& cell_facepos_ == other.cell_facepos_
&& node_coordinates_ == other.node_coordinates_
&& face_centroids_ == other.face_centroids_
&& face_areas_ == other.face_areas_
&& face_normals_ == other.face_normals_
&& cell_centroids_ == other.cell_centroids_
&& cell_volumes_ == other.cell_volumes_;
}
// make a grid which looks periodic but do not have 2 half faces for each
// periodic boundary
void makeQPeriodic(const std::vector<int>& hf_ind,const std::vector<int>& periodic_cells){
for(int i=0; i<int(hf_ind.size());++i){
//std::array<int,2> cells;
int& cell0=face_cells_[2*cell_faces_[ hf_ind[i] ]+0];
int& cell1=face_cells_[2*cell_faces_[ hf_ind[i] ]+1];
assert(periodic_cells[2*i+1]>=0);
if(periodic_cells[2*i+0] == cell0){
assert(cell1==-1);
cell1=periodic_cells[2*i+1];
}else{
assert(periodic_cells[2*i+0] == cell1);
assert(cell0==-1);
cell0=periodic_cells[2*i+1];
}
}
}
private:
UnstructuredGrid g_;
// Topology storage.
std::vector<int> face_nodes_;
std::vector<unsigned> face_nodepos_;
std::vector<int> face_cells_;
std::vector<int> cell_faces_;
std::vector<unsigned> cell_facepos_;
// Geometry storage.
std::vector<double> node_coordinates_;
std::vector<double> face_centroids_;
std::vector<double> face_areas_;
std::vector<double> face_normals_;
std::vector<double> cell_centroids_;
std::vector<double> cell_volumes_;
// The global cell information
std::vector<int> global_cell_;
void buildGlobalCell(const Dune::CpGrid& grid)
{
bool all_active=true;
int old_cell=-1;
global_cell_.resize(grid.numCells());
for(int c=0; c<grid.numCells(); ++c)
{
int new_cell=global_cell_[c]=grid.globalCell()[c];
all_active = all_active && (new_cell==old_cell+1);
old_cell=new_cell;
}
if(all_active){
g_.global_cell=0;
// really release memory
std::vector<int>().swap(global_cell_);
}
else
g_.global_cell=&(global_cell_[0]);
}
template<class Grid>
void buildGlobalCell(const Grid& /*grid*/)
{
g_.global_cell=0;
}
void copyCartDims(const Dune::CpGrid& grid)
{
for(int i=0; i<3; ++i)
g_.cartdims[i] = grid.logicalCartesianSize()[i];
}
template <class Grid>
void copyCartDims(const Grid& /*grid*/)
{}
template <class Grid>
void buildTopology(const Grid& grid)
{
// Face topology.
int num_cells = grid.numCells();
int num_faces = grid.numFaces();
face_nodepos_.resize(num_faces + 1);
int facenodecount = 0;
for (int f = 0; f < num_faces; ++f) {
face_nodepos_[f] = facenodecount;
facenodecount += grid.numFaceVertices(f);
}
face_nodepos_.back() = facenodecount;
face_nodes_.resize(facenodecount);
for (int f = 0; f < num_faces; ++f) {
for (int local = 0; local < grid.numFaceVertices(f); ++local) {
face_nodes_[face_nodepos_[f] + local] = grid.faceVertex(f, local);
}
}
face_cells_.resize(2*num_faces);
for (int f = 0; f < num_faces; ++f) {
face_cells_[2*f] = grid.faceCell(f, 0);
face_cells_[2*f + 1] = grid.faceCell(f, 1);
}
// Cell topology.
int cellfacecount = 0;
cell_facepos_.resize(num_cells + 1);
for (int c = 0; c < num_cells; ++c) {
cell_facepos_[c] = cellfacecount;
cellfacecount += grid.numCellFaces(c);
}
cell_facepos_.back() = cellfacecount;
cell_faces_.resize(cellfacecount);
for (int c = 0; c < num_cells; ++c) {
for (int local = 0; local < grid.numCellFaces(c); ++local) {
cell_faces_[cell_facepos_[c] + local] = grid.cellFace(c, local);
}
}
// Set C grid members.
g_.dimensions = Grid::dimension;
g_.number_of_cells = grid.numCells();
g_.number_of_faces = grid.numFaces();
g_.number_of_nodes = grid.numVertices();
g_.face_nodes = &face_nodes_[0];
g_.face_nodepos = &face_nodepos_[0];
g_.face_cells = &face_cells_[0];
g_.cell_faces = &cell_faces_[0];
g_.cell_facepos = &cell_facepos_[0];
}
template <class Grid>
void buildGeometry(const Grid& grid)
{
// Node geometry.
int num_cells = grid.numCells();
int num_nodes = grid.numVertices();
int num_faces = grid.numFaces();
int dim = Grid::dimension;
node_coordinates_.resize(dim*num_nodes);
for (int n = 0; n < num_nodes; ++n) {
for (int dd = 0; dd < dim; ++dd) {
node_coordinates_[dim*n + dd] = grid.vertexPosition(n)[dd];
}
}
// Face geometry.
face_centroids_.resize(dim*num_faces);
face_areas_.resize(num_faces);
face_normals_.resize(dim*num_faces);
for (int f = 0; f < num_faces; ++f) {
face_areas_[f] = grid.faceArea(f);
for (int dd = 0; dd < dim; ++dd) {
face_centroids_[dim*f + dd] = grid.faceCentroid(f)[dd];
face_normals_[dim*f + dd] = grid.faceNormal(f)[dd]*face_areas_[f];
}
}
// Cell geometry.
cell_centroids_.resize(dim*num_cells);
cell_volumes_.resize(num_cells);
for (int c = 0; c < num_cells; ++c) {
cell_volumes_[c] = grid.cellVolume(c);
for (int dd = 0; dd < dim; ++dd) {
cell_centroids_[dim*c + dd] = grid.cellCentroid(c)[dd];
}
}
// Set C grid members.
g_.node_coordinates = &node_coordinates_[0];
g_.face_centroids = &face_centroids_[0];
g_.face_areas = &face_areas_[0];
g_.face_normals = &face_normals_[0];
g_.cell_centroids = &cell_centroids_[0];
g_.cell_volumes = &cell_volumes_[0];
}
};
#endif // OPM_GRIDADAPTER_HEADER_INCLUDED
UnstructuredGrid.h
Main OPM-Core grid data structure along with helper functions for construction, destruction and readi...
Dune::CpGrid::globalCell
const std::vector< int > & globalCell() const
Retrieve mapping from internal ("compressed") active grid cells to external ("uncompressed") cells.
Definition CpGrid.cpp:680
Dune::CpGrid::logicalCartesianSize
const std::array< int, 3 > & logicalCartesianSize() const
The logical cartesian size of the global grid.
Definition CpGrid.cpp:663
Dune::CpGrid::numCells
int numCells(int level=-1) const
Get the number of cells.
Definition CpGrid.cpp:1147
GridAdapter
Definition GridAdapter.hpp:30
GridAdapter::init
void init(const Grid &grid)
Initialize the grid.
Definition GridAdapter.hpp:37
GridAdapter::c_grid
const UnstructuredGrid * c_grid() const
Access the underlying C grid.
Definition GridAdapter.hpp:50
GridAdapter::Vector
Definition GridAdapter.hpp:63
UnstructuredGrid
Data structure for an unstructured grid, unstructured meaning that any cell may have an arbitrary num...
Definition UnstructuredGrid.h:101
UnstructuredGrid::number_of_cells
int number_of_cells
The number of cells in the grid.
Definition UnstructuredGrid.h:111
UnstructuredGrid::number_of_faces
int number_of_faces
The number of faces in the grid.
Definition UnstructuredGrid.h:113