dune-grid  2.3.1
hierarchicsearch.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 
4 #ifndef DUNE_GRID_HIERARCHICSEARCH_HH
5 #define DUNE_GRID_HIERARCHICSEARCH_HH
6 
13 #include <cstddef>
14 #include <sstream>
15 #include <string>
16 
17 #include <dune/common/classname.hh>
18 #include <dune/common/exceptions.hh>
19 #include <dune/common/fvector.hh>
20 
21 #include <dune/grid/common/grid.hh>
23 
24 namespace Dune
25 {
26 
30  template<class Grid, class IS>
32  {
34  enum {dim=Grid::dimension};
35 
37  enum {dimw=Grid::dimensionworld};
38 
40  typedef typename Grid::ctype ct;
41 
43  typedef typename Grid::template Codim<0>::Entity Entity;
44 
46  typedef typename Grid::template Codim<0>::EntityPointer EntityPointer;
47 
49  typedef typename Grid::HierarchicIterator HierarchicIterator;
50 
51  static std::string formatEntityInformation ( const Entity &e ) {
52  const typename Entity::Geometry &geo = e.geometry();
53  std::ostringstream info;
54  info << "level=" << e.level() << " "
55  << "partition=" << e.partitionType() << " "
56  << "center=(" << geo.center() << ") "
57  << "corners=[(" << geo.corner(0) << ")";
58  for(int i = 1; i < geo.corners(); ++i)
59  info << " (" << e.geometry().corner(i) << ")";
60  info << "]";
61  return info.str();
62  }
63 
74  EntityPointer hFindEntity ( const Entity &entity,
75  const FieldVector<ct,dimw>& global) const
76  {
77  // type of element geometry
78  typedef typename Entity::Geometry Geometry;
79  // type of local coordinate
80  typedef typename Geometry::LocalCoordinate LocalCoordinate;
81 
82  const int childLevel = entity.level()+1 ;
83  // loop over all child Entities
84  const HierarchicIterator end = entity.hend( childLevel );
85  for( HierarchicIterator it = entity.hbegin( childLevel ); it != end; ++it )
86  {
87  const Entity &child = *it;
88  const Geometry &geo = child.geometry();
89 
90  LocalCoordinate local = geo.local(global);
91  if (ReferenceElements<double, dim>::general( child.type() ).checkInside(local))
92  {
93  // return if we found the leaf, else search through the child entites
94  if( indexSet_.contains( child ) )
95  return EntityPointer( child );
96  else
97  return hFindEntity( child, global );
98  }
99  }
100  std::ostringstream children;
101  HierarchicIterator it = entity.hbegin( childLevel );
102  if(it != end) {
103  children << "{" << formatEntityInformation(*it) << "}";
104  for( ++it; it != end; ++it )
105  children << " {" << formatEntityInformation(*it) << "}";
106  }
107  DUNE_THROW(Exception, "{" << className(*this) << "} Unexpected "
108  "internal Error: none of the children of the entity "
109  "{" << formatEntityInformation(entity) << "} contains "
110  "coordinate (" << global << "). Children are: "
111  "[" << children.str() << "].");
112  }
113 
114  public:
118  HierarchicSearch(const Grid & g, const IS & is) : grid_(g), indexSet_(is) {}
119 
127  EntityPointer findEntity(const FieldVector<ct,dimw>& global) const
128  { return findEntity<All_Partition>(global); }
129 
137  template<PartitionIteratorType partition>
138  EntityPointer findEntity(const FieldVector<ct,dimw>& global) const
139  {
140  typedef typename Grid::template Partition<partition>::LevelGridView
141  LevelGV;
142  const LevelGV &gv = grid_.template levelGridView<partition>(0);
143 
145  typedef typename LevelGV::template Codim<0>::Iterator LevelIterator;
146 
147  // type of element geometry
148  typedef typename Entity::Geometry Geometry;
149  // type of local coordinate
150  typedef typename Geometry::LocalCoordinate LocalCoordinate;
151 
152  // loop over macro level
153  LevelIterator it = gv.template begin<0>();
154  LevelIterator end = gv.template end<0>();
155  for (; it != end; ++it)
156  {
157  const Entity &entity = *it;
158  const Geometry &geo = entity.geometry();
159 
160  LocalCoordinate local = geo.local( global );
161  if( !ReferenceElements< double, dim >::general( geo.type() ).checkInside( local ) )
162  continue;
163 
164  if( (int(dim) != int(dimw)) && ((geo.global( local ) - global).two_norm() > 1e-8) )
165  continue;
166 
167  // return if we found the leaf, else search through the child entites
168  if( indexSet_.contains( entity ) )
169  return EntityPointer( entity );
170  else
171  return hFindEntity( entity, global );
172  }
173  DUNE_THROW( GridError, "Coordinate " << global << " is outside the grid." );
174  }
175 
176  private:
177  const Grid& grid_;
178  const IS& indexSet_;
179  };
180 
181 } // end namespace Dune
182 
183 #endif // DUNE_GRID_HIERARCHICSEARCH_HH