My Project
gridfactory.hh
1// -*- mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=2 sw=2 sts=2:
3#ifndef DUNE_POLYHEDRALGRID_GRIDFACTORY_HH
4#define DUNE_POLYHEDRALGRID_GRIDFACTORY_HH
5
6#include <algorithm>
7#include <numeric>
8
9#include <dune/common/typetraits.hh>
10#include <dune/common/version.hh>
11
12#include <dune/grid/common/gridfactory.hh>
13#include <opm/grid/polyhedralgrid/grid.hh>
14
15namespace Dune
16{
17
18
19 // GridFactory for PolyhedralGrid
20 // ---------------------------------
21
22 template< int dim, int dimworld, class coord_t >
23 class GridFactory< PolyhedralGrid< dim, dimworld, coord_t > >
24 : public GridFactoryInterface< PolyhedralGrid< dim, dimworld, coord_t > >
25 {
26 public:
28
29 const static int dimension = Grid::dimension;
30 const static int dimensionworld = Grid::dimensionworld;
31 typedef typename Grid::ctype ctype;
32
33 typedef MPIHelper::MPICommunicator MPICommunicatorType;
34 typedef typename Grid::template Codim<0>::Entity Element;
35 typedef typename Grid::template Codim<dimension>::Entity Vertex;
36
37 typedef Dune::FieldVector<ctype,dimensionworld> CoordinateType;
38 typedef CoordinateType Coordinate;
39
40#if DUNE_VERSION_GTE(DUNE_GRID, 2, 7)
41#if DUNE_VERSION_LT(DUNE_GRID, 2, 8)
42 typedef ToUniquePtr<Grid> UniquePtrType;
43#else
44 using UniquePtrType = std::unique_ptr<Grid>;
45#endif
46#else // #if DUNE_VERSION_GTE(DUNE_GRID, 2, 7)
47 typedef Grid* UniquePtrType;
48#endif // #else // #if DUNE_VERSION_GTE(DUNE_GRID, 2, 7)
49
50
52 explicit GridFactory ( const MPICommunicatorType& = MPIHelper::getCommunicator() )
53 : nodes_(),
54 faces_(),
55 cells_()
56 {}
57
58 virtual void insertVertex(const CoordinateType& pos)
59 {
60 nodes_.push_back( pos );
61 }
62
72 virtual void insertElement(const GeometryType& type,
73 const std::vector<unsigned int>& items)
74 {
75 if( type.isNone() )
76 {
77 // copy into vector of integers
78 std::vector< int > numbers( items.size() );
79 std::copy( items.begin(), items.end(), numbers.begin() );
80
81 if( type.dim() == dimension-1 )
82 {
83 faces_.push_back( numbers );
84 }
85 else if( type.dim() == dimension )
86 {
87 // note vertices holds the face
88 // numbers in this case
89 cells_.push_back( numbers );
90 }
91 else
92 {
93 DUNE_THROW(Dune::NotImplemented,"insertElement not implemented for type " << type );
94 }
95 }
96 else // use ReferenceElement to insert faces
97 {
98
99
100 }
101 }
102
103 void insertBoundarySegment(const std::vector<unsigned int>&)
104 {
105 DUNE_THROW(NotImplemented,"yet");
106 }
107
108 UniquePtrType createGrid()
109 {
110 std::vector< CoordinateType >& nodes = nodes_;
111 std::vector< std::vector< int > >& faces = faces_;
112 std::vector< std::vector< int > >& cells = cells_;
113
114 if( cells.empty() )
115 {
116 DUNE_THROW( GridError, "No cells found for PolyhedralGrid" );
117 }
118
119 const auto sumSize = [] ( std::size_t s, const std::vector< int > &v ) { return s + v.size(); };
120 const std::size_t numFaceNodes = std::accumulate( faces.begin(), faces.end(), std::size_t( 0 ), sumSize );
121 const std::size_t numCellFaces = std::accumulate( cells.begin(), cells.end(), std::size_t( 0 ), sumSize );
122
123 typename Grid::UnstructuredGridPtr ug =
124 Grid::allocateGrid( cells.size(), faces.size(), numFaceNodes, numCellFaces, nodes.size() );
125
126 // copy faces
127 {
128#ifndef NDEBUG
129 std::map< std::vector< int >, std::vector< int > > faceMap;
130#endif
131
132 const int nFaces = faces.size();
133 // set all face_cells values to -2 as default
134 std::fill( ug->face_cells, ug->face_cells + 2*nFaces, -1 );
135
136 int facepos = 0;
137 std::vector< int > faceVertices;
138 faceVertices.reserve( 30 );
139 for( int face = 0; face < nFaces; ++face )
140 {
141 //std::cout << "face " << face << ": ";
142 faceVertices.clear();
143 ug->face_nodepos[ face ] = facepos;
144 const int nVertices = faces[ face ].size();
145 for( int vx = 0; vx < nVertices; ++vx, ++facepos )
146 {
147 //std::cout << " " << faces[ face ][ vx ];
148 ug->face_nodes[ facepos ] = faces[ face ][ vx ];
149 faceVertices.push_back( faces[ face ][ vx ] );
150 }
151 //std::cout << std::endl;
152
153#ifndef NDEBUG
154 // sort vertices
155 std::sort( faceVertices.begin(), faceVertices.end() );
156 // make sure each face only exists once
157 faceMap[ faceVertices ].push_back( face );
158 assert( faceMap[ faceVertices ].size() == 1 );
159#endif
160 }
161 ug->face_nodepos[ nFaces ] = facepos ;
162 }
163
164 // copy cells
165 {
166 const int nCells = cells.size();
167 int cellpos = 0;
168 for( int cell = 0; cell < nCells; ++cell )
169 {
170 //std::cout << "Cell " << cell << ": ";
171 ug->cell_facepos[ cell ] = cellpos;
172 const int nFaces = cells[ cell ].size();
173 for( int f = 0; f < nFaces; ++f, ++cellpos )
174 {
175 const int face = cells[ cell ][ f ];
176 // std::cout << " " << face ;
177 ug->cell_faces[ cellpos ] = face;
178
179 // TODO find cells for each face
180 if( ug->face_cells[ 2*face ] == -1 )
181 {
182 ug->face_cells[ 2*face ] = cell;
183 }
184 else // if ( ug->face_cells[ 2*face+1 ] == -1 )
185 {
186 //assert( ug->face_cells[ 2*face+1 ] == -1 );
187 ug->face_cells[ 2*face+1 ] = cell;
188 }
189 }
190 //std::cout << std::endl;
191 }
192 ug->cell_facepos[ nCells ] = cellpos ;
193 }
194
195 // copy node coordinates
196 {
197 const int nNodes = nodes.size();
198 int nodepos = 0;
199 for( int vx = 0 ; vx < nNodes; ++vx )
200 {
201 for( int d=0; d<dim; ++d, ++nodepos )
202 ug->node_coordinates[ nodepos ] = nodes[ vx ][ d ];
203 }
204 }
205
206 /*
207 for( int i=0; i<int(faces.size() ); ++i)
208 {
209 std::cout << "face "<< i<< " connects to " << ug->face_cells[ 2*i ] << " " <<
210 ug->face_cells[ 2*i+1] << std::endl;
211 }
212 */
213
214 // free cell face tag since it's not a cartesian grid
215 if( ug->cell_facetag )
216 {
217 std::free( ug->cell_facetag );
218 ug->cell_facetag = nullptr ;
219 for( int i=0; i<3; ++i ) ug->cartdims[ i ] = 0;
220 }
221
222 // compute geometric quantities like cell volume and face normals
223 Grid::computeGeometry( ug );
224
225 // check normal direction
226 {
227 for( int face = 0 ; face < ug->number_of_faces; ++face )
228 {
229 const int a = ug->face_cells[ 2*face ];
230 const int b = ug->face_cells[ 2*face + 1 ];
231 if( a < 0 || b < 0 )
232 continue ;
233
234 Coordinate centerDiff( 0 );
235 Coordinate normal( 0 );
236 //std::cout << "Cell center " << a << " " << b << std::endl;
237 for( int d=0; d<dim; ++d )
238 {
239 //std::cout << ug->cell_centroids[ a*dim + d ] << " " << ug->cell_centroids[ b*dim + d ] << std::endl;
240 centerDiff[ d ] = ug->cell_centroids[ b*dim + d ] - ug->cell_centroids[ a*dim + d ];
241 normal[ d ] = ug->face_normals[ face*dim + d ];
242 }
243
244 // if diff and normal point in different direction, flip faces
245 if( centerDiff * normal > 0 )
246 {
247 ug->face_cells[ 2*face ] = b;
248 ug->face_cells[ 2*face + 1 ] = a;
249 }
250 }
251 }
252
253 return UniquePtrType( new Grid( std::move( ug ) ));
254 }
255
256 protected:
257 std::vector< CoordinateType > nodes_;
258 std::vector< std::vector< int > > faces_;
259 std::vector< std::vector< int > > cells_;
260 };
261
262} // namespace Dune
263
264#endif // #ifndef DUNE_POLYHEDRALGRID_DGFPARSER_HH
virtual void insertElement(const GeometryType &type, const std::vector< unsigned int > &items)
Insert an element into the coarse grid.
Definition: gridfactory.hh:72
GridFactory(const MPICommunicatorType &=MPIHelper::getCommunicator())
Default constructor.
Definition: gridfactory.hh:52
identical grid wrapper
Definition: grid.hh:163
Traits::ctype ctype
type of vector coordinates (e.g., double)
Definition: grid.hh:312
Copyright 2019 Equinor AS.
Definition: CartesianIndexMapper.hpp:10