|
ViennaGrid - The Vienna Grid Library
2.1.0
|
00001 #ifndef VIENNAGRID_ALGORITHM_VORONOI_HPP 00002 #define VIENNAGRID_ALGORITHM_VORONOI_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 <iostream> 00017 #include <sstream> 00018 #include <string> 00019 #include <stdexcept> 00020 00021 #include "viennagrid/forwards.hpp" 00022 #include "viennagrid/mesh/segmentation.hpp" 00023 #include "viennagrid/mesh/coboundary_iteration.hpp" 00024 #include "viennagrid/algorithm/circumcenter.hpp" 00025 #include "viennagrid/algorithm/spanned_volume.hpp" 00026 #include "viennagrid/algorithm/volume.hpp" 00027 #include "viennagrid/algorithm/inner_prod.hpp" 00028 00029 00035 namespace viennagrid 00036 { 00037 namespace result_of 00038 { 00043 template <typename ConstCellHandleT> 00044 struct voronoi_cell_contribution 00045 { 00046 typedef std::vector< std::pair<ConstCellHandleT, double> > type; 00047 }; 00048 } 00049 00050 00051 namespace detail 00052 { 00055 template <typename ContainerT, typename PairT> 00056 void voronoi_unique_quantity_update(ContainerT & container, 00057 PairT const & update_pair) 00058 { 00059 bool found = false; 00060 for (typename ContainerT::iterator it = container.begin(); 00061 it != container.end(); 00062 ++it) 00063 { 00064 if (it->first == update_pair.first) 00065 { 00066 it->second += update_pair.second; 00067 found = true; 00068 break; 00069 } 00070 } 00071 00072 if (!found) 00073 { 00074 container.push_back(update_pair); 00075 } 00076 } 00077 00078 00080 template <typename CellTag, 00081 typename MeshT, 00082 typename InterfaceAreaAccessorT, 00083 typename InterfaceAreaCellContributionAccessorT, 00084 typename VertexBoxVolumeAccessorT, 00085 typename VertexBoxVolumeCellContributionAccessorT, 00086 typename EdgeBoxVolumeAccessorT, 00087 typename EdgeBoxVolumeCellContributionAccessorT> 00088 void write_voronoi_info(MeshT const & mesh_obj, 00089 InterfaceAreaAccessorT interface_area_accessor, 00090 InterfaceAreaCellContributionAccessorT interface_area_cell_contribution_accessor, 00091 VertexBoxVolumeAccessorT vertex_box_volume_accessor, 00092 VertexBoxVolumeCellContributionAccessorT vertex_box_volume_cell_contribution_accessor, 00093 EdgeBoxVolumeAccessorT edge_box_volume_accessor, 00094 EdgeBoxVolumeCellContributionAccessorT edge_box_volume_cell_contribution_accessor, 00095 viennagrid::simplex_tag<1>) 00096 { 00097 typedef typename viennagrid::result_of::element<MeshT, CellTag>::type CellType; 00098 00099 typedef typename viennagrid::result_of::element<MeshT, vertex_tag>::type VertexType; 00100 00101 typedef typename viennagrid::result_of::const_element_range<MeshT, CellTag>::type CellRange; 00102 typedef typename viennagrid::result_of::iterator<CellRange>::type CellIterator; 00103 00104 00105 typedef typename viennagrid::result_of::const_element_range<CellType, vertex_tag>::type VertexOnCellRange; 00106 typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type VertexOnCellIterator; 00107 00108 // 00109 // Write Voronoi information 00110 // 00111 00112 CellRange cells = viennagrid::elements<CellTag>(mesh_obj); 00113 for (CellIterator cit = cells.begin(); 00114 cit != cells.end(); 00115 ++cit) 00116 { 00117 interface_area_accessor(*cit) = 1; 00118 interface_area_cell_contribution_accessor(*cit).push_back( std::make_pair( cit.handle(), 1) ); 00119 00120 double edge_contribution = 0; 00121 VertexOnCellRange vertices_on_cell = viennagrid::elements<VertexType>(*cit); 00122 for (VertexOnCellIterator vocit = vertices_on_cell.begin(); 00123 vocit != vertices_on_cell.end(); 00124 ++vocit) 00125 { 00126 double contribution = volume(*cit) / 2.0; 00127 edge_contribution += contribution; 00128 00129 vertex_box_volume_accessor(*vocit) += contribution; 00130 vertex_box_volume_cell_contribution_accessor(*vocit).push_back( std::make_pair( cit.handle(), contribution) ); 00131 } 00132 00133 edge_box_volume_accessor(*cit) = edge_contribution; 00134 edge_box_volume_cell_contribution_accessor(*cit).push_back( std::make_pair( cit.handle(), edge_contribution) ); 00135 00136 } 00137 } 00138 00140 template <typename CellTag, 00141 typename MeshT, 00142 typename InterfaceAreaAccessorT, 00143 typename InterfaceAreaCellContributionAccessorT, 00144 typename VertexBoxVolumeAccessorT, 00145 typename VertexBoxVolumeCellContributionAccessorT, 00146 typename EdgeBoxVolumeAccessorT, 00147 typename EdgeBoxVolumeCellContributionAccessorT> 00148 void write_voronoi_info(MeshT const & mesh_obj, 00149 InterfaceAreaAccessorT interface_area_accessor, 00150 InterfaceAreaCellContributionAccessorT interface_area_cell_contribution_accessor, 00151 VertexBoxVolumeAccessorT vertex_box_volume_accessor, 00152 VertexBoxVolumeCellContributionAccessorT vertex_box_volume_cell_contribution_accessor, 00153 EdgeBoxVolumeAccessorT edge_box_volume_accessor, 00154 EdgeBoxVolumeCellContributionAccessorT edge_box_volume_cell_contribution_accessor, 00155 viennagrid::hypercube_tag<1>) 00156 { 00157 write_voronoi_info<CellTag>(mesh_obj, 00158 interface_area_accessor, interface_area_cell_contribution_accessor, 00159 vertex_box_volume_accessor, vertex_box_volume_cell_contribution_accessor, 00160 edge_box_volume_accessor, edge_box_volume_cell_contribution_accessor, 00161 viennagrid::simplex_tag<1>()); 00162 } 00163 00164 // 00165 // Voronoi information in two (topological) dimensions 00166 // 00168 template <typename CellTag, 00169 typename MeshT, 00170 typename InterfaceAreaAccessorT, 00171 typename InterfaceAreaCellContributionAccessorT, 00172 typename VertexBoxVolumeAccessorT, 00173 typename VertexBoxVolumeCellContributionAccessorT, 00174 typename EdgeBoxVolumeAccessorT, 00175 typename EdgeBoxVolumeCellContributionAccessorT> 00176 void write_voronoi_info(MeshT const & mesh_obj, 00177 InterfaceAreaAccessorT interface_area_accessor, 00178 InterfaceAreaCellContributionAccessorT interface_area_cell_contribution_accessor, 00179 VertexBoxVolumeAccessorT vertex_box_volume_accessor, 00180 VertexBoxVolumeCellContributionAccessorT vertex_box_volume_cell_contribution_accessor, 00181 EdgeBoxVolumeAccessorT edge_box_volume_accessor, 00182 EdgeBoxVolumeCellContributionAccessorT edge_box_volume_cell_contribution_accessor, 00183 viennagrid::quadrilateral_tag) 00184 { 00185 typedef typename viennagrid::result_of::element<MeshT, CellTag>::type CellType; 00186 00187 typedef typename viennagrid::result_of::point<MeshT>::type PointType; 00188 typedef typename viennagrid::result_of::element<MeshT, vertex_tag>::type VertexType; 00189 typedef typename viennagrid::result_of::element<MeshT, line_tag>::type EdgeType; 00190 00191 typedef typename viennagrid::result_of::const_element_range<MeshT, CellTag>::type CellRange; 00192 typedef typename viennagrid::result_of::iterator<CellRange>::type CellIterator; 00193 00194 typedef typename viennagrid::result_of::const_element_range<CellType, line_tag>::type EdgeOnCellRange; 00195 typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type EdgeOnCellIterator; 00196 00197 typedef typename viennagrid::result_of::const_element_range<EdgeType, vertex_tag>::type VertexOnEdgeRange; 00198 typedef typename viennagrid::result_of::iterator<VertexOnEdgeRange>::type VertexOnEdgeIterator; 00199 00200 00201 // 00202 // Algorithm: Iterate over all cells, compute circumcenter and add interface area to edge, box volume to vertex. 00203 // 00204 00205 CellRange cells = viennagrid::elements<CellTag>(mesh_obj); 00206 for (CellIterator cit = cells.begin(); 00207 cit != cells.end(); 00208 ++cit) 00209 { 00210 PointType circ_center = circumcenter(*cit); 00211 00212 //iterate over edges: 00213 EdgeOnCellRange edges_on_cell = viennagrid::elements<EdgeType>(*cit); 00214 for (EdgeOnCellIterator eocit = edges_on_cell.begin(); 00215 eocit != edges_on_cell.end(); 00216 ++eocit) 00217 { 00218 PointType edge_midpoint = circumcenter(*eocit); 00219 00220 // interface contribution: 00221 double interface_contribution = spanned_volume(circ_center, edge_midpoint); 00222 00223 interface_area_accessor(*eocit) += interface_contribution; 00224 interface_area_cell_contribution_accessor(*eocit).push_back( std::make_pair( cit.handle(), interface_contribution) ); 00225 00226 //box volume contribution: 00227 double edge_contribution = 0; 00228 VertexOnEdgeRange vertices_on_edge = viennagrid::elements<VertexType>(*eocit); 00229 for (VertexOnEdgeIterator voeit = vertices_on_edge.begin(); 00230 voeit != vertices_on_edge.end(); 00231 ++voeit) 00232 { 00233 double contribution = spanned_volume(circ_center, edge_midpoint, viennagrid::point(*voeit)); 00234 edge_contribution += contribution; 00235 00236 vertex_box_volume_accessor(*voeit) += contribution; 00237 vertex_box_volume_cell_contribution_accessor(*voeit).push_back( std::make_pair( cit.handle(), contribution) ); 00238 } 00239 00240 edge_box_volume_accessor(*eocit) += edge_contribution; 00241 edge_box_volume_cell_contribution_accessor(*eocit).push_back( std::make_pair( cit.handle(), edge_contribution) ); 00242 00243 } //for edges on cells 00244 00245 } //for cells 00246 00247 } //write_voronoi_info(triangle_tag) 00248 00249 00251 //template <typename PointType, typename ConfigType> 00252 //PointType point_to_local_coordinates(PointType const & p, viennagrid::element<ConfigType, viennagrid::simplex_tag<2> > const & triangle) 00253 template<typename TriangleType> 00254 typename viennagrid::result_of::point<TriangleType>::type point_to_local_coordinates( typename viennagrid::result_of::point<TriangleType>::type const & p, 00255 const TriangleType & triangle) 00256 { 00257 typedef TriangleType CellType; 00258 typedef typename viennagrid::result_of::const_element_range<CellType, vertex_tag>::type VertexRange; 00259 typedef typename viennagrid::result_of::point<TriangleType>::type PointType; 00260 00261 VertexRange vertices = viennagrid::elements<viennagrid::vertex_tag>(triangle); 00262 00263 PointType const & a = viennagrid::point(vertices[0]); 00264 PointType const & b = viennagrid::point(vertices[1]); 00265 PointType const & c = viennagrid::point(vertices[2]); 00266 00267 PointType v0 = b - a; 00268 PointType v1 = c - a; 00269 PointType vp = p - a; 00270 00271 double det = v0[0] * v1[1] - v0[1] * v1[0]; 00272 00273 return PointType( (v1[1] * vp[0] - v1[0] * vp[1]) / det, 00274 (v0[0] * vp[1] - v0[1] * vp[0]) / det ); 00275 } 00276 00277 00279 template <typename PointType> 00280 PointType line_intersection(PointType const & a, PointType const & b, 00281 PointType const & c, PointType const & d) 00282 { 00283 double det = (a[0] - b[0])*(c[1] - d[1]) - (a[1] - b[1])*(c[0] - d[0]); 00284 00285 if (std::fabs(det) < 1e-10 * viennagrid::norm(b-a)) 00286 std::cerr << "viennagrid::detail::line_intersection(): Warning: Determinant is close to zero!" << std::endl; 00287 00288 double px = (a[0]*b[1] - a[1]*b[0]) * (c[0] - d[0]) 00289 - (a[0] - b[0]) * (c[0]*d[1] - c[1]*d[0]); 00290 double py = (a[0]*b[1] - a[1]*b[0]) * (c[1] - d[1]) 00291 - (a[1] - b[1]) * (c[0]*d[1] - c[1]*d[0]); 00292 00293 return PointType(px / det, py / det); 00294 } 00295 00296 00298 template <typename CellTag, 00299 typename MeshT, 00300 typename InterfaceAreaAccessorT, 00301 typename InterfaceAreaCellContributionAccessorT, 00302 typename VertexBoxVolumeAccessorT, 00303 typename VertexBoxVolumeCellContributionAccessorT, 00304 typename EdgeBoxVolumeAccessorT, 00305 typename EdgeBoxVolumeCellContributionAccessorT> 00306 void write_voronoi_info(MeshT const & mesh_obj, 00307 InterfaceAreaAccessorT interface_area_accessor, 00308 InterfaceAreaCellContributionAccessorT interface_area_cell_contribution_accessor, 00309 VertexBoxVolumeAccessorT vertex_box_volume_accessor, 00310 VertexBoxVolumeCellContributionAccessorT vertex_box_volume_cell_contribution_accessor, 00311 EdgeBoxVolumeAccessorT edge_box_volume_accessor, 00312 EdgeBoxVolumeCellContributionAccessorT edge_box_volume_cell_contribution_accessor, 00313 viennagrid::triangle_tag) 00314 { 00315 typedef typename viennagrid::result_of::element<MeshT, CellTag>::type CellType; 00316 typedef typename viennagrid::result_of::const_handle<MeshT, CellTag>::type ConstCellHandleType; 00317 00318 typedef typename viennagrid::result_of::point<MeshT>::type PointType; 00319 typedef typename viennagrid::result_of::element<MeshT, vertex_tag>::type VertexType; 00320 typedef typename viennagrid::result_of::element<MeshT, line_tag>::type EdgeType; 00321 00322 typedef typename viennagrid::result_of::const_element_range<MeshT, CellTag>::type CellRange; 00323 typedef typename viennagrid::result_of::iterator<CellRange>::type CellIterator; 00324 00325 typedef typename viennagrid::result_of::const_element_range<CellType, vertex_tag>::type VertexOnCellRange; 00326 00327 typedef typename viennagrid::result_of::const_coboundary_range<MeshT, EdgeType, CellTag>::type CellOnEdgeRange; 00328 00329 // typedef typename viennagrid::result_of::const_element_range<EdgeType, CellTag>::type CellOnEdgeRange; 00330 typedef typename viennagrid::result_of::iterator<CellOnEdgeRange>::type CellOnEdgeIterator; 00331 00332 typedef typename viennagrid::result_of::const_element_range<CellType, line_tag>::type EdgeOnCellRange; 00333 typedef typename viennagrid::result_of::iterator<EdgeOnCellRange>::type EdgeOnCellIterator; 00334 00335 typedef typename viennagrid::result_of::const_element_range<EdgeType, vertex_tag>::type VertexOnEdgeRange; 00336 typedef typename viennagrid::result_of::iterator<VertexOnEdgeRange>::type VertexOnEdgeIterator; 00337 00338 CellRange cells = viennagrid::elements<CellTag>(mesh_obj); 00339 for (CellIterator cit = cells.begin(); 00340 cit != cells.end(); 00341 ++cit) 00342 { 00343 PointType circ_center = circumcenter(*cit); 00344 00345 00346 PointType circ_center_local = point_to_local_coordinates(circ_center, *cit); 00347 00348 if ( (circ_center_local[0] < 0) 00349 || (circ_center_local[1] < 0) 00350 || (circ_center_local[0] + circ_center_local[1] > 1.0) ) //circumcenter is outside triangle - things get complicated 00351 { 00352 EdgeOnCellRange edges_on_cell = viennagrid::elements<EdgeType>(*cit); 00353 VertexOnCellRange vertices_on_cell = viennagrid::elements<VertexType>(*cit); 00354 00355 // 00356 // Step 1: Get intersected edge, find opposite vertex and the 'other' triangle 00357 // 00358 EdgeType const * intersected_edge_ptr = NULL; 00359 VertexType const * opposite_vertex_ptr = NULL; 00360 if (circ_center_local[1] < 0) 00361 { 00362 intersected_edge_ptr = &(edges_on_cell[0]); 00363 opposite_vertex_ptr = &(vertices_on_cell[2]); 00364 } 00365 else if (circ_center_local[0] < 0) 00366 { 00367 intersected_edge_ptr = &(edges_on_cell[1]); 00368 opposite_vertex_ptr = &(vertices_on_cell[1]); 00369 } 00370 else 00371 { 00372 intersected_edge_ptr = &(edges_on_cell[2]); 00373 opposite_vertex_ptr = &(vertices_on_cell[0]); 00374 } 00375 00376 //find 'other' triangle 00377 ConstCellHandleType other_cell; 00378 00379 // CellType const * other_cell = NULL; 00380 CellOnEdgeRange other_cells = viennagrid::coboundary_elements<EdgeType, CellTag>(mesh_obj, viennagrid::handle(mesh_obj, *intersected_edge_ptr) ); 00381 detail::set_handle_invalid( other_cells, other_cell ); 00382 00383 for (CellOnEdgeIterator coeit = other_cells.begin(); 00384 coeit != other_cells.end(); 00385 ++coeit) 00386 { 00387 if ( coeit->id() != cit->id() ) 00388 { 00389 other_cell = coeit.handle(); 00390 break; 00391 } 00392 } 00393 00394 // 00395 // Step 2: Precompute intersection of opposite vertex with intersected edge 00396 // 00397 VertexOnEdgeRange vertices_on_intersected_edge = viennagrid::elements<VertexType>(*intersected_edge_ptr); 00398 PointType opposite_vertex_edge_intersection = line_intersection( viennagrid::point(*opposite_vertex_ptr), circ_center, 00399 viennagrid::point(vertices_on_intersected_edge[0]), 00400 viennagrid::point(vertices_on_intersected_edge[1])); 00401 00402 // 00403 // Step 3: Compute contributions 00404 // 00405 for (EdgeOnCellIterator eocit = edges_on_cell.begin(); 00406 eocit != edges_on_cell.end(); 00407 ++eocit) 00408 { 00409 PointType edge_midpoint = circumcenter(*eocit); 00410 00411 if ( intersected_edge_ptr != &(*eocit) ) // non-intersected edge: Split contributions into contributions from cell itself and contribution from other cell 00412 { 00413 PointType edge_intersection = line_intersection( edge_midpoint, circ_center, 00414 viennagrid::point(vertices_on_intersected_edge[0]), 00415 viennagrid::point(vertices_on_intersected_edge[1])); 00416 00417 // 00418 // Interface contributions (two contributions: From edge midpoint to intersected edge, and from intersected edge to circumcenter) 00419 // 00420 00421 double interface_contribution = 0; 00422 00423 // contribution from cell: 00424 interface_contribution = spanned_volume(edge_intersection, edge_midpoint); 00425 00426 interface_area_accessor(*eocit) += interface_contribution; 00427 voronoi_unique_quantity_update(interface_area_cell_contribution_accessor(*eocit), std::make_pair( cit.handle(), interface_contribution) ); 00428 00429 // contribution from other cell containing the circumcenter: 00430 interface_contribution = spanned_volume(edge_intersection, circ_center); 00431 interface_area_accessor(*eocit) += interface_contribution; 00432 //if ( other_cell != NULL) 00433 if (!detail::is_handle_invalid( other_cells, other_cell )) 00434 voronoi_unique_quantity_update(interface_area_cell_contribution_accessor(*eocit), std::make_pair( other_cell, interface_contribution) ); 00435 00436 00437 // 00438 // Box volume contribution: 00439 // 00440 double edge_contribution = 0; 00441 00442 // contribution from cell itself: 00443 VertexOnEdgeRange vertices_on_edge = viennagrid::elements<VertexType>(*eocit); 00444 for (VertexOnEdgeIterator voeit = vertices_on_edge.begin(); 00445 voeit != vertices_on_edge.end(); 00446 ++voeit) 00447 { 00448 double contribution = spanned_volume(edge_intersection, edge_midpoint, viennagrid::point(*voeit)); 00449 edge_contribution += contribution; 00450 00451 vertex_box_volume_accessor(*voeit) += contribution; 00452 voronoi_unique_quantity_update(vertex_box_volume_cell_contribution_accessor(*voeit), std::make_pair( cit.handle(), contribution) ); 00453 00454 if ( &(*voeit) != opposite_vertex_ptr ) // non-splitted contribution 00455 { 00456 // if (other_cell != NULL) 00457 if (!detail::is_handle_invalid( other_cells, other_cell )) 00458 { 00459 double contribution_other = spanned_volume(circ_center, edge_intersection, viennagrid::point(*voeit)); 00460 vertex_box_volume_accessor(*voeit) += contribution_other; 00461 voronoi_unique_quantity_update(vertex_box_volume_cell_contribution_accessor(*voeit), 00462 std::make_pair(other_cell, contribution_other) ); 00463 edge_box_volume_accessor(*eocit) += contribution_other; 00464 voronoi_unique_quantity_update(edge_box_volume_cell_contribution_accessor(*eocit), 00465 std::make_pair(other_cell, contribution_other) ); 00466 } 00467 } 00468 else //splitted contribution around center of the edge. Contributes to opposite vertex 00469 { 00470 double contribution_cell = spanned_volume(opposite_vertex_edge_intersection, edge_intersection, viennagrid::point(*voeit)); 00471 vertex_box_volume_accessor(*voeit) += contribution_cell; 00472 voronoi_unique_quantity_update(vertex_box_volume_cell_contribution_accessor(*voeit), 00473 std::make_pair( cit.handle(), contribution_cell) ); 00474 edge_box_volume_accessor(*eocit) += contribution_cell; 00475 voronoi_unique_quantity_update(edge_box_volume_cell_contribution_accessor(*eocit), 00476 std::make_pair( cit.handle(), contribution_cell) ); 00477 00478 // if (other_cell != NULL) 00479 if (!detail::is_handle_invalid( other_cells, other_cell )) 00480 { 00481 double contribution_other = spanned_volume(circ_center, edge_intersection, opposite_vertex_edge_intersection); 00482 vertex_box_volume_accessor(*voeit) += contribution_other; 00483 voronoi_unique_quantity_update(vertex_box_volume_cell_contribution_accessor(*voeit), 00484 std::make_pair(other_cell, contribution_other) ); 00485 edge_box_volume_accessor(*eocit) += contribution_other; 00486 voronoi_unique_quantity_update(edge_box_volume_cell_contribution_accessor(*eocit), 00487 std::make_pair(other_cell, contribution_other) ); 00488 } 00489 } 00490 00491 00492 } 00493 edge_box_volume_accessor(*eocit) += edge_contribution; 00494 voronoi_unique_quantity_update(edge_box_volume_cell_contribution_accessor(*eocit), 00495 std::make_pair( cit.handle(), edge_contribution) ); 00496 00497 } 00498 // else if (other_cell != NULL)// intersected edge: Write negative contributions to other cell 00499 else if (!detail::is_handle_invalid( other_cells, other_cell )) 00500 { 00501 // interface contribution: 00502 double interface_contribution = spanned_volume(circ_center, edge_midpoint); 00503 interface_area_accessor(*eocit) -= interface_contribution; 00504 voronoi_unique_quantity_update(interface_area_cell_contribution_accessor(*eocit), 00505 std::make_pair(other_cell, -1.0 * interface_contribution) ); 00506 00507 00508 //box volume contribution: 00509 double edge_contribution = 0; 00510 VertexOnEdgeRange vertices_on_edge = viennagrid::elements<VertexType>(*eocit); 00511 for (VertexOnEdgeIterator voeit = vertices_on_edge.begin(); 00512 voeit != vertices_on_edge.end(); 00513 ++voeit) 00514 { 00515 double contribution = spanned_volume(circ_center, edge_midpoint, viennagrid::point(*voeit)); 00516 edge_contribution += contribution; 00517 vertex_box_volume_accessor(*voeit) -= contribution; 00518 voronoi_unique_quantity_update(vertex_box_volume_cell_contribution_accessor(*voeit), 00519 std::make_pair(other_cell, -1.0 * contribution) ); 00520 } 00521 edge_box_volume_accessor(*eocit) -= edge_contribution; 00522 voronoi_unique_quantity_update(edge_box_volume_cell_contribution_accessor(*eocit), 00523 std::make_pair(other_cell, -1.0 * edge_contribution) ); 00524 } 00525 // else {} //nothing to do for a boundary edge 00526 00527 } //for edges on cells 00528 00529 } 00530 else // circumcenter is inside triangle (everything is simple... ;-) ) 00531 { 00532 00533 //iterate over edges: 00534 EdgeOnCellRange edges_on_cell = viennagrid::elements<EdgeType>(*cit); 00535 for (EdgeOnCellIterator eocit = edges_on_cell.begin(); 00536 eocit != edges_on_cell.end(); 00537 ++eocit) 00538 { 00539 PointType edge_midpoint = circumcenter(*eocit); 00540 00541 // interface contribution: 00542 double interface_contribution = spanned_volume(circ_center, edge_midpoint); 00543 interface_area_accessor(*eocit) += interface_contribution; 00544 voronoi_unique_quantity_update(interface_area_cell_contribution_accessor(*eocit), 00545 std::make_pair( cit.handle(), interface_contribution) ); 00546 00547 00548 //box volume contribution: 00549 double edge_contribution = 0; 00550 VertexOnEdgeRange vertices_on_edge = viennagrid::elements<VertexType>(*eocit); 00551 for (VertexOnEdgeIterator voeit = vertices_on_edge.begin(); 00552 voeit != vertices_on_edge.end(); 00553 ++voeit) 00554 { 00555 double contribution = spanned_volume(circ_center, edge_midpoint, viennagrid::point(*voeit)); 00556 edge_contribution += contribution; 00557 vertex_box_volume_accessor(*voeit) += contribution; 00558 voronoi_unique_quantity_update(vertex_box_volume_cell_contribution_accessor(*voeit), 00559 std::make_pair( cit.handle(), contribution) ); 00560 } 00561 edge_box_volume_accessor(*eocit) += edge_contribution; 00562 voronoi_unique_quantity_update(edge_box_volume_cell_contribution_accessor(*eocit), 00563 std::make_pair( cit.handle(), edge_contribution) ); 00564 } //for edges on cells 00565 00566 } 00567 00568 } //for cells 00569 00570 } 00571 00572 00573 // 00574 // Voronoi information in three dimensions 00575 // 00576 00585 template <typename CellTag, 00586 typename MeshT, 00587 typename InterfaceAreaAccessorT, 00588 typename InterfaceAreaCellContributionAccessorT, 00589 typename VertexBoxVolumeAccessorT, 00590 typename VertexBoxVolumeCellContributionAccessorT, 00591 typename EdgeBoxVolumeAccessorT, 00592 typename EdgeBoxVolumeCellContributionAccessorT> 00593 void write_voronoi_info(MeshT const & mesh_obj, 00594 InterfaceAreaAccessorT interface_area_accessor, 00595 InterfaceAreaCellContributionAccessorT interface_area_cell_contribution_accessor, 00596 VertexBoxVolumeAccessorT vertex_box_volume_accessor, 00597 VertexBoxVolumeCellContributionAccessorT vertex_box_volume_cell_contribution_accessor, 00598 EdgeBoxVolumeAccessorT edge_box_volume_accessor, 00599 EdgeBoxVolumeCellContributionAccessorT edge_box_volume_cell_contribution_accessor, 00600 viennagrid::tetrahedron_tag) 00601 { 00602 typedef typename viennagrid::result_of::element<MeshT, CellTag>::type CellType; 00603 typedef typename viennagrid::result_of::const_handle<MeshT, CellTag>::type ConstCellHandleType; 00604 00605 typedef typename viennagrid::result_of::point<MeshT>::type PointType; 00606 typedef typename viennagrid::result_of::element<MeshT, vertex_tag>::type VertexType; 00607 typedef typename viennagrid::result_of::element<MeshT, line_tag>::type EdgeType; 00608 typedef typename viennagrid::result_of::const_handle<MeshT, line_tag>::type ConstEdgeHandleType; 00609 typedef typename viennagrid::result_of::element<MeshT, triangle_tag>::type FacetType; 00610 typedef typename viennagrid::result_of::const_handle<MeshT, triangle_tag>::type ConstFacetHandleType; 00611 00612 typedef typename viennagrid::result_of::coord<PointType>::type CoordType; 00613 00614 typedef typename viennagrid::result_of::const_element_range<MeshT, CellTag>::type CellRange; 00615 typedef typename viennagrid::result_of::iterator<CellRange>::type CellIterator; 00616 00617 typedef typename viennagrid::result_of::const_element_range<MeshT, typename CellTag::facet_tag>::type FacetRange; 00618 typedef typename viennagrid::result_of::iterator<FacetRange>::type FacetIterator; 00619 00620 typedef typename viennagrid::result_of::const_element_range<MeshT, line_tag>::type EdgeRange; 00621 typedef typename viennagrid::result_of::iterator<EdgeRange>::type EdgeIterator; 00622 00623 typedef typename viennagrid::result_of::const_element_range<CellType, typename CellTag::facet_tag>::type FacetOnCellRange; 00624 typedef typename viennagrid::result_of::iterator<FacetOnCellRange>::type FacetOnCellIterator; 00625 00626 typedef typename viennagrid::result_of::const_element_range<FacetType, line_tag>::type EdgeOnFacetRange; 00627 typedef typename viennagrid::result_of::iterator<EdgeOnFacetRange>::type EdgeOnFacetIterator; 00628 00629 typedef typename viennagrid::result_of::const_element_range<EdgeType, vertex_tag>::type VertexOnEdgeRange; 00630 typedef typename viennagrid::result_of::iterator<VertexOnEdgeRange>::type VertexOnEdgeIterator; 00631 00632 typedef std::pair<PointType, ConstCellHandleType> PointWithCellInfo; 00633 typedef std::pair<std::pair<PointType, PointType>, 00634 ConstCellHandleType> EdgePointsWithCellInfo; 00635 typedef std::vector< PointWithCellInfo > CircumcenterContainer; 00636 00637 00638 // 00639 // Step one: Write circumcenters to facets 00640 // 00641 00642 viennagrid::detail::dereference_handle_comparator<MeshT> comp(mesh_obj); 00643 00644 std::map< ConstFacetHandleType, CircumcenterContainer, viennagrid::detail::dereference_handle_comparator<MeshT> > 00645 circumcenters_on_facets( comp ); 00646 //std::map< EdgeType const *, 00647 std::map< ConstEdgeHandleType, 00648 std::vector< EdgePointsWithCellInfo >, 00649 viennagrid::detail::dereference_handle_comparator<MeshT> 00650 > interface_boundaries_on_edges(comp); 00651 00652 00653 CellRange cells = viennagrid::elements<CellType>(mesh_obj); 00654 for (CellIterator cit = cells.begin(); 00655 cit != cells.end(); 00656 ++cit) 00657 { 00658 const CellType & cell = *cit; 00659 PointType circ_center = circumcenter(cell); 00660 00661 FacetOnCellRange facets_on_cell = viennagrid::elements<FacetType>(cell); 00662 for (FacetOnCellIterator focit = facets_on_cell.begin(); 00663 focit != facets_on_cell.end(); 00664 ++focit) 00665 { 00666 circumcenters_on_facets[focit.handle()].push_back( PointWithCellInfo(circ_center, cit.handle()) ); 00667 } //for edges on cells 00668 } //for cells 00669 00670 00671 // 00672 // Step two: Write lines connecting circumcenters to edges 00673 // 00674 FacetRange facets = viennagrid::elements<FacetType>(mesh_obj); 00675 for (FacetIterator fit = facets.begin(); 00676 fit != facets.end(); 00677 ++fit) 00678 { 00679 CircumcenterContainer & circ_centers = circumcenters_on_facets[fit.handle()]; 00680 const FacetType & facet = *fit; 00681 00682 EdgeOnFacetRange edges_on_facet = viennagrid::elements<EdgeType>( facet ); 00683 for (EdgeOnFacetIterator eofit = edges_on_facet.begin(); 00684 eofit != edges_on_facet.end(); 00685 ++eofit) 00686 { 00687 const EdgeType & edge = *eofit; 00688 00689 if (circ_centers.size() == 1) 00690 { 00691 interface_boundaries_on_edges[eofit.handle()].push_back(std::make_pair( std::make_pair(circ_centers[0].first, circumcenter(facet)), 00692 circ_centers[0].second) 00693 ); 00694 interface_boundaries_on_edges[eofit.handle()].push_back(std::make_pair( std::make_pair(circumcenter(edge), circumcenter(facet)), 00695 circ_centers[0].second) 00696 ); 00697 } 00698 else if (circ_centers.size() == 2) 00699 { 00700 PointType edge_mid = circ_centers[0].first + circ_centers[1].first; 00701 edge_mid /= 2.0; 00702 00703 interface_boundaries_on_edges[eofit.handle()].push_back(std::make_pair( std::make_pair(circ_centers[0].first, edge_mid), 00704 circ_centers[0].second) 00705 ); 00706 interface_boundaries_on_edges[eofit.handle()].push_back(std::make_pair( std::make_pair(edge_mid, circ_centers[1].first), 00707 circ_centers[1].second) 00708 ); 00709 } 00710 else 00711 { 00712 std::cerr << "circ_centers.size() = " << circ_centers.size() << std::endl; 00713 std::cerr << "*fit: " << facet << std::endl; 00714 throw std::runtime_error("More than two circumcenters for a facet in three dimensions!"); 00715 } 00716 00717 } //for edges on cells 00718 } 00719 00720 00721 // 00722 // Step three: Compute Voronoi information: 00723 // 00724 EdgeRange edges = viennagrid::elements<EdgeType>(mesh_obj); 00725 for (EdgeIterator eit = edges.begin(); 00726 eit != edges.end(); 00727 ++eit) 00728 { 00729 const EdgeType & edge = *eit; 00730 00731 //get vertices of edge: 00732 VertexOnEdgeRange vertices_on_edge = viennagrid::elements<VertexType>(edge); 00733 VertexOnEdgeIterator voeit = vertices_on_edge.begin(); 00734 00735 00736 00737 VertexType const & v0 = *voeit; 00738 ++voeit; 00739 VertexType const & v1 = *voeit; 00740 00741 double edge_length = spanned_volume( viennagrid::point(mesh_obj, v0), viennagrid::point(mesh_obj, v1)); 00742 00743 std::vector< EdgePointsWithCellInfo > & interface_segments = interface_boundaries_on_edges[eit.handle()]; 00744 00745 // 00746 // determine inner point of convex interface polygon: 00747 // 00748 PointType inner_point = (interface_segments[0].first.first + interface_segments[0].first.second) / 2.0; 00749 for (std::size_t i=1; i<interface_segments.size(); ++i) 00750 { 00751 inner_point += (interface_segments[i].first.first + interface_segments[i].first.second) / 2.0; 00752 } 00753 inner_point /= static_cast<CoordType>(interface_segments.size()); 00754 00755 // 00756 // compute interface area 00757 // 00758 double interface_area = 0.0; 00759 for (std::size_t i=0; i<interface_segments.size(); ++i) 00760 { 00761 // interface area: 00762 double interface_contribution = spanned_volume(interface_segments[i].first.first, interface_segments[i].first.second, inner_point); 00763 if (interface_contribution > 0) 00764 { 00765 voronoi_unique_quantity_update(interface_area_cell_contribution_accessor( edge ), 00766 std::make_pair(interface_segments[i].second, interface_contribution) ); 00767 interface_area += interface_contribution; 00768 00769 // box volume: 00770 double volume_contribution = interface_contribution * edge_length / 6.0; 00771 voronoi_unique_quantity_update(edge_box_volume_cell_contribution_accessor( edge ), 00772 std::make_pair(interface_segments[i].second, 2.0 * volume_contribution) ); //volume contribution of both box volumes associated with the edge 00773 voronoi_unique_quantity_update(vertex_box_volume_cell_contribution_accessor(v0), 00774 std::make_pair( interface_segments[i].second, volume_contribution) ); 00775 voronoi_unique_quantity_update(vertex_box_volume_cell_contribution_accessor(v1), 00776 std::make_pair( interface_segments[i].second, volume_contribution) ); 00777 } 00778 } 00779 00780 // 00781 // Write Voronoi info: 00782 // 00783 00784 //std::cout << "Interface area: " << interface_area << std::endl; 00785 interface_area_accessor( edge ) = interface_area; 00786 double volume_contribution = interface_area * edge_length / 6.0; 00787 edge_box_volume_accessor(edge) = 2.0 * volume_contribution; //volume contribution of both box volumes associated with the edge 00788 vertex_box_volume_accessor(v0) += volume_contribution; 00789 vertex_box_volume_accessor(v1) += volume_contribution; 00790 00791 } //for edges 00792 00793 } //write_voronoi_info(tetrahedron_tag) 00794 00795 00796 00798 template <typename CellTag, 00799 typename MeshT, 00800 typename InterfaceAreaAccessorT, 00801 typename InterfaceAreaCellContributionAccessorT, 00802 typename VertexBoxVolumeAccessorT, 00803 typename VertexBoxVolumeCellContributionAccessorT, 00804 typename EdgeBoxVolumeAccessorT, 00805 typename EdgeBoxVolumeCellContributionAccessorT> 00806 void write_voronoi_info(MeshT const & mesh_obj, 00807 InterfaceAreaAccessorT interface_area_accessor, 00808 InterfaceAreaCellContributionAccessorT interface_area_cell_contribution_accessor, 00809 VertexBoxVolumeAccessorT vertex_box_volume_accessor, 00810 VertexBoxVolumeCellContributionAccessorT vertex_box_volume_cell_contribution_accessor, 00811 EdgeBoxVolumeAccessorT edge_box_volume_accessor, 00812 EdgeBoxVolumeCellContributionAccessorT edge_box_volume_cell_contribution_accessor, 00813 viennagrid::hexahedron_tag) 00814 { 00815 //std::cout << "Warning: Voronoi info for hexahedron is only correct when having regular cuboids only." << std::endl; 00816 00817 typedef typename viennagrid::result_of::element<MeshT, CellTag>::type CellType; 00818 00819 typedef typename viennagrid::result_of::point<MeshT>::type PointType; 00820 typedef typename viennagrid::result_of::element<MeshT, vertex_tag>::type VertexType; 00821 typedef typename viennagrid::result_of::element<MeshT, line_tag>::type EdgeType; 00822 typedef typename viennagrid::result_of::element<MeshT, quadrilateral_tag>::type FacetType; 00823 00824 typedef typename viennagrid::result_of::const_element_range<MeshT, CellTag>::type CellRange; 00825 typedef typename viennagrid::result_of::iterator<CellRange>::type CellIterator; 00826 00827 typedef typename viennagrid::result_of::const_element_range<CellType, quadrilateral_tag>::type FacetOnCellRange; 00828 typedef typename viennagrid::result_of::iterator<FacetOnCellRange>::type FacetOnCellIterator; 00829 00830 typedef typename viennagrid::result_of::const_element_range<FacetType, line_tag>::type EdgeOnFacetRange; 00831 typedef typename viennagrid::result_of::iterator<EdgeOnFacetRange>::type EdgeOnFacetIterator; 00832 00833 typedef typename viennagrid::result_of::const_element_range<EdgeType, vertex_tag>::type VertexOnEdgeRange; 00834 typedef typename viennagrid::result_of::iterator<VertexOnEdgeRange>::type VertexOnEdgeIterator; 00835 00836 // 00837 // Algorithm: Iterate over all cells, compute circumcenter and add interface area to edge, box volume to vertex. 00838 // 00839 00840 CellRange cells = viennagrid::elements<CellTag>(mesh_obj); 00841 for (CellIterator cit = cells.begin(); 00842 cit != cells.end(); 00843 ++cit) 00844 { 00845 const CellType & cell = *cit; 00846 PointType cell_center = circumcenter(cell); 00847 00848 FacetOnCellRange facets_on_cell = viennagrid::elements<viennagrid::quadrilateral_tag>(cell); 00849 for (FacetOnCellIterator focit = facets_on_cell.begin(); 00850 focit != facets_on_cell.end(); 00851 ++focit) 00852 { 00853 PointType facet_center = circumcenter(*focit); 00854 00855 //iterate over edges: 00856 EdgeOnFacetRange edges_on_facet = viennagrid::elements<EdgeType>(*focit); 00857 for (EdgeOnFacetIterator eocit = edges_on_facet.begin(); 00858 eocit != edges_on_facet.end(); 00859 ++eocit) 00860 { 00861 PointType edge_midpoint = viennagrid::circumcenter(*eocit); 00862 00863 // interface contribution: 00864 double interface_contribution = spanned_volume(cell_center, facet_center, edge_midpoint); 00865 interface_area_cell_contribution_accessor(*eocit).push_back(std::make_pair(cit.handle(), interface_contribution) ); //Note: Due to iteration over cells there is no need to use a voronoi_unique_quantity_update() here 00866 interface_area_accessor(*eocit) += interface_contribution; 00867 00868 //box volume contribution: 00869 double edge_contribution = 0; 00870 VertexOnEdgeRange vertices_on_edge = viennagrid::elements<VertexType>(*eocit); 00871 for (VertexOnEdgeIterator voeit = vertices_on_edge.begin(); 00872 voeit != vertices_on_edge.end(); 00873 ++voeit) 00874 { 00875 //double contribution = spanned_volume(cell_center, facet_center, edge_midpoint, voeit->point()); 00876 double contribution = spanned_volume(cell_center, facet_center, edge_midpoint, viennagrid::point(mesh_obj, *voeit)); 00877 vertex_box_volume_cell_contribution_accessor(*voeit).push_back(std::make_pair( cit.handle(), contribution) ); 00878 vertex_box_volume_accessor(*voeit) += contribution; 00879 edge_contribution += contribution; 00880 } 00881 edge_box_volume_cell_contribution_accessor(*eocit).push_back(std::make_pair( cit.handle(), edge_contribution) ); 00882 edge_box_volume_accessor(*eocit) += edge_contribution; 00883 } //for edges on facet 00884 00885 } //for facets on cell 00886 00887 } //for cells 00888 00889 } 00890 00891 } //namespace detail 00892 00893 // 00894 // The public interface 00895 // 00896 00908 template <typename ElementTypeOrTagT, 00909 typename MeshT, 00910 typename InterfaceAreaAccessorT, 00911 typename InterfaceAreaCellContributionAccessorT, 00912 typename VertexBoxVolumeAccessorT, 00913 typename VertexBoxVolumeCellContributionAccessorT, 00914 typename EdgeBoxVolumeAccessorT, 00915 typename EdgeBoxVolumeCellContributionAccessorT> 00916 void apply_voronoi(MeshT const & mesh_obj, 00917 InterfaceAreaAccessorT interface_area_accessor, 00918 InterfaceAreaCellContributionAccessorT interface_area_cell_contribution_accessor, 00919 VertexBoxVolumeAccessorT vertex_box_volume_accessor, 00920 VertexBoxVolumeCellContributionAccessorT vertex_box_volume_cell_contribution_accessor, 00921 EdgeBoxVolumeAccessorT edge_box_volume_accessor, 00922 EdgeBoxVolumeCellContributionAccessorT edge_box_volume_cell_contribution_accessor) 00923 { 00924 typedef typename viennagrid::result_of::element_tag<ElementTypeOrTagT>::type ElementTag; 00925 00926 detail::write_voronoi_info<ElementTag>(mesh_obj, 00927 interface_area_accessor, 00928 interface_area_cell_contribution_accessor, 00929 vertex_box_volume_accessor, 00930 vertex_box_volume_cell_contribution_accessor, 00931 edge_box_volume_accessor, 00932 edge_box_volume_cell_contribution_accessor, 00933 ElementTag()); 00934 } 00935 00936 } //namespace viennagrid 00937 #endif
1.7.6.1