Loading...
Searching...
No Matches
gridfactory.hh
// -*- mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=2 sw=2 sts=2:
#ifndef DUNE_POLYHEDRALGRID_GRIDFACTORY_HH
#define DUNE_POLYHEDRALGRID_GRIDFACTORY_HH
#include <algorithm>
#include <numeric>
#include <dune/common/typetraits.hh>
#include <dune/common/version.hh>
#include <dune/grid/common/gridfactory.hh>
#include <opm/grid/polyhedralgrid/grid.hh>
namespace Dune
{
// GridFactory for PolyhedralGrid
// ---------------------------------
template< int dim, int dimworld, class coord_t >
class GridFactory< PolyhedralGrid< dim, dimworld, coord_t > >
: public GridFactoryInterface< PolyhedralGrid< dim, dimworld, coord_t > >
{
public:
typedef PolyhedralGrid< dim, dimworld, coord_t > Grid;
const static int dimension = Grid::dimension;
const static int dimensionworld = Grid::dimensionworld;
typedef typename Grid::ctype ctype;
typedef MPIHelper::MPICommunicator MPICommunicatorType;
typedef typename Grid::template Codim<0>::Entity Element;
typedef typename Grid::template Codim<dimension>::Entity Vertex;
typedef Dune::FieldVector<ctype,dimensionworld> CoordinateType;
typedef CoordinateType Coordinate;
using UniquePtrType = std::unique_ptr<Grid>;
explicit GridFactory ( const MPICommunicatorType& = MPIHelper::getCommunicator() )
: nodes_(),
faces_(),
cells_()
{}
virtual void insertVertex(const CoordinateType& pos)
{
nodes_.push_back( pos );
}
virtual void insertElement(const GeometryType& type,
const std::vector<unsigned int>& items)
{
if( type.isNone() )
{
// copy into vector of integers
std::vector< int > numbers( items.size() );
std::copy( items.begin(), items.end(), numbers.begin() );
if( type.dim() == dimension-1 )
{
faces_.push_back( numbers );
}
else if( type.dim() == dimension )
{
// note vertices holds the face
// numbers in this case
cells_.push_back( numbers );
}
else
{
DUNE_THROW(Dune::NotImplemented,"insertElement not implemented for type " << type );
}
}
else // use ReferenceElement to insert faces
{
}
}
void insertBoundarySegment(const std::vector<unsigned int>&)
{
DUNE_THROW(NotImplemented,"yet");
}
UniquePtrType createGrid()
{
std::vector< CoordinateType >& nodes = nodes_;
std::vector< std::vector< int > >& faces = faces_;
std::vector< std::vector< int > >& cells = cells_;
if( cells.empty() )
{
DUNE_THROW( GridError, "No cells found for PolyhedralGrid" );
}
const auto sumSize = [] ( std::size_t s, const std::vector< int > &v ) { return s + v.size(); };
const std::size_t numFaceNodes = std::accumulate( faces.begin(), faces.end(), std::size_t( 0 ), sumSize );
const std::size_t numCellFaces = std::accumulate( cells.begin(), cells.end(), std::size_t( 0 ), sumSize );
typename Grid::UnstructuredGridPtr ug =
Grid::allocateGrid( cells.size(), faces.size(), numFaceNodes, numCellFaces, nodes.size() );
// copy faces
{
#ifndef NDEBUG
std::map< std::vector< int >, std::vector< int > > faceMap;
#endif
const int nFaces = faces.size();
// set all face_cells values to -2 as default
std::fill( ug->face_cells, ug->face_cells + 2*nFaces, -1 );
int facepos = 0;
std::vector< int > faceVertices;
faceVertices.reserve( 30 );
for( int face = 0; face < nFaces; ++face )
{
//std::cout << "face " << face << ": ";
faceVertices.clear();
ug->face_nodepos[ face ] = facepos;
const int nVertices = faces[ face ].size();
for( int vx = 0; vx < nVertices; ++vx, ++facepos )
{
//std::cout << " " << faces[ face ][ vx ];
ug->face_nodes[ facepos ] = faces[ face ][ vx ];
faceVertices.push_back( faces[ face ][ vx ] );
}
//std::cout << std::endl;
#ifndef NDEBUG
// sort vertices
std::sort( faceVertices.begin(), faceVertices.end() );
// make sure each face only exists once
faceMap[ faceVertices ].push_back( face );
assert( faceMap[ faceVertices ].size() == 1 );
#endif
}
ug->face_nodepos[ nFaces ] = facepos ;
}
// copy cells
{
const int nCells = cells.size();
int cellpos = 0;
for( int cell = 0; cell < nCells; ++cell )
{
//std::cout << "Cell " << cell << ": ";
ug->cell_facepos[ cell ] = cellpos;
const int nFaces = cells[ cell ].size();
for( int f = 0; f < nFaces; ++f, ++cellpos )
{
const int face = cells[ cell ][ f ];
// std::cout << " " << face ;
ug->cell_faces[ cellpos ] = face;
// TODO find cells for each face
if( ug->face_cells[ 2*face ] == -1 )
{
ug->face_cells[ 2*face ] = cell;
}
else // if ( ug->face_cells[ 2*face+1 ] == -1 )
{
//assert( ug->face_cells[ 2*face+1 ] == -1 );
ug->face_cells[ 2*face+1 ] = cell;
}
}
//std::cout << std::endl;
}
ug->cell_facepos[ nCells ] = cellpos ;
}
// copy node coordinates
{
const int nNodes = nodes.size();
int nodepos = 0;
for( int vx = 0 ; vx < nNodes; ++vx )
{
for( int d=0; d<dim; ++d, ++nodepos )
ug->node_coordinates[ nodepos ] = nodes[ vx ][ d ];
}
}
/*
for( int i=0; i<int(faces.size() ); ++i)
{
std::cout << "face "<< i<< " connects to " << ug->face_cells[ 2*i ] << " " <<
ug->face_cells[ 2*i+1] << std::endl;
}
*/
// free cell face tag since it's not a cartesian grid
if( ug->cell_facetag )
{
std::free( ug->cell_facetag );
ug->cell_facetag = nullptr ;
for( int i=0; i<3; ++i ) ug->cartdims[ i ] = 0;
}
// compute geometric quantities like cell volume and face normals
Grid::computeGeometry( ug );
// check normal direction
{
for( int face = 0 ; face < ug->number_of_faces; ++face )
{
const int a = ug->face_cells[ 2*face ];
const int b = ug->face_cells[ 2*face + 1 ];
if( a < 0 || b < 0 )
continue ;
Coordinate centerDiff( 0 );
Coordinate normal( 0 );
//std::cout << "Cell center " << a << " " << b << std::endl;
for( int d=0; d<dim; ++d )
{
//std::cout << ug->cell_centroids[ a*dim + d ] << " " << ug->cell_centroids[ b*dim + d ] << std::endl;
centerDiff[ d ] = ug->cell_centroids[ b*dim + d ] - ug->cell_centroids[ a*dim + d ];
normal[ d ] = ug->face_normals[ face*dim + d ];
}
// if diff and normal point in different direction, flip faces
if( centerDiff * normal > 0 )
{
ug->face_cells[ 2*face ] = b;
ug->face_cells[ 2*face + 1 ] = a;
}
}
}
return UniquePtrType( new Grid( std::move( ug ) ));
}
protected:
std::vector< CoordinateType > nodes_;
std::vector< std::vector< int > > faces_;
std::vector< std::vector< int > > cells_;
};
} // namespace Dune
#endif // #ifndef DUNE_POLYHEDRALGRID_DGFPARSER_HH
Dune::GridFactory< PolyhedralGrid< dim, dimworld, coord_t > >::insertElement
virtual void insertElement(const GeometryType &type, const std::vector< unsigned int > &items)
Insert an element into the coarse grid.
Definition gridfactory.hh:63
Dune::GridFactory< PolyhedralGrid< dim, dimworld, coord_t > >::GridFactory
GridFactory(const MPICommunicatorType &=MPIHelper::getCommunicator())
Default constructor.
Definition gridfactory.hh:43
Dune::PolyhedralGrid
identical grid wrapper
Definition grid.hh:159
Dune::PolyhedralGrid::ctype
Traits::ctype ctype
type of vector coordinates (e.g., double)
Definition grid.hh:308
Dune
The namespace Dune is the main namespace for all Dune code.
Definition CartesianIndexMapper.hpp:10