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 
15 namespace 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:158
Traits::ctype ctype
type of vector coordinates (e.g., double)
Definition: grid.hh:307
Copyright 2019 Equinor AS.
Definition: CartesianIndexMapper.hpp:10