ViennaGrid - The Vienna Grid Library: viennagrid/algorithm/detail/refine

Go to the documentation of this file.
00001 #ifndef VIENNAGRID_ALGORITHM_DETAIL_REFINE_TRI_HPP
00002 #define VIENNAGRID_ALGORITHM_DETAIL_REFINE_TRI_HPP
00003 
00004 /* =======================================================================
00005    Copyright (c) 2011-2014, Institute for Microelectronics,
00006                             Institute for Analysis and Scientific Computing,
00007                             TU Wien.
00008 
00009                             -----------------
00010                      ViennaGrid - The Vienna Grid Library
00011                             -----------------
00012 
00013    License:      MIT (X11), see file LICENSE in the base directory
00014 ======================================================================= */
00015 
00016 #include "viennagrid/forwards.hpp"
00017 #include "viennagrid/mesh/segmentation.hpp"
00018 #include "viennagrid/topology/vertex.hpp"
00019 #include "viennagrid/topology/line.hpp"
00020 #include "viennagrid/algorithm/norm.hpp"
00021 
00026 namespace viennagrid
00027 {
00028   namespace detail
00029   {
00030 
00031     template<typename ElementsVerticesHandleContainerT, typename VertexHandleContainer>
00032     void add_refinement_element(ElementsVerticesHandleContainerT & elements_vertices,
00033                                 VertexHandleContainer vertex_handle_container,
00034                                 unsigned int i0, unsigned int i1, unsigned int i2)
00035     {
00036       elements_vertices.resize( elements_vertices.size()+1 );
00037 
00038       elements_vertices.back().resize(3);
00039       elements_vertices.back()[0] = *viennagrid::advance(vertex_handle_container.begin(), i0);
00040       elements_vertices.back()[1] = *viennagrid::advance(vertex_handle_container.begin(), i1);
00041       elements_vertices.back()[2] = *viennagrid::advance(vertex_handle_container.begin(), i2);
00042     }
00043 
00044 
00046     template<>
00047     struct element_refinement<triangle_tag>
00048     {
00049 
00051       template<typename ElementType, typename MeshT,
00052                typename ElementsVerticesHandleContainerT,
00053                typename VertexCopyMapT,
00054                typename EdgeRefinementFlagAccessor, typename EdgeToVertexHandleAccessor>
00055       static void apply0(ElementType const & element_in, MeshT const &,
00056                          ElementsVerticesHandleContainerT & elements_vertices,
00057                          VertexCopyMapT & vertex_copy_map_,
00058                          EdgeRefinementFlagAccessor const &, EdgeToVertexHandleAccessor const &)
00059       {
00060         typedef typename viennagrid::result_of::const_element_range<ElementType, viennagrid::vertex_tag>::type            VertexOnCellRange;
00061         typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type         VertexOnCellIterator;
00062 
00063         typedef typename viennagrid::result_of::vertex_handle<MeshT>::type             VertexHandleType;
00064 
00065         static_array<VertexHandleType, boundary_elements<triangle_tag, vertex_tag>::num> vertex_handles;
00066 
00067         //
00068         // Step 1: Get vertices on the new mesh
00069         //
00070 
00071         //grab existing vertices:
00072         VertexOnCellRange vertices_on_cell = viennagrid::elements<viennagrid::vertex_tag>(element_in);
00073         VertexOnCellIterator vocit = vertices_on_cell.begin();
00074         vertex_handles[0] = vertex_copy_map_(*vocit); ++vocit;
00075         vertex_handles[1] = vertex_copy_map_(*vocit); ++vocit;
00076         vertex_handles[2] = vertex_copy_map_(*vocit);
00077 
00078         //
00079         // Step 2: Add new cells to new mesh:
00080         //
00081 
00082         add_refinement_element( elements_vertices, vertex_handles, 0, 1, 2);
00083 
00084       } //apply0()
00085 
00086 
00096       template<typename ElementType, typename MeshT,
00097                typename ElementsVerticesHandleContainerT,
00098                typename VertexCopyMapT,
00099                typename EdgeRefinementFlagAccessor, typename EdgeToVertexHandleAccessor>
00100       static void apply1(ElementType const & element_in, MeshT const &,
00101                          ElementsVerticesHandleContainerT & elements_vertices,
00102                          VertexCopyMapT & vertex_copy_map_,
00103                          EdgeRefinementFlagAccessor const & edge_refinement_flag_accessor, EdgeToVertexHandleAccessor const & edge_to_vertex_handle_accessor)
00104       {
00105         typedef typename viennagrid::result_of::const_element_range<ElementType, viennagrid::vertex_tag>::type            VertexOnCellRange;
00106         typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type         VertexOnCellIterator;
00107         typedef typename viennagrid::result_of::const_element_range<ElementType, viennagrid::line_tag>::type            EdgeOnCellRange;
00108         typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type           EdgeOnCellIterator;
00109 
00110         typedef typename viennagrid::result_of::handle<MeshT, viennagrid::vertex_tag>::type             VertexHandleType;
00111         typedef typename viennagrid::result_of::element<MeshT, viennagrid::line_tag>::type             EdgeType;
00112 
00113         const unsigned int num_vertices = boundary_elements<triangle_tag, vertex_tag>::num;
00114         static_array<VertexHandleType, num_vertices+1> vertex_handles;
00115 
00116 
00117         //
00118         // Step 1: Get vertices on the new mesh
00119         //
00120 
00121         //grab existing vertices:
00122         VertexOnCellRange vertices_on_cell(element_in);
00123         VertexOnCellIterator vocit = vertices_on_cell.begin();
00124         vertex_handles[0] = vertex_copy_map_(*vocit); ++vocit;
00125         vertex_handles[1] = vertex_copy_map_(*vocit); ++vocit;
00126         vertex_handles[2] = vertex_copy_map_(*vocit);
00127 
00128         //add vertices from edge
00129         EdgeOnCellRange edges_on_cell(element_in);
00130         std::size_t offset = 0;
00131         EdgeOnCellIterator eocit = edges_on_cell.begin();
00132         EdgeType const & e0 = *eocit; ++eocit;
00133         EdgeType const & e1 = *eocit; ++eocit;
00134         EdgeType const & e2 = *eocit;
00135 
00136         if ( edge_refinement_flag_accessor(e0) )
00137         {
00138           vertex_handles[3] = edge_to_vertex_handle_accessor(e0);
00139           offset = 0;
00140         }
00141         else if ( edge_refinement_flag_accessor(e1) )
00142         {
00143           vertex_handles[3] = edge_to_vertex_handle_accessor(e1);
00144           offset = 2;
00145         }
00146         else if ( edge_refinement_flag_accessor(e2) )
00147         {
00148           vertex_handles[3] = edge_to_vertex_handle_accessor(e2);
00149           offset = 1;
00150         }
00151         else
00152         {
00153           assert( (2 < 2) && "Logic error: Triangle does not have an edges for bisection!");
00154         }
00155 
00156         //
00157         // Step 2: Add new cells to new mesh:
00158         //
00159         //0-3-2
00160         add_refinement_element( elements_vertices, vertex_handles, (offset + 0) % num_vertices, 3, (offset + 2) % num_vertices );
00161 
00162         //3-1-2:
00163         add_refinement_element( elements_vertices, vertex_handles, 3, (offset + 1) % num_vertices, (offset + 2) % num_vertices );
00164       } //apply1()
00165 
00166 
00177       template<typename ElementType, typename MeshT,
00178                typename ElementsVerticesHandleContainerT,
00179                typename VertexCopyMapT,
00180                typename EdgeRefinementFlagAccessor, typename EdgeToVertexHandleAccessor>
00181       static void apply2(ElementType const & element_in, MeshT const & mesh,
00182                          ElementsVerticesHandleContainerT & elements_vertices,
00183                          VertexCopyMapT & vertex_copy_map_,
00184                          EdgeRefinementFlagAccessor const & edge_refinement_flag_accessor, EdgeToVertexHandleAccessor const & edge_to_vertex_handle_accessor)
00185       {
00186         typedef typename viennagrid::result_of::const_element_range<ElementType, viennagrid::vertex_tag>::type            VertexOnCellRange;
00187         typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type         VertexOnCellIterator;
00188         typedef typename viennagrid::result_of::const_element_range<ElementType, viennagrid::line_tag>::type            EdgeOnCellRange;
00189         typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type           EdgeOnCellIterator;
00190 
00191         typedef typename viennagrid::result_of::handle<MeshT, viennagrid::vertex_tag>::type             VertexHandleType;
00192         typedef typename viennagrid::result_of::element<MeshT, viennagrid::line_tag>::type             EdgeType;
00193 
00194         const unsigned int num_vertices = boundary_elements<triangle_tag, vertex_tag>::num;
00195         static_array<VertexHandleType, num_vertices+2> vertex_handles;
00196 
00197 
00198         //
00199         // Step 1: Get vertices on the new mesh
00200         //
00201 
00202         //grab existing vertices:
00203         VertexOnCellRange vertices_on_cell = viennagrid::elements<viennagrid::vertex_tag>(element_in);
00204         VertexOnCellIterator vocit = vertices_on_cell.begin();
00205         vertex_handles[0] = vertex_copy_map_(*vocit); ++vocit;
00206         vertex_handles[1] = vertex_copy_map_(*vocit); ++vocit;
00207         vertex_handles[2] = vertex_copy_map_(*vocit);
00208 
00209         //Find rotation offset such that first two edges are to be refined
00210         EdgeOnCellRange edges_on_cell = viennagrid::elements<viennagrid::line_tag>(element_in);
00211         std::size_t offset = 0;
00212 
00213         EdgeOnCellIterator eocit = edges_on_cell.begin();
00214         EdgeType const & e0 = *eocit; ++eocit;
00215         EdgeType const & e1 = *eocit; ++eocit;
00216         EdgeType const & e2 = *eocit;
00217 
00218         if ( edge_refinement_flag_accessor(e0) && edge_refinement_flag_accessor(e1) )
00219         {
00220           vertex_handles[3] = edge_to_vertex_handle_accessor(e1);
00221           vertex_handles[4] = edge_to_vertex_handle_accessor(e0);
00222           offset = 2;
00223         }
00224         else if ( edge_refinement_flag_accessor(e0) && edge_refinement_flag_accessor(e2) )
00225         {
00226           vertex_handles[3] = edge_to_vertex_handle_accessor(e0);
00227           vertex_handles[4] = edge_to_vertex_handle_accessor(e2);
00228           offset = 0;
00229         }
00230         else if ( edge_refinement_flag_accessor(e1) && edge_refinement_flag_accessor(e2) )
00231         {
00232           vertex_handles[3] = edge_to_vertex_handle_accessor(e2);
00233           vertex_handles[4] = edge_to_vertex_handle_accessor(e1);
00234           offset = 1;
00235         }
00236         else
00237         {
00238           assert( (2 < 2) && "Logic error: Triangle does not have two edges for bisection!");
00239         }
00240 
00241         //
00242         // Step 2: Add new cells to new mesh:
00243         //
00244 
00245         //3-1-4:
00246         add_refinement_element( elements_vertices, vertex_handles, 3, (offset + 1) % num_vertices, 4 );
00247 
00248         //split second-longest edge
00249         VertexHandleType vh0 = vertex_handles[(offset + 0) % num_vertices];
00250         VertexHandleType vh1 = vertex_handles[(offset + 1) % num_vertices];
00251         VertexHandleType vh2 = vertex_handles[(offset + 2) % num_vertices];
00252         double len_edge1 = viennagrid::norm( viennagrid::point(mesh, vh1) - viennagrid::point(mesh, vh0) );
00253         double len_edge2 = viennagrid::norm( viennagrid::point(mesh, vh2) - viennagrid::point(mesh, vh1) );
00254 
00255         if (len_edge1 > len_edge2) //split edge [v0, v1] again
00256         {
00257           //0-3-2:
00258           add_refinement_element( elements_vertices, vertex_handles, (offset + 0) % num_vertices, 3, (offset + 2) % num_vertices );
00259 
00260           //2-3-4:
00261           add_refinement_element( elements_vertices, vertex_handles, (offset + 2) % num_vertices, 3, 4 );
00262         }
00263         else //split edge [v1, v2]
00264         {
00265           //0-3-4:
00266           add_refinement_element( elements_vertices, vertex_handles, (offset + 0) % num_vertices, 3, 4 );
00267 
00268           //0-4-2:
00269           add_refinement_element( elements_vertices, vertex_handles, (offset + 0) % num_vertices, 4, (offset + 2) % num_vertices );
00270         }
00271 
00272 
00273       } //apply2()
00274 
00275 
00276 
00277 
00279       template<typename ElementType, typename MeshT,
00280                typename ElementsVerticesHandleContainerT,
00281                typename VertexCopyMapT,
00282                typename EdgeRefinementFlagAccessor, typename EdgeToVertexHandleAccessor>
00283       static void apply3(ElementType const & element_in, MeshT const &,
00284                          ElementsVerticesHandleContainerT & elements_vertices,
00285                          VertexCopyMapT & vertex_copy_map_,
00286                          EdgeRefinementFlagAccessor const &, EdgeToVertexHandleAccessor const & edge_to_vertex_handle_accessor)
00287       {
00288         typedef typename viennagrid::result_of::const_element_range<ElementType, viennagrid::vertex_tag>::type            VertexOnCellRange;
00289         typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type         VertexOnCellIterator;
00290         typedef typename viennagrid::result_of::const_element_range<ElementType, viennagrid::line_tag>::type            EdgeOnCellRange;
00291         typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type           EdgeOnCellIterator;
00292 
00293         typedef typename viennagrid::result_of::handle<MeshT, viennagrid::vertex_tag>::type             VertexHandleType;
00294 
00295         const unsigned int num_vertices = boundary_elements<triangle_tag, vertex_tag>::num;
00296         const unsigned int num_lines = boundary_elements<triangle_tag, line_tag>::num;
00297 
00298         static_array<VertexHandleType, num_vertices+num_lines> vertex_handles;
00299 
00300         //
00301         // Step 1: Get vertices on the new mesh
00302         //
00303 
00304         //grab existing vertices:
00305         VertexOnCellRange vertices_on_cell = viennagrid::elements<viennagrid::vertex_tag>(element_in);
00306         VertexOnCellIterator vocit = vertices_on_cell.begin();
00307         vertex_handles[0] = vertex_copy_map_(*vocit); ++vocit;
00308         vertex_handles[1] = vertex_copy_map_(*vocit); ++vocit;
00309         vertex_handles[2] = vertex_copy_map_(*vocit);
00310 
00311         //add vertices from edge
00312         EdgeOnCellRange edges_on_cell = viennagrid::elements<viennagrid::line_tag>(element_in);
00313         EdgeOnCellIterator eocit = edges_on_cell.begin();
00314         vertex_handles[3] = edge_to_vertex_handle_accessor(*eocit); ++eocit;
00315         vertex_handles[4] = edge_to_vertex_handle_accessor(*eocit); ++eocit;
00316         vertex_handles[5] = edge_to_vertex_handle_accessor(*eocit);
00317 
00318         //
00319         // Step 2: Add new cells to new mesh:
00320         //
00321 
00322         //0-3-4:
00323         add_refinement_element( elements_vertices, vertex_handles, 0, 3, 4 );
00324 
00325         //3-1-5:
00326         add_refinement_element( elements_vertices, vertex_handles, 3, 1, 5 );
00327 
00328         //4-5-2:
00329         add_refinement_element( elements_vertices, vertex_handles, 4, 5, 2 );
00330 
00331         //4-3-5:
00332         add_refinement_element( elements_vertices, vertex_handles, 4, 3, 5 );
00333 
00334       } //apply3()
00335 
00336 
00346       template<typename ElementT, typename MeshT,
00347                typename ElementsVerticesHandleContainerT, typename VertexCopyMapT,
00348                typename EdgeRefinementFlagAccessor, typename EdgeToVertexHandleAccessor>
00349       static void apply(ElementT const & element_in, MeshT const & mesh,
00350                         ElementsVerticesHandleContainerT & elements_vertices,
00351                         VertexCopyMapT & vertex_copy_map_,
00352                         EdgeRefinementFlagAccessor const & edge_refinement_flag_accessor,
00353                         EdgeToVertexHandleAccessor const & edge_to_vertex_handle_accessor)
00354       {
00355         typedef typename viennagrid::result_of::const_element_range<ElementT, viennagrid::line_tag>::type            EdgeOnCellRange;
00356         typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type                 EdgeOnCellIterator;
00357 
00358         std::size_t edges_to_refine = 0;
00359         EdgeOnCellRange edges_on_cell(element_in);
00360         for (EdgeOnCellIterator eocit = edges_on_cell.begin();
00361                                 eocit != edges_on_cell.end();
00362                               ++eocit)
00363         {
00364           if ( edge_refinement_flag_accessor(*eocit) )
00365             ++edges_to_refine;
00366         }
00367 
00368         switch (edges_to_refine)
00369         {
00370           case 0: apply0(element_in, mesh, elements_vertices, vertex_copy_map_, edge_refinement_flag_accessor,  edge_to_vertex_handle_accessor); break;
00371           case 1: apply1(element_in, mesh, elements_vertices, vertex_copy_map_, edge_refinement_flag_accessor, edge_to_vertex_handle_accessor); break;
00372           case 2: apply2(element_in, mesh, elements_vertices, vertex_copy_map_, edge_refinement_flag_accessor, edge_to_vertex_handle_accessor); break;
00373           case 3: apply3(element_in, mesh, elements_vertices, vertex_copy_map_, edge_refinement_flag_accessor, edge_to_vertex_handle_accessor); break;
00374           default: //nothing to do...
00375                   break;
00376         }
00377       } //apply()
00378     };
00379 
00380   } // namespace detail
00381 }
00382 
00383 #endif
00384