dune-grid  2.2.0
hierarchicsearch.hh
Go to the documentation of this file.
00001 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
00002 // vi: set et ts=4 sw=2 sts=2:
00003 
00004 #ifndef DUNE_GRID_HIERARCHICSEARCH_HH
00005 #define DUNE_GRID_HIERARCHICSEARCH_HH
00006 
00013 #include <cstddef>
00014 #include <sstream>
00015 #include <string>
00016 
00017 #include <dune/common/classname.hh>
00018 #include <dune/common/exceptions.hh>
00019 #include <dune/common/fvector.hh>
00020 
00021 #include <dune/grid/common/grid.hh>
00022 #include <dune/grid/common/gridenums.hh>
00023 
00024 namespace Dune
00025 {
00026 
00030   template<class Grid, class IS>
00031   class HierarchicSearch
00032   {
00034         enum {dim=Grid::dimension};
00035 
00037         enum {dimw=Grid::dimensionworld};
00038 
00040         typedef typename Grid::ctype ct;
00041 
00043         typedef typename Grid::template Codim<0>::Entity Entity;
00044 
00046         typedef typename Grid::template Codim<0>::EntityPointer EntityPointer;
00047 
00049         typedef typename Grid::HierarchicIterator HierarchicIterator;
00050 
00051     static std::string formatEntityInformation ( const Entity &e ) {
00052       const typename Entity::Geometry &geo = e.geometry();
00053       std::ostringstream info;
00054       info << "level=" << e.level() << " "
00055            << "partition=" << e.partitionType() << " "
00056            << "center=(" << geo.center() << ") "
00057            << "corners=[(" << geo.corner(0) << ")";
00058       for(int i = 1; i < geo.corners(); ++i)
00059         info << " (" << e.geometry().corner(i) << ")";
00060       info << "]";
00061       return info.str();
00062     }
00063 
00074     EntityPointer hFindEntity ( const Entity &e,
00075                                 const FieldVector<ct,dimw>& global) const
00076       {
00077         // loop over all child Entities
00078         const HierarchicIterator end = e.hend(e.level()+1);
00079         for( HierarchicIterator it = e.hbegin( e.level()+1 ); it != end; ++it )
00080         {
00081           FieldVector<ct,dim> local = it->geometry().local(global);
00082           if (GenericReferenceElements<double, dim>::general(it->type()).checkInside(local))
00083           {
00084             // return if we found the leaf, else search through the child entites
00085             if( is.contains( *it ) )
00086               return EntityPointer( it );
00087             else
00088               return hFindEntity( *it, global );
00089           }
00090         }
00091         std::ostringstream children;
00092         HierarchicIterator it = e.hbegin( e.level()+1 );
00093         if(it != end) {
00094           children << "{" << formatEntityInformation(*it) << "}";
00095           for( ++it; it != end; ++it )
00096             children << " {" << formatEntityInformation(*it) << "}";
00097         }
00098         DUNE_THROW(Exception, "{" << className(*this) << "} Unexpected "
00099                    "internal Error: none of the children of the entity "
00100                    "{" << formatEntityInformation(e) << "} contains "
00101                    "coordinate (" << global << ").  Children are: "
00102                    "[" << children.str() << "].");
00103       }
00104 
00105   public:
00109     HierarchicSearch(const Grid & _g, const IS & _is) : g(_g), is(_is) {};
00110 
00118     EntityPointer findEntity(const FieldVector<ct,dimw>& global) const
00119     { return findEntity<All_Partition>(global); }
00120 
00128     template<PartitionIteratorType partition>
00129     EntityPointer findEntity(const FieldVector<ct,dimw>& global) const
00130       {
00131         typedef typename Grid::template Partition<partition>::LevelGridView
00132           LevelGV;
00133         const LevelGV &gv = g.template levelView<partition>(0);
00134 
00136         typedef typename LevelGV::template Codim<0>::Iterator LevelIterator;
00137 
00138         // loop over macro level
00139         LevelIterator it = gv.template begin<0>();
00140         LevelIterator end = gv.template end<0>();
00141         for (; it != end; ++it)
00142         {
00143           const Entity &e = *it;
00144           const typename Entity::Geometry &geo = e.geometry();
00145 
00146           FieldVector< ct, dim > local = geo.local( global );
00147           if( !GenericReferenceElements< double, dim >::general( geo.type() ).checkInside( local ) )
00148             continue;
00149 
00150           if( (int(dim) != int(dimw)) && ((geo.global( local ) - global).two_norm() > 1e-8) )
00151             continue;
00152 
00153           // return if we found the leaf, else search through the child entites
00154           if( is.contains( *it ) )
00155             return EntityPointer( it );
00156           else
00157             return hFindEntity( *it, global );
00158         }
00159         DUNE_THROW( GridError, "Coordinate " << global << " is outside the grid." );
00160       }
00161     
00162   private:
00163     const Grid& g;
00164     const IS& is;
00165   };
00166   
00167 } // end namespace Dune
00168 
00169 #endif // DUNE_GRID_HIERARCHICSEARCH_HH
00170