00001 #ifndef VIENNAGRID_ALGORITHM_VORONOI_HPP
00002 #define VIENNAGRID_ALGORITHM_VORONOI_HPP
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include <iostream>
00022 #include <sstream>
00023 #include <string>
00024 #include <stdexcept>
00025
00026 #include "viennagrid/forwards.h"
00027 #include "viennagrid/algorithm/circumcenter.hpp"
00028 #include "viennagrid/algorithm/spanned_volume.hpp"
00029 #include "viennagrid/algorithm/volume.hpp"
00030
00031 #include "viennadata/api.hpp"
00032
00038 namespace viennagrid
00039 {
00040 namespace detail
00041 {
00043 template <typename DomainType,
00044 typename InterfaceAreaKey,
00045 typename BoxVolumeKey>
00046 void write_voronoi_info(DomainType const & domain,
00047 InterfaceAreaKey const & interface_key,
00048 BoxVolumeKey const & box_volume_key,
00049 viennagrid::simplex_tag<1>)
00050 {
00051 typedef typename DomainType::config_type Config;
00052 typedef typename Config::cell_tag CellTag;
00053
00054 typedef typename viennagrid::result_of::point<Config>::type PointType;
00055 typedef typename viennagrid::result_of::ncell<Config, 0>::type VertexType;
00056 typedef typename viennagrid::result_of::ncell<Config, CellTag::dim>::type CellType;
00057
00058 typedef typename viennagrid::result_of::const_ncell_range<DomainType, CellTag::dim>::type CellRange;
00059 typedef typename viennagrid::result_of::iterator<CellRange>::type CellIterator;
00060
00061
00062 typedef typename viennagrid::result_of::const_ncell_range<CellType, 0>::type VertexOnCellRange;
00063 typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type VertexOnCellIterator;
00064
00065 CellRange cells = viennagrid::ncells<CellTag::dim>(domain);
00066 for (CellIterator cit = cells.begin();
00067 cit != cells.end();
00068 ++cit)
00069 {
00070 viennadata::access<InterfaceAreaKey, double>(interface_key)(*cit) = 1;
00071 viennadata::access<BoxVolumeKey, double>(box_volume_key)(*cit) += volume(*cit) / 2.0;
00072
00073 VertexOnCellRange vertices_on_cell = viennagrid::ncells<0>(*cit);
00074 for (VertexOnCellIterator vocit = vertices_on_cell.begin();
00075 vocit != vertices_on_cell.end();
00076 ++vocit)
00077 {
00078 viennadata::access<BoxVolumeKey, double>(box_volume_key)(*vocit) += volume(*cit) / 2.0;
00079 }
00080 }
00081 }
00082
00084 template <typename DomainType,
00085 typename InterfaceAreaKey,
00086 typename BoxVolumeKey>
00087 void write_voronoi_info(DomainType const & domain,
00088 InterfaceAreaKey const & interface_key,
00089 BoxVolumeKey const & box_volume_key,
00090 viennagrid::hypercube_tag<1>)
00091 {
00092 write_voronoi_info(domain, interface_key, box_volume_key, viennagrid::simplex_tag<1>());
00093 }
00094
00095
00096
00097
00099 template <typename DomainType,
00100 typename InterfaceAreaKey,
00101 typename BoxVolumeKey>
00102 void write_voronoi_info(DomainType const & domain,
00103 InterfaceAreaKey const & interface_key,
00104 BoxVolumeKey const & box_volume_key,
00105 viennagrid::quadrilateral_tag)
00106 {
00107
00108 typedef typename DomainType::config_type Config;
00109 typedef typename Config::cell_tag CellTag;
00110 typedef typename viennagrid::result_of::point<Config>::type PointType;
00111 typedef typename viennagrid::result_of::ncell<Config, 0>::type VertexType;
00112 typedef typename viennagrid::result_of::ncell<Config, 1>::type EdgeType;
00113 typedef typename viennagrid::result_of::ncell<Config, CellTag::dim>::type CellType;
00114
00115 typedef typename viennagrid::result_of::const_ncell_range<DomainType, CellTag::dim>::type CellRange;
00116 typedef typename viennagrid::result_of::iterator<CellRange>::type CellIterator;
00117
00118 typedef typename viennagrid::result_of::const_ncell_range<CellType, 1>::type EdgeOnCellRange;
00119 typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type EdgeOnCellIterator;
00120
00121 typedef typename viennagrid::result_of::const_ncell_range<EdgeType, 0>::type VertexOnEdgeRange;
00122 typedef typename viennagrid::result_of::iterator<VertexOnEdgeRange>::type VertexOnEdgeIterator;
00123
00124
00125
00126
00127
00128 CellRange cells = viennagrid::ncells<CellTag::dim>(domain);
00129 for (CellIterator cit = cells.begin();
00130 cit != cells.end();
00131 ++cit)
00132 {
00133 PointType circ_center = circumcenter(*cit);
00134
00135
00136 EdgeOnCellRange edges_on_cell = viennagrid::ncells<1>(*cit);
00137 for (EdgeOnCellIterator eocit = edges_on_cell.begin();
00138 eocit != edges_on_cell.end();
00139 ++eocit)
00140 {
00141 PointType edge_midpoint = circumcenter(*eocit);
00142
00143
00144 viennadata::access<InterfaceAreaKey, double>(interface_key)(*eocit) += spanned_volume(circ_center, edge_midpoint);
00145
00146
00147 double edge_contribution = 0;
00148 VertexOnEdgeRange vertices_on_edge = viennagrid::ncells<0>(*eocit);
00149 for (VertexOnEdgeIterator voeit = vertices_on_edge.begin();
00150 voeit != vertices_on_edge.end();
00151 ++voeit)
00152 {
00153 double contribution = spanned_volume(circ_center, edge_midpoint, voeit->point());
00154 edge_contribution += contribution;
00155 viennadata::access<BoxVolumeKey, double>(box_volume_key)(*voeit) += contribution;
00156 }
00157 viennadata::access<BoxVolumeKey, double>(box_volume_key)(*eocit) += edge_contribution;
00158 }
00159
00160 }
00161
00162 }
00163
00164
00166 template <typename DomainType,
00167 typename InterfaceAreaKey,
00168 typename BoxVolumeKey>
00169 void write_voronoi_info(DomainType const & domain,
00170 InterfaceAreaKey const & interface_key,
00171 BoxVolumeKey const & box_volume_key,
00172 viennagrid::triangle_tag)
00173 {
00174 typedef typename DomainType::config_type Config;
00175 typedef typename Config::cell_tag CellTag;
00176 typedef typename viennagrid::result_of::point<Config>::type PointType;
00177 typedef typename viennagrid::result_of::ncell<Config, 0>::type VertexType;
00178 typedef typename viennagrid::result_of::ncell<Config, 1>::type EdgeType;
00179 typedef typename viennagrid::result_of::ncell<Config, CellTag::dim>::type CellType;
00180
00181 typedef typename viennagrid::result_of::const_ncell_range<DomainType, CellTag::dim>::type CellRange;
00182 typedef typename viennagrid::result_of::iterator<CellRange>::type CellIterator;
00183
00184 typedef typename viennagrid::result_of::const_ncell_range<DomainType, 1>::type EdgeRange;
00185 typedef typename viennagrid::result_of::iterator<EdgeRange>::type EdgeIterator;
00186
00187 typedef typename viennagrid::result_of::const_ncell_range<CellType, 1>::type EdgeOnCellRange;
00188 typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type EdgeOnCellIterator;
00189
00190 typedef typename viennagrid::result_of::const_ncell_range<EdgeType, 0>::type VertexOnEdgeRange;
00191 typedef typename viennagrid::result_of::iterator<VertexOnEdgeRange>::type VertexOnEdgeIterator;
00192
00193
00194
00195
00196
00197 CellRange cells = viennagrid::ncells<CellTag::dim>(domain);
00198 for (CellIterator cit = cells.begin();
00199 cit != cells.end();
00200 ++cit)
00201 {
00202 PointType circ_center = circumcenter(*cit);
00203
00204
00205 EdgeOnCellRange edges_on_cell = viennagrid::ncells<1>(*cit);
00206 for (EdgeOnCellIterator eocit = edges_on_cell.begin();
00207 eocit != edges_on_cell.end();
00208 ++eocit)
00209 {
00210 viennadata::access<InterfaceAreaKey, std::vector<PointType> >(interface_key)(*eocit).push_back(circ_center);
00211 }
00212 }
00213
00214
00215
00216
00217
00218 EdgeRange edges = viennagrid::ncells<1>(domain);
00219 for (EdgeIterator eit = edges.begin();
00220 eit != edges.end();
00221 ++eit)
00222 {
00223 double edge_length = volume(*eit);
00224
00225 {
00226 std::vector<PointType> & circ_centers = viennadata::access<InterfaceAreaKey, std::vector<PointType> >(interface_key)(*eit);
00227
00228 double interface_area = -42.0;
00229 if (circ_centers.size() == 1)
00230 interface_area = spanned_volume(circ_centers[0], circumcenter(*eit));
00231 else if (circ_centers.size() == 2)
00232 interface_area = spanned_volume(circ_centers[0], circ_centers[1]);
00233 else
00234 throw "More than two circum centers for an edge in two dimensions!";
00235
00236 viennadata::access<InterfaceAreaKey, double>(interface_key)(*eit) = interface_area;
00237 double volume_contribution = interface_area * edge_length / 4.0;
00238 viennadata::access<BoxVolumeKey, double>(box_volume_key)(*eit) = 2.0 * volume_contribution;
00239 viennadata::access<BoxVolumeKey, double>(box_volume_key)(viennagrid::ncells<0>(*eit)[0]) += volume_contribution;
00240 viennadata::access<BoxVolumeKey, double>(box_volume_key)(viennagrid::ncells<0>(*eit)[1]) += volume_contribution;
00241 }
00242
00243
00244 viennadata::erase<InterfaceAreaKey, std::vector<PointType> >()(*eit);
00245 }
00246
00247 }
00248
00249
00250
00251
00252
00253
00255 template <typename DomainType,
00256 typename InterfaceAreaKey,
00257 typename BoxVolumeKey>
00258 void write_voronoi_info(DomainType const & domain,
00259 InterfaceAreaKey const & interface_key,
00260 BoxVolumeKey const & box_volume_key,
00261 viennagrid::tetrahedron_tag)
00262 {
00263 typedef typename DomainType::config_type Config;
00264 typedef typename Config::cell_tag CellTag;
00265 typedef typename viennagrid::result_of::point<Config>::type PointType;
00266 typedef typename viennagrid::result_of::ncell<Config, 0>::type VertexType;
00267 typedef typename viennagrid::result_of::ncell<Config, 1>::type EdgeType;
00268 typedef typename viennagrid::result_of::ncell<Config, 2>::type FacetType;
00269 typedef typename viennagrid::result_of::ncell<Config, CellTag::dim>::type CellType;
00270
00271 typedef typename viennagrid::result_of::const_ncell_range<DomainType, CellTag::dim>::type CellRange;
00272 typedef typename viennagrid::result_of::iterator<CellRange>::type CellIterator;
00273
00274 typedef typename viennagrid::result_of::const_ncell_range<DomainType, CellTag::dim-1>::type FacetRange;
00275 typedef typename viennagrid::result_of::iterator<FacetRange>::type FacetIterator;
00276
00277 typedef typename viennagrid::result_of::const_ncell_range<DomainType, 1>::type EdgeRange;
00278 typedef typename viennagrid::result_of::iterator<EdgeRange>::type EdgeIterator;
00279
00280 typedef typename viennagrid::result_of::const_ncell_range<CellType, CellTag::dim-1>::type FacetOnCellRange;
00281 typedef typename viennagrid::result_of::iterator<FacetOnCellRange>::type FacetOnCellIterator;
00282
00283 typedef typename viennagrid::result_of::const_ncell_range<FacetType, 1>::type EdgeOnFacetRange;
00284 typedef typename viennagrid::result_of::iterator<EdgeOnFacetRange>::type EdgeOnFacetIterator;
00285
00286 typedef typename viennagrid::result_of::const_ncell_range<EdgeType, 0>::type VertexOnEdgeRange;
00287 typedef typename viennagrid::result_of::iterator<VertexOnEdgeRange>::type VertexOnEdgeIterator;
00288
00289
00290
00291
00292
00293 CellRange cells = viennagrid::ncells<CellTag::dim>(domain);
00294 for (CellIterator cit = cells.begin();
00295 cit != cells.end();
00296 ++cit)
00297 {
00298 PointType circ_center = circumcenter(*cit);
00299
00300 FacetOnCellRange facets_on_cell = viennagrid::ncells<CellTag::dim-1>(*cit);
00301 for (FacetOnCellIterator focit = facets_on_cell.begin();
00302 focit != facets_on_cell.end();
00303 ++focit)
00304 {
00305 viennadata::access<InterfaceAreaKey, std::vector<PointType> >(interface_key)(*focit).push_back(circ_center);
00306 }
00307 }
00308
00309
00310
00311
00312
00313 FacetRange facets = viennagrid::ncells<CellTag::dim-1>(domain);
00314 for (FacetIterator fit = facets.begin();
00315 fit != facets.end();
00316 ++fit)
00317 {
00318 std::vector<PointType> & circ_centers = viennadata::access<InterfaceAreaKey, std::vector<PointType> >(interface_key)(*fit);
00319
00320 EdgeOnFacetRange edges_on_facet = viennagrid::ncells<1>(*fit);
00321 for (EdgeOnFacetIterator eofit = edges_on_facet.begin();
00322 eofit != edges_on_facet.end();
00323 ++eofit)
00324 {
00325 viennadata::access<InterfaceAreaKey, std::vector<PointType> >(interface_key)(*eofit).push_back(circ_centers[0]);
00326
00327 if (circ_centers.size() == 1)
00328 {
00329 viennadata::access<InterfaceAreaKey,
00330 std::vector<std::pair<PointType, PointType> >
00331 >(interface_key)(*eofit).push_back(std::make_pair(circ_centers[0], circumcenter(*fit)));
00332 viennadata::access<InterfaceAreaKey,
00333 std::vector<std::pair<PointType, PointType> >
00334 >(interface_key)(*eofit).push_back(std::make_pair(circumcenter(*eofit), circumcenter(*fit)));
00335 }
00336 else if (circ_centers.size() == 2)
00337 viennadata::access<InterfaceAreaKey,
00338 std::vector<std::pair<PointType, PointType> >
00339 >(interface_key)(*eofit).push_back(std::make_pair(circ_centers[0], circ_centers[1]));
00340 else
00341 {
00342 std::cerr << "circ_centers.size() = " << circ_centers.size() << std::endl;
00343 std::cerr << "*fit: " << *fit << std::endl;
00344 throw "More than two circumcenters for a facet in three dimensions!";
00345 }
00346
00347 }
00348 }
00349
00350
00351
00352
00353
00354 EdgeRange edges = viennagrid::ncells<1>(domain);
00355 for (EdgeIterator eit = edges.begin();
00356 eit != edges.end();
00357 ++eit)
00358 {
00359
00360 VertexOnEdgeRange vertices_on_edge = viennagrid::ncells<0>(*eit);
00361 VertexOnEdgeIterator voeit = vertices_on_edge.begin();
00362
00363 VertexType const & v0 = *voeit;
00364 ++voeit;
00365 VertexType const & v1 = *voeit;
00366
00367 PointType const & p0 = v0.point();
00368 PointType const & p1 = v1.point();
00369
00370
00371 double edge_length = spanned_volume(p0, p1);
00372
00373 std::vector<std::pair<PointType, PointType> > & interface_segments =
00374 viennadata::access<InterfaceAreaKey, std::vector<std::pair<PointType, PointType> > >(interface_key)(*eit);
00375
00376
00377 PointType inner_point = (interface_segments[0].first + interface_segments[0].second) / 2.0;
00378 for (size_t i=1; i<interface_segments.size(); ++i)
00379 {
00380 inner_point += (interface_segments[i].first + interface_segments[i].second) / 2.0;
00381 }
00382 inner_point /= interface_segments.size();
00383
00384
00385
00386 double interface_area = 0.0;
00387 for (size_t i=0; i<interface_segments.size(); ++i)
00388 {
00389
00390 interface_area += spanned_volume(interface_segments[i].first, interface_segments[i].second, inner_point);
00391 }
00392
00393 viennadata::access<InterfaceAreaKey, double>(interface_key)(*eit) = interface_area;
00394 double volume_contribution = interface_area * edge_length / 6.0;
00395
00396 viennadata::access<BoxVolumeKey, double>(box_volume_key)(*eit) = 2.0 * volume_contribution;
00397 viennadata::access<BoxVolumeKey, double>(box_volume_key)(v0) += volume_contribution;
00398 viennadata::access<BoxVolumeKey, double>(box_volume_key)(v1) += volume_contribution;
00399 }
00400
00401 }
00402
00403
00404
00406 template <typename DomainType,
00407 typename InterfaceAreaKey,
00408 typename BoxVolumeKey>
00409 void write_voronoi_info(DomainType const & domain,
00410 InterfaceAreaKey const & interface_key,
00411 BoxVolumeKey const & box_volume_key,
00412 viennagrid::hexahedron_tag)
00413 {
00414
00415 typedef typename DomainType::config_type Config;
00416 typedef typename Config::cell_tag CellTag;
00417 typedef typename viennagrid::result_of::point<Config>::type PointType;
00418 typedef typename viennagrid::result_of::ncell<Config, 0>::type VertexType;
00419 typedef typename viennagrid::result_of::ncell<Config, 1>::type EdgeType;
00420 typedef typename viennagrid::result_of::ncell<Config, 2>::type FacetType;
00421 typedef typename viennagrid::result_of::ncell<Config, CellTag::dim>::type CellType;
00422
00423 typedef typename viennagrid::result_of::const_ncell_range<DomainType, CellTag::dim>::type CellRange;
00424 typedef typename viennagrid::result_of::iterator<CellRange>::type CellIterator;
00425
00426 typedef typename viennagrid::result_of::const_ncell_range<CellType, 2>::type FacetOnCellRange;
00427 typedef typename viennagrid::result_of::iterator<FacetOnCellRange>::type FacetOnCellIterator;
00428
00429 typedef typename viennagrid::result_of::const_ncell_range<FacetType, 1>::type EdgeOnFacetRange;
00430 typedef typename viennagrid::result_of::iterator<EdgeOnFacetRange>::type EdgeOnFacetIterator;
00431
00432 typedef typename viennagrid::result_of::const_ncell_range<EdgeType, 0>::type VertexOnEdgeRange;
00433 typedef typename viennagrid::result_of::iterator<VertexOnEdgeRange>::type VertexOnEdgeIterator;
00434
00435
00436
00437
00438
00439 CellRange cells = viennagrid::ncells<CellTag::dim>(domain);
00440 for (CellIterator cit = cells.begin();
00441 cit != cells.end();
00442 ++cit)
00443 {
00444 PointType cell_center = circumcenter(*cit);
00445
00446 FacetOnCellRange facets_on_cell = viennagrid::ncells<2>(*cit);
00447 for (FacetOnCellIterator focit = facets_on_cell.begin();
00448 focit != facets_on_cell.end();
00449 ++focit)
00450 {
00451 PointType facet_center = circumcenter(*focit);
00452
00453
00454 EdgeOnFacetRange edges_on_facet = viennagrid::ncells<1>(*focit);
00455 for (EdgeOnFacetIterator eocit = edges_on_facet.begin();
00456 eocit != edges_on_facet.end();
00457 ++eocit)
00458 {
00459 PointType edge_midpoint = viennagrid::circumcenter(*eocit);
00460
00461
00462 viennadata::access<InterfaceAreaKey, double>(interface_key)(*eocit) += spanned_volume(cell_center, facet_center, edge_midpoint);
00463
00464
00465 double edge_contribution = 0;
00466 VertexOnEdgeRange vertices_on_edge = viennagrid::ncells<0>(*eocit);
00467 for (VertexOnEdgeIterator voeit = vertices_on_edge.begin();
00468 voeit != vertices_on_edge.end();
00469 ++voeit)
00470 {
00471 double contribution = spanned_volume(cell_center, facet_center, edge_midpoint, voeit->point());
00472 edge_contribution += contribution;
00473 viennadata::access<BoxVolumeKey, double>(box_volume_key)(*voeit) += contribution;
00474 }
00475 viennadata::access<BoxVolumeKey, double>(box_volume_key)(*eocit) += edge_contribution;
00476 }
00477
00478 }
00479
00480 }
00481
00482 }
00483
00484 }
00485
00486
00487
00488
00489
00496 template <typename DomainType,
00497 typename InterfaceAreaKey,
00498 typename BoxVolumeKey>
00499 void apply_voronoi(DomainType const & domain,
00500 InterfaceAreaKey const & interface_area_key = viennagrid::voronoi_interface_area_key(),
00501 BoxVolumeKey const & box_volume_key = viennagrid::voronoi_box_volume_key())
00502 {
00503 detail::write_voronoi_info(domain,
00504 interface_area_key,
00505 box_volume_key,
00506 typename DomainType::config_type::cell_tag());
00507 }
00508
00510 template <typename DomainType>
00511 void apply_voronoi(DomainType const & domain)
00512 {
00513 apply_voronoi(domain, viennagrid::voronoi_interface_area_key(), viennagrid::voronoi_box_volume_key());
00514 }
00515
00516 }
00517 #endif