Skip to content

Bug fixes PolyhedralGrid #513

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 5 commits into from
Feb 3, 2021
Merged

Bug fixes PolyhedralGrid #513

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
a172f54
[Feature][PolyhedralGrid] Make polyhedral grid work in 1d
Feb 1, 2021
558fdbc
[bugfix][PolyhedralGrid] fix geometryType and global coordinates
Feb 2, 2021
e27e001
Give PolyhedralgridIntersection an empty constructor
Feb 2, 2021
129dcb0
Added grid tests to PolyhedralGrid
Feb 2, 2021
52bb4b8
Fixed formating in polyhedralgrid/grid.hh
Feb 2, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
189 changes: 66 additions & 123 deletions opm/grid/polyhedralgrid/grid.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -179,7 +179,8 @@ namespace Dune
static UnstructuredGridPtr
allocateGrid ( std::size_t nCells, std::size_t nFaces, std::size_t nFaceNodes, std::size_t nCellFaces, std::size_t nNodes )
{
UnstructuredGridType *grid = allocate_grid( dim, nCells, nFaces, nFaceNodes, nCellFaces, nNodes );
// Note that we here assign a grid of dimension dimworld in order to obtain global coordinates in the correct dimension
UnstructuredGridType *grid = allocate_grid( dimworld, nCells, nFaces, nFaceNodes, nCellFaces, nNodes );
if( !grid )
DUNE_THROW( GridError, "Unable to allocate grid" );
return UnstructuredGridPtr( grid );
Expand Down Expand Up @@ -972,60 +973,41 @@ namespace Dune
{
const int codim = EntitySeed :: codimension;
const int index = seed.index();
switch (codim)
{
case 0:
{
return cellVertices_[ index ].size();
}
case 1:
{
//return grid_.cell_facepos[ index+1 ] - grid_.cell_facepos[ index ];
return grid_.face_nodepos[ index+1 ] - grid_.face_nodepos[ index ];
}
case dim:
{
return 1;
}
default:
{
return 0;
}
}
if (codim==0)
return cellVertices_[ index ].size();
if (codim==1)
return grid_.face_nodepos[ index+1 ] - grid_.face_nodepos[ index ];
if (codim==dim)
return 1;
return 0;
Copy link
Contributor Author

@rbe051 rbe051 Feb 2, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This change is necessary to make PolyhedralGrid work for 1D grids (then dim = 1, and two cases in the switch are equal). In this case faces and nodes have the same co-dimension. With the current change we always return the face properties for codim=1. If the grid dimension is 1, should we rather return the node properties? I'm not sure if it matters as there is (always?) a one-to-one mapping between faces and nodes in 1D.

}

template <class EntitySeed>
GlobalCoordinate
corner ( const EntitySeed& seed, const int i ) const
{
const int codim = EntitySeed :: codimension;
switch (codim)
if (codim==0)
{
case 0:
{
const int coordIndex = GlobalCoordinate :: dimension * cellVertices_[ seed.index() ][ i ];
return copyToGlobalCoordinate( grid_.node_coordinates + coordIndex );
}
case 1:
{
// for faces we need to swap vertices in 3d since in UnstructuredGrid
// those are ordered counter clockwise, for 2d this does not matter
// TODO: Improve this for performance reasons
const int crners = corners( seed );
const int crner = (crners == 4 && EntitySeed :: dimension == 3 && i > 1 ) ? 5 - i : i;
const int faceVertex = grid_.face_nodes[ grid_.face_nodepos[seed.index() ] + crner ];
return copyToGlobalCoordinate( grid_.node_coordinates + GlobalCoordinate :: dimension * faceVertex );
}
case dim:
{
const int coordIndex = GlobalCoordinate :: dimension * seed.index();
return copyToGlobalCoordinate( grid_.node_coordinates + coordIndex );
}
default:
{
return GlobalCoordinate( 0 );
}
const int coordIndex = GlobalCoordinate :: dimension * cellVertices_[ seed.index() ][ i ];
return copyToGlobalCoordinate( grid_.node_coordinates + coordIndex );
}
if (codim==1)
{
// for faces we need to swap vertices in 3d since in UnstructuredGrid
// those are ordered counter clockwise, for 2d this does not matter
// TODO: Improve this for performance reasons
const int crners = corners( seed );
const int crner = (crners == 4 && EntitySeed :: dimension == 3 && i > 1 ) ? 5 - i : i;
const int faceVertex = grid_.face_nodes[ grid_.face_nodepos[seed.index() ] + crner ];
return copyToGlobalCoordinate( grid_.node_coordinates + GlobalCoordinate :: dimension * faceVertex );
}
if (codim==dim)
{
const int coordIndex = GlobalCoordinate :: dimension * seed.index();
return copyToGlobalCoordinate( grid_.node_coordinates + coordIndex );
}
return GlobalCoordinate( 0 );
}

template <class EntitySeed>
Expand All @@ -1034,25 +1016,19 @@ namespace Dune
const int index = seed.index();
if( seed.codimension == 0 )
{
switch (codim)
{
case 0:
return 1;
case 1:
return grid_.cell_facepos[ index+1 ] - grid_.cell_facepos[ index ];
case dim:
return cellVertices_[ index ].size();
}
if (codim==0)
return 1;
if (codim==1)
return grid_.cell_facepos[ index+1 ] - grid_.cell_facepos[ index ];
if (codim==dim)
return cellVertices_[ index ].size();
}
else if( seed.codimension == 1 )
{
switch (codim)
{
case 1:
return 1;
case dim:
return grid_.face_nodepos[ index+1 ] - grid_.face_nodepos[ index ];
}
if (codim==1)
return 1;
if (codim==dim)
return grid_.face_nodepos[ index+1 ] - grid_.face_nodepos[ index ];
}
else if ( seed.codimension == dim )
{
Expand Down Expand Up @@ -1532,7 +1508,6 @@ namespace Dune

// check face normals
{
typedef Dune::FieldVector< double, dim > Coordinate;
const int faces = grid_.number_of_faces;
for( int face = 0 ; face < faces; ++face )
{
Expand All @@ -1544,41 +1519,41 @@ namespace Dune
if( grid_.face_areas[ face ] < 0 )
std::abort();

Coordinate centerDiff( 0 );
GlobalCoordinate centerDiff( 0 );
if( b >= 0 )
{
for( int d=0; d<dim; ++d )
for( int d=0; d<dimworld; ++d )
{
centerDiff[ d ] = grid_.cell_centroids[ b*dim + d ];
centerDiff[ d ] = grid_.cell_centroids[ b*dimworld + d ];
}
}
else
{
for( int d=0; d<dim; ++d )
for( int d=0; d<dimworld; ++d )
{
centerDiff[ d ] = grid_.face_centroids[ face*dim + d ];
centerDiff[ d ] = grid_.face_centroids[ face*dimworld + d ];
}
}

if( a >= 0 )
{
for( int d=0; d<dim; ++d )
for( int d=0; d<dimworld; ++d )
{
centerDiff[ d ] -= grid_.cell_centroids[ a*dim + d ];
centerDiff[ d ] -= grid_.cell_centroids[ a*dimworld + d ];
}
}
else
{
for( int d=0; d<dim; ++d )
for( int d=0; d<dimworld; ++d )
{
centerDiff[ d ] -= grid_.face_centroids[ face*dim + d ];
centerDiff[ d ] -= grid_.face_centroids[ face*dimworld + d ];
}
}

Coordinate normal( 0 );
for( int d=0; d<dim; ++d )
GlobalCoordinate normal( 0 );
for( int d=0; d<dimworld; ++d )
{
normal[ d ] = grid_.face_normals[ face*dim + d ];
normal[ d ] = grid_.face_normals[ face*dimworld + d ];
}

if( centerDiff.two_norm() < 1e-10 )
Expand All @@ -1593,9 +1568,11 @@ namespace Dune
}
}

// if no face_tag is available we set the reference element based
// on the number of nodes in a cell.
// By default set all types to None. This corresponds to hasPolygon
GeometryType tmp;
tmp = Dune::GeometryTypes::none(dim);
// by default set all types to None
cellGeomTypes_.resize( numCells );
std::fill( cellGeomTypes_.begin(), cellGeomTypes_.end(), tmp );

Expand All @@ -1622,45 +1599,22 @@ namespace Dune
hasPolyhedron = true;
}
}

