00001 #ifndef VIENNAGRID_ALGORITHM_REFINE_HPP
00002 #define VIENNAGRID_ALGORITHM_REFINE_HPP
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
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
00073
00074 while (something_changed)
00075 {
00076 something_changed = false;
00077
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
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
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
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 }
00135
00136
00137 }
00138 }
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
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
00209
00210 cell_refinement_to_edge_refinement(domain_in);
00211 ensure_longest_edge_refinement(domain_in);
00212
00213
00214
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
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
00246
00247
00248
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
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
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
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
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
00337
00338
00339
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
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
00355 }
00356
00357
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 }
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