• Main Page
  • Namespaces
  • Data Structures
  • Files
  • File List
  • Globals

/export/development/ViennaGrid/release/ViennaGrid-1.0.0/viennagrid/algorithm/refine.hpp

Go to the documentation of this file.
00001 #ifndef VIENNAGRID_ALGORITHM_REFINE_HPP
00002 #define VIENNAGRID_ALGORITHM_REFINE_HPP
00003 
00004 /* =======================================================================
00005    Copyright (c) 2011, Institute for Microelectronics,
00006                        Institute for Analysis and Scientific Computing,
00007                        TU Wien.
00008 
00009                             -----------------
00010                      ViennaGrid - The Vienna Grid Library
00011                             -----------------
00012 
00013    Authors:      Karl Rupp                           rupp@iue.tuwien.ac.at
00014                  Josef Weinbub                    weinbub@iue.tuwien.ac.at
00015                
00016    (A list of additional contributors can be found in the PDF manual)
00017 
00018    License:      MIT (X11), see file LICENSE in the base directory
00019 ======================================================================= */
00020 
00021 #include "viennagrid/forwards.h"
00022 #include "viennagrid/domain.hpp"
00023 #include "viennagrid/algorithm/centroid.hpp"
00024 #include "viennagrid/algorithm/norm.hpp"
00025 #include "viennadata/api.hpp"
00026 
00031 namespace viennagrid 
00032 {
00033   namespace detail
00034   {
00036     template <typename DomainType,
00037               typename VertexIDHandler = typename viennagrid::result_of::element_id_handler<typename DomainType::config_type,
00038                                                                                             viennagrid::point_tag>::type>
00039     struct refinement_vertex_id_requirement
00040     {
00041       typedef long        type;
00042     };
00043     
00044     template <typename DomainType>
00045     struct refinement_vertex_id_requirement< DomainType, pointer_id>
00046     {
00047       typedef typename DomainType::ERROR_VERTEX_IDs_FOR_REFINEMENT_REQUIRED   type;
00048     };
00049     
00050     
00051     
00053     template <typename ConfigTypeIn>
00054     void ensure_longest_edge_refinement(domain_t<ConfigTypeIn> const & domain_in)
00055     {
00056       typedef typename viennagrid::result_of::domain<ConfigTypeIn>::type                      DomainTypeIn;
00057       typedef typename ConfigTypeIn::cell_tag           CellTagIn;   
00058       
00059       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 0>::type                                      VertexType;
00060       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 1>::type                                      EdgeType;
00061       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, CellTagIn::dim>::type              CellType;
00062       
00063       typedef typename viennagrid::result_of::const_ncell_range<DomainTypeIn, CellTagIn::dim>::type   CellRange;          
00064       typedef typename viennagrid::result_of::iterator<CellRange>::type                                          CellIterator;         
00065       typedef typename viennagrid::result_of::const_ncell_range<CellType, 1>::type                               EdgeOnCellRange;
00066       typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type                                    EdgeOnCellIterator;            
00067       typedef typename viennagrid::result_of::const_ncell_range<EdgeType, 0>::type                               VertexOnEdgeRange;
00068       typedef typename viennagrid::result_of::iterator<VertexOnEdgeRange>::type                                  VertexOnEdgeIterator;            
00069       
00070       
00071       bool something_changed = true;
00072       //std::size_t iter_cnt = 0;
00073       
00074       while (something_changed)
00075       {
00076         something_changed = false;
00077         //std::cout << "Longest edge bisection, run " << ++iter_cnt << std::endl;
00078         
00079       
00080         CellRange cells = viennagrid::ncells<CellTagIn::dim>(domain_in);
00081         for (CellIterator cit  = cells.begin();
00082                           cit != cells.end();
00083                         ++cit)
00084         {
00085           //
00086           // check if any edges are tagged for refinement
00087           //
00088           bool has_refinement = false;
00089           EdgeOnCellRange edges_on_cell = viennagrid::ncells<1>(*cit);
00090           for (EdgeOnCellIterator eocit = edges_on_cell.begin();
00091                                   eocit != edges_on_cell.end();
00092                                 ++eocit)
00093           {
00094             if (viennadata::access<refinement_key, bool>(refinement_key())(*eocit) == true)
00095             {
00096               has_refinement = true;
00097               break;
00098             }
00099           }
00100           
00101           if (has_refinement)
00102           {
00103             //
00104             // Find longest edge
00105             //
00106             EdgeType const * longest_edge_ptr = NULL;
00107             double longest_edge_len = 0;
00108             
00109             for (EdgeOnCellIterator eocit = edges_on_cell.begin();
00110                                     eocit != edges_on_cell.end();
00111                                   ++eocit)
00112             {
00113               VertexOnEdgeRange vertices_on_edge = viennagrid::ncells<0>(*eocit);
00114               VertexOnEdgeIterator voeit = vertices_on_edge.begin();
00115               VertexType const & v0 = *voeit; ++voeit;
00116               VertexType const & v1 = *voeit;
00117               
00118               double len = viennagrid::norm(v0.point() - v1.point());
00119               if (len > longest_edge_len)
00120               {
00121                 longest_edge_len = len;
00122                 longest_edge_ptr = &(*eocit);
00123               }
00124             }
00125             
00126             //tag longest edge:
00127             if (viennadata::access<refinement_key, bool>(refinement_key())(*longest_edge_ptr) == false)
00128             {
00129               viennadata::access<refinement_key, bool>(refinement_key())(*longest_edge_ptr) = true;
00130               something_changed = true;
00131             }
00132           }
00133           
00134         } //for cells
00135 
00136 
00137       } //while
00138     } //ensure_longest_edge_refinement
00139     
00141     template <typename ConfigTypeIn>
00142     void cell_refinement_to_edge_refinement(domain_t<ConfigTypeIn> const & domain_in)
00143     {
00144       typedef typename viennagrid::result_of::domain<ConfigTypeIn>::type                      DomainTypeIn;
00145       typedef typename ConfigTypeIn::cell_tag           CellTagIn;   
00146       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, CellTagIn::dim>::type              CellType;
00147       
00148       typedef typename viennagrid::result_of::const_ncell_range<DomainTypeIn, CellTagIn::dim>::type   CellRange;          
00149       typedef typename viennagrid::result_of::iterator<CellRange>::type                                          CellIterator;         
00150       typedef typename viennagrid::result_of::const_ncell_range<CellType, 1>::type                               EdgeOnCellRange;
00151       typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type                                    EdgeOnCellIterator;            
00152       
00153       std::size_t cells_for_refinement = 0;
00154       CellRange cells = viennagrid::ncells<CellTagIn::dim>(domain_in);
00155       for (CellIterator cit  = cells.begin();
00156                         cit != cells.end();
00157                       ++cit)
00158       {
00159         if (viennadata::access<refinement_key, bool>(refinement_key())(*cit) == true)
00160         {
00161           ++cells_for_refinement;
00162           
00163           EdgeOnCellRange edges_on_cell = viennagrid::ncells<1>(*cit);
00164           for (EdgeOnCellIterator eocit = edges_on_cell.begin();
00165                                   eocit != edges_on_cell.end();
00166                                 ++eocit)
00167           {
00168             viennadata::access<refinement_key, bool>(refinement_key())(*eocit) = true;
00169           }
00170         }
00171       }
00172       
00173       //std::cout << "No. of cells tagged for refinement: " << cells_for_refinement << std::endl;
00174     }
00175     
00176     
00177 
00179     template <typename ConfigTypeIn, typename ConfigTypeOut>
00180     void refine_impl(domain_t<ConfigTypeIn> const & domain_in,
00181                     domain_t<ConfigTypeOut> & domain_out,
00182                     local_refinement_tag)
00183     {
00184       typedef domain_t<ConfigTypeIn>                      DomainTypeIn;
00185       typedef typename ConfigTypeIn::cell_tag           CellTagIn;   
00186       typedef typename ConfigTypeIn::numeric_type       NumericType;
00187 
00188       typedef typename DomainTypeIn::segment_type                                                                    SegmentTypeIn; 
00189       typedef typename viennagrid::result_of::point<ConfigTypeIn>::type                                         PointType;
00190       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 0>::type                                      VertexType;
00191       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 1>::type                                      EdgeType;
00192       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, CellTagIn::dim>::type              CellType;
00193       
00194       typedef typename viennagrid::result_of::const_ncell_range<DomainTypeIn, 0>::type                           VertexRange;          
00195       typedef typename viennagrid::result_of::iterator<VertexRange>::type                                        VertexIterator;         
00196       typedef typename viennagrid::result_of::const_ncell_range<DomainTypeIn, 1>::type                           EdgeRange;          
00197       typedef typename viennagrid::result_of::iterator<EdgeRange>::type                                          EdgeIterator;         
00198       typedef typename viennagrid::result_of::const_ncell_range<SegmentTypeIn, CellTagIn::dim>::type  CellRange;          
00199       typedef typename viennagrid::result_of::iterator<CellRange>::type                                          CellIterator;         
00200       typedef typename viennagrid::result_of::const_ncell_range<CellType, 0>::type                               VertexOnCellRange;
00201       typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type                                  VertexOnCellIterator;            
00202       typedef typename viennagrid::result_of::const_ncell_range<EdgeType, 0>::type                               VertexOnEdgeRange;
00203       typedef typename viennagrid::result_of::iterator<VertexOnEdgeRange>::type                                  VertexOnEdgeIterator;            
00204 
00205       typedef typename detail::refinement_vertex_id_requirement<domain_t<ConfigTypeIn> >::type   checked_type;
00206       
00207       //
00208       // Step 1: Write tags from cells to edges:
00209       //
00210       cell_refinement_to_edge_refinement(domain_in);
00211       ensure_longest_edge_refinement(domain_in);
00212       
00213       //
00214       // Step 2: Write old vertices to new domain
00215       //
00216       std::size_t num_vertices = 0;
00217       VertexRange vertices = viennagrid::ncells<0>(domain_in);
00218       for (VertexIterator vit  = vertices.begin();
00219                           vit != vertices.end();
00220                         ++vit)
00221       {
00222         domain_out.push_back(*vit);
00223         ++num_vertices;
00224       }
00225 
00226       //
00227       // Step 3: Each tagged edge in old domain results in a new vertex (temporarily store new vertex IDs on old domain)
00228       //
00229       EdgeRange edges = viennagrid::ncells<1>(domain_in);
00230       for (EdgeIterator eit = edges.begin();
00231                         eit != edges.end();
00232                       ++eit)
00233       {
00234         if (viennadata::access<refinement_key, bool>()(*eit) == true)
00235         {
00236           VertexType v;
00237           v.point() = viennagrid::centroid(*eit);
00238           
00239           v.id(num_vertices);
00240           domain_out.push_back(v);
00241           viennadata::access<refinement_key, std::size_t>()(*eit) = num_vertices;
00242           ++num_vertices;
00243         }
00244       }
00245       //std::cout << "Number of vertices in new domain: " << num_vertices << std::endl;
00246 
00247       //
00248       // Step 4: Now write new cells to new domain
00249       //
00250       domain_out.segments().resize(domain_in.segments().size());
00251       for (std::size_t i=0; i<domain_in.segments().size(); ++i)
00252       {
00253         //std::cout << "Working on segment " << i << std::endl;
00254         CellRange cells = viennagrid::ncells<CellTagIn::dim>(domain_in.segments()[i]);
00255         for (CellIterator cit  = cells.begin();
00256                           cit != cells.end();
00257                         ++cit)
00258         {
00259           element_refinement<CellTagIn>::apply(*cit, domain_out.segments()[i]);
00260         }
00261       }
00262 
00263       // Clean up refinement flags:
00264       for (EdgeIterator eit = edges.begin();
00265                         eit != edges.end();
00266                       ++eit)
00267       {
00268         viennadata::erase<refinement_key, std::size_t>()(*eit);
00269         viennadata::erase<refinement_key, bool>()(*eit);
00270       }
00271 
00272     }
00273     
00274     
00275     
00276     
00278     template <typename ConfigTypeIn, typename ConfigTypeOut>
00279     void refine_impl(domain_t<ConfigTypeIn> const & domain_in,
00280                     domain_t<ConfigTypeOut> & domain_out,
00281                     uniform_refinement_tag)
00282     {
00283       typedef domain_t<ConfigTypeIn>                    DomainTypeIn;
00284       typedef typename ConfigTypeIn::cell_tag           CellTagIn;   
00285       typedef typename ConfigTypeIn::numeric_type       NumericType;
00286 
00287       typedef typename DomainTypeIn::segment_type                                                                    SegmentTypeIn; 
00288       typedef typename viennagrid::result_of::point<ConfigTypeIn>::type                                         PointType;
00289       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 0>::type                                      VertexType;
00290       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, 1>::type                                      EdgeType;
00291       typedef typename viennagrid::result_of::ncell<ConfigTypeIn, CellTagIn::dim>::type              CellType;
00292       
00293       typedef typename viennagrid::result_of::const_ncell_range<DomainTypeIn, 0>::type                           VertexRange;          
00294       typedef typename viennagrid::result_of::iterator<VertexRange>::type                                        VertexIterator;         
00295       typedef typename viennagrid::result_of::const_ncell_range<DomainTypeIn, 1>::type                           EdgeRange;          
00296       typedef typename viennagrid::result_of::iterator<EdgeRange>::type                                          EdgeIterator;         
00297       typedef typename viennagrid::result_of::const_ncell_range<SegmentTypeIn, CellTagIn::dim>::type  CellRange;          
00298       typedef typename viennagrid::result_of::iterator<CellRange>::type                                          CellIterator;         
00299       typedef typename viennagrid::result_of::const_ncell_range<EdgeType, 0>::type                               VertexOnEdgeRange;
00300       typedef typename viennagrid::result_of::iterator<VertexOnEdgeRange>::type                                  VertexOnEdgeIterator;            
00301       typedef typename viennagrid::result_of::const_ncell_range<CellType, 0>::type                               VertexOnCellRange;
00302       typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type                                  VertexOnCellIterator;            
00303       
00304       typedef typename detail::refinement_vertex_id_requirement<domain_t<ConfigTypeIn> >::type   checked_type;
00305       
00306       //
00307       // Step 1: Write old vertices to new domain
00308       //
00309       std::size_t num_vertices = 0;
00310       VertexRange vertices = viennagrid::ncells<0>(domain_in);
00311       for (VertexIterator vit  = vertices.begin();
00312                           vit != vertices.end();
00313                         ++vit)
00314       {
00315         domain_out.push_back(*vit);
00316         ++num_vertices;
00317       }
00318 
00319       //
00320       // Step 2: Each edge in old domain results in a new vertex (temporarily store new vertex IDs on old domain)
00321       //
00322       EdgeRange edges = viennagrid::ncells<1>(domain_in);
00323       for (EdgeIterator eit = edges.begin();
00324                         eit != edges.end();
00325                       ++eit)
00326       {
00327         VertexType v;
00328         v.point() = viennagrid::centroid(*eit);
00329 
00330         v.id(num_vertices);
00331         domain_out.push_back(v);
00332         viennadata::access<refinement_key, std::size_t>()(*eit) = num_vertices;
00333         viennadata::access<refinement_key, bool>()(*eit) = true;
00334         ++num_vertices;
00335       }
00336       //std::cout << "Number of vertices in new domain: " << num_vertices << std::endl;
00337 
00338       //
00339       // Step 3: Now write new cells to new domain
00340       //
00341       domain_out.segments().resize(domain_in.segments().size());
00342       for (std::size_t i=0; i<domain_in.segments().size(); ++i)
00343       {
00344         //std::cout << "Working on segment " << i << std::endl;
00345         std::size_t counter = 0;
00346         CellRange cells = viennagrid::ncells<CellTagIn::dim>(domain_in.segments()[i]);
00347         for (CellIterator cit  = cells.begin();
00348                           cit != cells.end();
00349                         ++cit)
00350         {
00351           element_refinement<CellTagIn>::apply(*cit, domain_out.segments()[i]);
00352           ++counter;
00353         }
00354         //std::cout << "Refined cells in segment: " << counter << std::endl;
00355       }
00356 
00357       // Clean up refinement flags:
00358       for (EdgeIterator eit = edges.begin();
00359                         eit != edges.end();
00360                       ++eit)
00361       {
00362         viennadata::erase<refinement_key, std::size_t>()(*eit);
00363         viennadata::erase<refinement_key, bool>()(*eit);
00364       }
00365     }
00366     
00367   } //namespace detail
00368   
00370   template <typename DomainType, typename RefinementTag>
00371   class refinement_proxy
00372   {
00373     public:
00374       refinement_proxy(DomainType const & domain_in,
00375                        RefinementTag const & refinement_tag) : domain_in_(domain_in), tag_(refinement_tag) {}
00376       
00377       DomainType const & get() const 
00378       {
00379         return domain_in_;
00380       }
00381       
00382       RefinementTag const & tag() const { return tag_; }
00383       
00384     private:
00385       DomainType const & domain_in_;
00386       RefinementTag const & tag_;
00387   };
00388   
00389   
00391   template <typename ConfigTypeIn, typename RefinementTag>
00392   refinement_proxy< domain_t<ConfigTypeIn>, RefinementTag >
00393   refine(domain_t<ConfigTypeIn> const & domain_in, RefinementTag const & tag)
00394   {
00395     typedef typename detail::refinement_vertex_id_requirement<domain_t<ConfigTypeIn> >::type   checked_type;
00396     return refinement_proxy< domain_t<ConfigTypeIn>, RefinementTag >(domain_in, tag);
00397   }
00398   
00400   template <typename ConfigTypeIn>
00401   refinement_proxy< domain_t<ConfigTypeIn>, uniform_refinement_tag >
00402   refine_uniformly(domain_t<ConfigTypeIn> const & domain_in)
00403   {
00404     typedef typename detail::refinement_vertex_id_requirement<domain_t<ConfigTypeIn> >::type   checked_type;
00405     return refinement_proxy< domain_t<ConfigTypeIn>, uniform_refinement_tag >(domain_in, uniform_refinement_tag());
00406   }
00407   
00409   template <typename ConfigTypeIn>
00410   refinement_proxy< domain_t<ConfigTypeIn>, local_refinement_tag >
00411   refine_locally(domain_t<ConfigTypeIn> const & domain_in)
00412   {
00413     typedef typename detail::refinement_vertex_id_requirement<domain_t<ConfigTypeIn> >::type   checked_type;
00414     return refinement_proxy< domain_t<ConfigTypeIn>, local_refinement_tag >(domain_in, local_refinement_tag());
00415   }
00416   
00417 } 
00418 
00419 #endif

Generated on Wed Sep 14 2011 19:21:30 for ViennaGrid by  doxygen 1.7.1