// if no face_tag is available we assume that no reference element can be
// assigned to the elements
// Propogate the cell geometry type to all codimensions
geomTypes_.resize(dim + 1);
tmp = Dune::GeometryTypes::cube(0);
geomTypes_[ dim ].push_back( tmp );
tmp = Dune::GeometryTypes::cube(1);
geomTypes_[ dim-1 ].push_back( tmp );

if( hasSimplex )
{
tmp = Dune::GeometryTypes::simplex(dim);
geomTypes_[ 0 ].push_back( tmp );
}

if( hasCube )
{
tmp = Dune::GeometryTypes::cube(dim);
geomTypes_[ 0 ].push_back( tmp );
}

if( hasPolyhedron )
{
tmp = Dune::GeometryTypes::none(dim);
geomTypes_[ 0 ].push_back( tmp );
}

if( dim > 2 )
for (int codim = 0; codim <= dim; ++codim)
{
for( const auto& elemType : geomTypes_[ 0 ] )
{
if( elemType.isSimplex() )
tmp = Dune::GeometryTypes::simplex(2);
else if ( elemType.isCube() )
tmp = Dune::GeometryTypes::cube(2);
if( hasSimplex ){
tmp = Dune::GeometryTypes::simplex(dim - codim);
geomTypes_[ codim ].push_back( tmp );
}
else if ( hasCube ){
tmp = Dune::GeometryTypes::cube(dim - codim);
geomTypes_[ codim ].push_back( tmp );}
else if (hasPolyhedron){
tmp = Dune::GeometryTypes::none(dim - codim);
geomTypes_[ codim ].push_back( tmp );}
Copy link
Member

@blattms blattms Feb 2, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nitpick: Formatting is messed up (not enough indentation and } on the same line) which makes this hard to read.

else
tmp = Dune::GeometryTypes::none(2);
geomTypes_[ 1 ].push_back( tmp );
}
OPM_THROW(std::runtime_error, "Grid error, unkown geometry type.");
}

} // end else of ( grid_.cell_facetag )
Expand Down Expand Up @@ -1691,17 +1645,6 @@ namespace Dune
}
}
}

/*
for( int i=0; i<= dim ; ++ i )
{
for( const auto& geomType : geomTypes_[ i ] )
{
std::cout << "Codim " << i << " type = " << geomType << std::endl;
}
}
*/
//print( std::cout, grid_ );
}

void print( std::ostream& out, const UnstructuredGridType& grid ) const
Expand Down
6 changes: 6 additions & 0 deletions opm/grid/polyhedralgrid/intersection.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,12 @@ namespace Dune
intersectionIdx_( -1 )
{}

PolyhedralGridIntersection ( )
: data_( ),
seed_(),
intersectionIdx_( -1 )
{}

PolyhedralGridIntersection ( ExtraData data, const EntitySeed& seed, const int intersectionIdx )
: data_( data ),
seed_( seed ),
Expand Down