21 #ifndef EWOMS_POLYHEDRALGRIDCONVERTER_HH
22 #define EWOMS_POLYHEDRALGRIDCONVERTER_HH
25 #include <dune/grid/polyhedralgrid.hh>
27 #include <dune/grid/CpGrid.hpp>
36 struct FaceKey :
public std::pair< int, int >
40 FaceKey(
const int inside,
const int outside )
41 : BaseType( inside < outside ?
std::make_pair(inside,outside) :
std::make_pair(outside,inside) )
47 template <
class Gr
idView,
class CartesianIndexMapper>
48 inline UnstructuredGrid*
50 const CartesianIndexMapper& cartesianIndexMapper,
52 const bool onlyInterior =
true )
55 if( gridView.grid().size( 0 ) == 0 )
59 typedef typename GridView :: template Codim< 0 > :: template Partition<
60 Dune :: All_Partition > :: Iterator Iterator;
61 typedef typename GridView :: template Codim< 0 > :: Entity Element;
62 typedef typename Element :: Geometry ElementGeometry;
63 typedef typename GridView :: IntersectionIterator IntersectionIterator;
64 typedef typename IntersectionIterator :: Intersection Intersection;
65 typedef typename Intersection :: Geometry IntersectionGeometry;
66 typedef typename GridView :: IndexSet IndexSet;
68 typedef typename ElementGeometry :: GlobalCoordinate GlobalCoordinate;
70 const IndexSet& indexSet = gridView.indexSet();
71 const int dimension = GridView::Grid::dimension;
73 const int numCells = indexSet.size( 0 );
74 const int numNodes = indexSet.size( dimension );
76 const int maxNumVerticesPerFace = 4;
77 const int maxNumFacesPerCell = 6;
79 int maxFaceIdx = indexSet.size( 1 );
80 const bool validFaceIndexSet = maxFaceIdx > 0;
84 std::map< FaceKey, int > faceIndexSet;
86 if( ! validFaceIndexSet )
89 const Iterator end = gridView.template end<0, Dune::All_Partition> ();
90 for( Iterator it = gridView.template begin<0, Dune::All_Partition> (); it != end; ++it )
92 const auto& element = *it;
93 const int elIndex = indexSet.index( element );
94 const IntersectionIterator endiit = gridView.iend( element );
95 for( IntersectionIterator iit = gridView.ibegin( element ); iit != endiit; ++iit)
97 const Intersection& intersection = *iit;
101 if( intersection.neighbor() )
103 #if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
104 const Element& neighbor = intersection.outside();
106 auto ep = intersection.outside();
107 const Element& neighbor = *ep;
109 if( ! ( onlyInterior && neighbor.partitionType() != Dune::InteriorEntity ) )
110 nbIndex = indexSet.index( neighbor );
113 FaceKey faceKey( elIndex, nbIndex );
114 if( faceIndexSet.find( faceKey ) == faceIndexSet.end() )
115 faceIndexSet[ faceKey ] = maxFaceIdx++;
119 const int numFaces = maxFaceIdx ;
122 UnstructuredGrid* ug = allocate_grid( dimension, numCells, numFaces,
123 numFaces*maxNumVerticesPerFace,
124 numCells * maxNumFacesPerCell,
127 std::fill( ug->face_cells, ug->face_cells+(numCells * maxNumFacesPerCell), -1 );
129 for(
int d=0; d<dimension; ++d )
130 ug->cartdims[ d ] = cartesianIndexMapper.cartesianDimensions()[ d ];
132 assert( ug->number_of_cells > 0 );
134 if( ! ug->global_cell )
135 ug->global_cell = (
int *) std::malloc( ug->number_of_cells *
sizeof(
int) );
140 const Iterator end = gridView.template end<0, Dune::All_Partition> ();
141 for( Iterator it = gridView.template begin<0, Dune::All_Partition> (); it != end; ++it, ++count )
143 const Element& element = *it;
144 const ElementGeometry geometry = element.geometry();
149 const int elIndex = indexSet.index( element );
150 assert( indexSet.index( element ) == elIndex );
152 const bool isGhost = element.partitionType() != Dune :: InteriorEntity ;
156 assert( count == elIndex );
159 ug->global_cell[ elIndex ] = cartesianIndexMapper.cartesianIndex( elIndex );
163 const GlobalCoordinate center = geometry.center();
164 int idx = elIndex * dimension;
165 for(
int d=0; d<dimension; ++d, ++idx )
166 ug->cell_centroids[ idx ] = center[ d ];
167 ug->cell_volumes[ elIndex ] = geometry.volume();
169 const int vertices = geometry.corners();
170 for(
int vx=0; vx<vertices; ++vx )
172 const GlobalCoordinate vertex = geometry.corner( vx );
173 int idx = indexSet.subIndex( element, vx, dimension ) * dimension;
174 for(
int d=0; d<dimension; ++d, ++idx )
175 ug->node_coordinates[ idx ] = vertex[ d ];
178 ug->cell_facepos[ elIndex ] = cellFace;
180 Dune::GeometryType geomType = element.type();
181 if( geomType.isNone() )
182 geomType = Dune::GeometryType( Dune::GeometryType::cube, dimension );
184 const Dune::ReferenceElement< ctype, dimension > &refElem
185 = Dune::ReferenceElements< ctype, dimension >::general( geomType );
188 const IntersectionIterator endiit = gridView.iend( element );
189 for( IntersectionIterator iit = gridView.ibegin( element ); iit != endiit; ++iit, ++faceCount )
191 const Intersection& intersection = *iit;
192 IntersectionGeometry intersectionGeometry = intersection.geometry();
193 const double faceVol = intersectionGeometry.volume();
195 const int localFace = intersection.indexInInside();
196 const int localFaceIdx = isGhost ? 0 : localFace;
198 int faceIndex = validFaceIndexSet ? indexSet.subIndex( element, localFace, 1 ) : -1;
199 if( ! validFaceIndexSet )
202 if( intersection.neighbor() )
204 #if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
205 const Element& neighbor = intersection.outside();
207 auto ep = intersection.outside();
208 const Element& neighbor = *ep;
210 if( ! ( onlyInterior && neighbor.partitionType() != Dune::InteriorEntity ) )
211 nbIndex = indexSet.index( neighbor );
213 FaceKey faceKey( elIndex, nbIndex );
214 faceIndex = faceIndexSet[ faceKey ];
217 maxFaceIdx = std::max( faceIndex, maxFaceIdx );
219 ug->face_areas[ faceIndex ] = faceVol;
222 const int vxSize = refElem.size( localFace, 1, dimension );
223 int faceIdx = faceIndex * maxNumVerticesPerFace ;
224 ug->face_nodepos[ faceIndex ] = faceIdx;
225 ug->face_nodepos[ faceIndex+1 ] = faceIdx + maxNumVerticesPerFace;
226 for(
int vx=0; vx<vxSize; ++vx, ++faceIdx )
228 const int localVx = refElem.subEntity( localFace, 1, vx, dimension );
229 const int vxIndex = indexSet.subIndex( element, localVx, dimension );
230 ug->face_nodes[ faceIdx ] = vxIndex ;
233 assert( vxSize <= maxNumVerticesPerFace );
234 assert( localFace < maxNumFacesPerCell );
237 ug->cell_faces [ cellFace + localFaceIdx ] = faceIndex;
241 ug->cell_facetag[ cellFace + localFaceIdx ] = localFaceIdx;
244 GlobalCoordinate normal = intersection.centerUnitOuterNormal();
248 if( intersection.neighbor() )
250 #if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
251 const Element& neighbor = intersection.outside();
253 auto ep = intersection.outside();
254 const Element& neighbor = *ep;
258 if( ! ( onlyInterior && neighbor.partitionType() != Dune::InteriorEntity ) )
259 nbIndex = indexSet.index( neighbor );
261 if( nbIndex == -1 || elIndex < nbIndex )
263 ug->face_cells[ 2*faceIndex ] = elIndex;
264 ug->face_cells[ 2*faceIndex + 1 ] = nbIndex;
268 ug->face_cells[ 2*faceIndex ] = nbIndex;
269 ug->face_cells[ 2*faceIndex + 1 ] = elIndex;
276 ug->face_cells[ 2*faceIndex ] = elIndex;
277 ug->face_cells[ 2*faceIndex + 1 ] = -1;
280 const GlobalCoordinate center = intersectionGeometry.center();
282 int idx = faceIndex * dimension;
283 for(
int d=0; d<dimension; ++d, ++idx )
285 ug->face_normals [ idx ] = normal[ d ];
286 ug->face_centroids[ idx ] = center[ d ];
289 if( faceCount > maxNumFacesPerCell )
290 OPM_THROW(std::logic_error,
"DuneGrid only supports conforming hexahedral currently");
291 cellFace += faceCount;
295 ug->cell_facepos[ numCells ] = cellFace;
297 ug->number_of_faces = maxFaceIdx+1;
FaceKey()
Definition: polyhedralgridconverter.hh:39
Definition: polyhedralgridconverter.hh:36
FaceKey(const int inside, const int outside)
Definition: polyhedralgridconverter.hh:40
Definition: baseauxiliarymodule.hh:35
std::pair< int, int > BaseType
Definition: polyhedralgridconverter.hh:38
UnstructuredGrid * dune2UnstructuredGrid(const GridView &gridView, const CartesianIndexMapper &cartesianIndexMapper, const bool faceTags, const bool onlyInterior=true)
Definition: polyhedralgridconverter.hh:49