dune-grid  2.2.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 &e,
75  const FieldVector<ct,dimw>& global) const
76  {
77  // loop over all child Entities
78  const HierarchicIterator end = e.hend(e.level()+1);
79  for( HierarchicIterator it = e.hbegin( e.level()+1 ); it != end; ++it )
80  {
81  FieldVector<ct,dim> local = it->geometry().local(global);
82  if (GenericReferenceElements<double, dim>::general(it->type()).checkInside(local))
83  {
84  // return if we found the leaf, else search through the child entites
85  if( is.contains( *it ) )
86  return EntityPointer( it );
87  else
88  return hFindEntity( *it, global );
89  }
90  }
91  std::ostringstream children;
92  HierarchicIterator it = e.hbegin( e.level()+1 );
93  if(it != end) {
94  children << "{" << formatEntityInformation(*it) << "}";
95  for( ++it; it != end; ++it )
96  children << " {" << formatEntityInformation(*it) << "}";
97  }
98  DUNE_THROW(Exception, "{" << className(*this) << "} Unexpected "
99  "internal Error: none of the children of the entity "
100  "{" << formatEntityInformation(e) << "} contains "
101  "coordinate (" << global << "). Children are: "
102  "[" << children.str() << "].");
103  }
104 
105  public:
109  HierarchicSearch(const Grid & _g, const IS & _is) : g(_g), is(_is) {};
110 
118  EntityPointer findEntity(const FieldVector<ct,dimw>& global) const
119  { return findEntity<All_Partition>(global); }
120 
128  template<PartitionIteratorType partition>
129  EntityPointer findEntity(const FieldVector<ct,dimw>& global) const
130  {
131  typedef typename Grid::template Partition<partition>::LevelGridView
132  LevelGV;
133  const LevelGV &gv = g.template levelView<partition>(0);
134 
136  typedef typename LevelGV::template Codim<0>::Iterator LevelIterator;
137 
138  // loop over macro level
139  LevelIterator it = gv.template begin<0>();
140  LevelIterator end = gv.template end<0>();
141  for (; it != end; ++it)
142  {
143  const Entity &e = *it;
144  const typename Entity::Geometry &geo = e.geometry();
145 
146  FieldVector< ct, dim > local = geo.local( global );
147  if( !GenericReferenceElements< double, dim >::general( geo.type() ).checkInside( local ) )
148  continue;
149 
150  if( (int(dim) != int(dimw)) && ((geo.global( local ) - global).two_norm() > 1e-8) )
151  continue;
152 
153  // return if we found the leaf, else search through the child entites
154  if( is.contains( *it ) )
155  return EntityPointer( it );
156  else
157  return hFindEntity( *it, global );
158  }
159  DUNE_THROW( GridError, "Coordinate " << global << " is outside the grid." );
160  }
161 
162  private:
163  const Grid& g;
164  const IS& is;
165  };
166 
167 } // end namespace Dune
168 
169 #endif // DUNE_GRID_HIERARCHICSEARCH_HH
170