|
ViennaGrid - The Vienna Grid Library
2.1.0
|
00001 #ifndef VIENNAGRID_ALGORITHM_CLOSEST_POINTS_HPP 00002 #define VIENNAGRID_ALGORITHM_CLOSEST_POINTS_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 #include <limits> 00021 #include <utility> 00022 00023 #include "viennagrid/forwards.hpp" 00024 #include "viennagrid/topology/all.hpp" 00025 #include "viennagrid/algorithm/norm.hpp" 00026 #include "viennagrid/algorithm/inner_prod.hpp" 00027 #include "viennagrid/algorithm/boundary.hpp" 00028 #include "viennagrid/mesh/mesh.hpp" 00029 00035 namespace viennagrid 00036 { 00037 namespace detail 00038 { 00039 00040 template <typename CoordT, typename CoordinateSystemT> 00041 std::pair<spatial_point<CoordT, CoordinateSystemT>, 00042 spatial_point<CoordT, CoordinateSystemT> > const & 00043 point_pair_with_shortest_distance( std::pair<spatial_point<CoordT, CoordinateSystemT>, 00044 spatial_point<CoordT, CoordinateSystemT> > const & pair_1, 00045 std::pair<spatial_point<CoordT, CoordinateSystemT>, 00046 spatial_point<CoordT, CoordinateSystemT> > const & pair_2) 00047 { 00048 if ( norm_2(pair_1.first - pair_1.second) <= norm_2(pair_2.first - pair_2.second) ) 00049 return pair_1; 00050 00051 return pair_2; 00052 } 00053 00054 template <typename CoordT, typename CoordinateSystemT> 00055 std::pair<spatial_point<CoordT, CoordinateSystemT>, 00056 spatial_point<CoordT, CoordinateSystemT> > const & 00057 point_pair_with_shortest_distance( std::pair<spatial_point<CoordT, CoordinateSystemT>, 00058 spatial_point<CoordT, CoordinateSystemT> > const & pair_1, 00059 std::pair<spatial_point<CoordT, CoordinateSystemT>, 00060 spatial_point<CoordT, CoordinateSystemT> > const & pair_2, 00061 std::pair<spatial_point<CoordT, CoordinateSystemT>, 00062 spatial_point<CoordT, CoordinateSystemT> > const & pair_3) 00063 { 00064 CoordT dist_pair_1 = norm_2(pair_1.first - pair_1.second); 00065 CoordT dist_pair_2 = norm_2(pair_2.first - pair_2.second); 00066 CoordT dist_pair_3 = norm_2(pair_3.first - pair_3.second); 00067 00068 if ( (dist_pair_1 <= dist_pair_2) 00069 && (dist_pair_1 <= dist_pair_3) ) 00070 return pair_1; 00071 00072 if ( (dist_pair_2 <= dist_pair_1) 00073 && (dist_pair_2 <= dist_pair_3) ) 00074 return pair_2; 00075 00076 return pair_3; 00077 } 00078 00079 00080 template <typename CoordT, typename CoordinateSystemT> 00081 std::pair<spatial_point<CoordT, CoordinateSystemT>, 00082 spatial_point<CoordT, CoordinateSystemT> > const & 00083 point_pair_with_shortest_distance( std::pair<spatial_point<CoordT, CoordinateSystemT>, 00084 spatial_point<CoordT, CoordinateSystemT> > const & pair_1, 00085 std::pair<spatial_point<CoordT, CoordinateSystemT>, 00086 spatial_point<CoordT, CoordinateSystemT> > const & pair_2, 00087 std::pair<spatial_point<CoordT, CoordinateSystemT>, 00088 spatial_point<CoordT, CoordinateSystemT> > const & pair_3, 00089 std::pair<spatial_point<CoordT, CoordinateSystemT>, 00090 spatial_point<CoordT, CoordinateSystemT> > const & pair_4) 00091 { 00092 // Note: Recursive use of point_pair_with_shortest_distance() has a bit of overhead. Feel free to improve this 00093 return point_pair_with_shortest_distance(point_pair_with_shortest_distance(pair_1, pair_2), 00094 point_pair_with_shortest_distance(pair_3, pair_4)); 00095 } 00096 00097 template <typename CoordT, typename CoordinateSystemT> 00098 std::pair<spatial_point<CoordT, CoordinateSystemT>, 00099 spatial_point<CoordT, CoordinateSystemT> > const & 00100 point_pair_with_shortest_distance( std::pair<spatial_point<CoordT, CoordinateSystemT>, 00101 spatial_point<CoordT, CoordinateSystemT> > const & pair_1, 00102 std::pair<spatial_point<CoordT, CoordinateSystemT>, 00103 spatial_point<CoordT, CoordinateSystemT> > const & pair_2, 00104 std::pair<spatial_point<CoordT, CoordinateSystemT>, 00105 spatial_point<CoordT, CoordinateSystemT> > const & pair_3, 00106 std::pair<spatial_point<CoordT, CoordinateSystemT>, 00107 spatial_point<CoordT, CoordinateSystemT> > const & pair_4, 00108 std::pair<spatial_point<CoordT, CoordinateSystemT>, 00109 spatial_point<CoordT, CoordinateSystemT> > const & pair_5) 00110 { 00111 // Note: Recursive use of point_pair_with_shortest_distance() has a bit of overhead. Feel free to improve this 00112 return point_pair_with_shortest_distance(point_pair_with_shortest_distance(pair_1, pair_2, pair_3), 00113 point_pair_with_shortest_distance(pair_4, pair_5)); 00114 } 00115 00116 template <typename CoordT, typename CoordinateSystemT> 00117 std::pair<spatial_point<CoordT, CoordinateSystemT>, 00118 spatial_point<CoordT, CoordinateSystemT> > const & 00119 point_pair_with_shortest_distance( std::pair<spatial_point<CoordT, CoordinateSystemT>, 00120 spatial_point<CoordT, CoordinateSystemT> > const & pair_1, 00121 std::pair<spatial_point<CoordT, CoordinateSystemT>, 00122 spatial_point<CoordT, CoordinateSystemT> > const & pair_2, 00123 std::pair<spatial_point<CoordT, CoordinateSystemT>, 00124 spatial_point<CoordT, CoordinateSystemT> > const & pair_3, 00125 std::pair<spatial_point<CoordT, CoordinateSystemT>, 00126 spatial_point<CoordT, CoordinateSystemT> > const & pair_4, 00127 std::pair<spatial_point<CoordT, CoordinateSystemT>, 00128 spatial_point<CoordT, CoordinateSystemT> > const & pair_5, 00129 std::pair<spatial_point<CoordT, CoordinateSystemT>, 00130 spatial_point<CoordT, CoordinateSystemT> > const & pair_6) 00131 { 00132 // Note: Recursive use of point_pair_with_shortest_distance() has a bit of overhead. Feel free to improve this 00133 return point_pair_with_shortest_distance(point_pair_with_shortest_distance(pair_1, pair_2, pair_3), 00134 point_pair_with_shortest_distance(pair_4, pair_5, pair_6)); 00135 } 00136 00137 00138 // 00139 // Closest points between points (trivial) 00140 // 00141 00142 // Closest points between two points: 00143 template <typename PointT> 00144 std::pair<PointT, PointT> 00145 closest_points_impl(PointT const & p1, 00146 PointT const & p2) 00147 { 00148 return std::make_pair(p1, p2); 00149 } 00150 00151 template <typename PointAccessorT, typename CoordT1, typename CoordinateSystemT1, typename CoordT2, typename CoordinateSystemT2> 00152 std::pair< spatial_point<CoordT1, CoordinateSystemT1>, spatial_point<CoordT2, CoordinateSystemT2> > 00153 closest_points_impl(PointAccessorT const, 00154 spatial_point<CoordT1, CoordinateSystemT1> const & p1, 00155 spatial_point<CoordT2, CoordinateSystemT2> const & p2) 00156 { 00157 return closest_points_impl(p1, p2); 00158 } 00159 00160 template <typename PointAccessorT, typename PointT, typename WrappedConfigT> 00161 std::pair<PointT, typename viennagrid::result_of::point<viennagrid::element< vertex_tag,WrappedConfigT> >::type> 00162 closest_points_impl(PointAccessorT const accessor, 00163 PointT const & p1, 00164 viennagrid::element<vertex_tag,WrappedConfigT> const & v2) 00165 { 00166 return closest_points_impl(p1, accessor(v2)); 00167 } 00168 00169 template <typename PointAccessorT, typename PointT, typename WrappedConfigT> 00170 std::pair<typename viennagrid::result_of::point<viennagrid::element< vertex_tag,WrappedConfigT> >::type, PointT> 00171 closest_points_impl(PointAccessorT const accessor, 00172 viennagrid::element<vertex_tag,WrappedConfigT> const & v1, 00173 PointT const & p2) 00174 { 00175 return closest_points_impl(accessor(v1), p2); 00176 } 00177 00178 // Closest points between vertices: 00179 template <typename PointAccessorT, typename PointT, typename WrappedConfigT> 00180 std::pair<typename viennagrid::result_of::point<viennagrid::element< vertex_tag,WrappedConfigT> >::type, typename viennagrid::result_of::point<viennagrid::element< vertex_tag,WrappedConfigT> >::type> 00181 closest_points_impl(PointAccessorT const accessor, 00182 viennagrid::element<vertex_tag,WrappedConfigT> const & v1, 00183 viennagrid::element<vertex_tag,WrappedConfigT> const & v2) 00184 { 00185 return closest_points_impl( accessor(v1), accessor(v2) ); 00186 } 00187 00188 00189 00190 // 00191 // Closest points between point and line segment 00192 // 00193 00194 // Implementation: Supposed to work for arbitrary dimensions 00195 template <typename PointT, typename LinePointT> 00196 std::pair<PointT, LinePointT> 00197 closest_points_point_line(PointT const & p, 00198 LinePointT const & line_p1, 00199 LinePointT const & line_p2) 00200 { 00201 //typedef spatial_point<CoordT, CoordinateSystemT> PointT; 00202 typedef typename viennagrid::result_of::coord< PointT >::type coord_type; 00203 00204 //compute t such that projection of p onto [line_p1, line_p2] is given by p' = line_p1 + t * (line_p2 - line_p1) 00205 coord_type t = viennagrid::inner_prod( (p - line_p1), (line_p2 - line_p1) ) / viennagrid::inner_prod(line_p2 - line_p1, line_p2 - line_p1); 00206 00207 //restrict t to line segment, i.e. t \in [0, 1] 00208 t = std::max<coord_type>(0, std::min<coord_type>(1, t)); 00209 00210 LinePointT p_prime = line_p1 + t * (line_p2 - line_p1); //closest point to p on line 00211 00212 return std::make_pair(p, p_prime); 00213 } 00214 00215 00216 //convenience overload: point and simplex line 00217 template <typename PointAccessorT, typename PointT, typename WrappedConfigT> 00218 std::pair< PointT, typename PointAccessorT::value_type > 00219 closest_points_impl(PointAccessorT const accessor, 00220 PointT const & p, 00221 viennagrid::element<simplex_tag<1>,WrappedConfigT> const & el) 00222 { 00223 return closest_points_point_line(p, 00224 accessor(vertices(el)[0]), 00225 accessor(vertices(el)[1])); 00226 } 00227 00228 //convenience overload: point and hypercube line 00229 template <typename PointAccessorT, typename PointT, typename WrappedConfigT> 00230 std::pair< PointT, typename PointAccessorT::value_type > 00231 closest_points_impl(PointAccessorT const accessor, 00232 PointT const & p, 00233 viennagrid::element<hypercube_tag<1>,WrappedConfigT> const & el) 00234 { 00235 return closest_points_point_line(p, 00236 accessor(vertices(el)[0]), 00237 accessor(vertices(el)[1])); 00238 } 00239 00240 //convenience overload: vertex and simplex line 00241 template <typename PointAccessorT, typename WrappedConfigT1, typename WrappedConfigT2> 00242 std::pair< typename PointAccessorT::value_type, typename PointAccessorT::value_type > 00243 closest_points_impl(PointAccessorT const accessor, 00244 viennagrid::element<vertex_tag,WrappedConfigT1> const & v, 00245 viennagrid::element<simplex_tag<1>,WrappedConfigT2> const & el) 00246 { 00247 return closest_points_point_line(accessor(v), 00248 accessor(vertices(el)[0]), 00249 accessor(vertices(el)[1])); 00250 } 00251 00252 //convenience overload: vertex and hypercube line 00253 template <typename PointAccessorT, typename PointT, typename WrappedConfigT1, typename WrappedConfigT2> 00254 std::pair< typename PointAccessorT::value_type, typename PointAccessorT::value_type > 00255 closest_points_impl(PointAccessorT const accessor, 00256 viennagrid::element<vertex_tag,WrappedConfigT1> const & v, 00257 viennagrid::element<simplex_tag<1>,WrappedConfigT2> const & el) 00258 { 00259 return closest_points_point_line(accessor(v), 00260 accessor(vertices(el)[0]), 00261 accessor(vertices(el)[1])); 00262 } 00263 00264 // 00265 // Distance between two line segments 00266 // 00267 00268 // Implementation: Supposed to work for arbitrary dimensions 00269 template <typename PointT> 00270 std::pair<PointT, PointT> 00271 closest_points_line_line(PointT const & v0, PointT const & v1, //endpoints of first line 00272 PointT const & w0, PointT const & w1) //endpoints of second line 00273 { 00274 //typedef spatial_point<CoordT, CoordinateSystemT> PointT; 00275 typedef typename viennagrid::result_of::coord< PointT >::type coord_type; 00276 00277 // write V(s) = v0 + s * (v1 - v0), s \in [0,1] 00278 // W(t) = w0 + t * (w1 - w0), t \in [0,1] 00279 00280 // compute s and t assuming V(s) and W(t) to be infinite lines: 00281 // cf. http://www.softsurfer.com/Archive/algorithm_0106/algorithm_0106.htm 00282 PointT dir_v = v1 - v0; //direction vector for line V(s) 00283 PointT dir_w = w1 - w0; //direction vector for line W(t) 00284 00285 coord_type v_in_v = viennagrid::inner_prod(dir_v, dir_v); 00286 coord_type v_in_w = viennagrid::inner_prod(dir_v, dir_w); 00287 coord_type w_in_w = viennagrid::inner_prod(dir_w, dir_w); 00288 00289 coord_type denominator = v_in_v * w_in_w - v_in_w * v_in_w; 00290 //std::cout << "denominator: " << denominator << std::endl; 00291 00292 if (denominator < 1e-6 * v_in_v * w_in_w) //lines parallel (up to round-off) 00293 { 00294 return point_pair_with_shortest_distance(closest_points_point_line(w0, v0, v1), 00295 closest_points_point_line(w1, v0, v1), 00296 closest_points_point_line(v0, w0, w1), 00297 closest_points_point_line(v1, w0, w1)); 00298 } 00299 00300 //Lines are not parallel: Compute minimizers s, t: 00301 PointT dir_distance = v0 - w0; //any vector connecting two points on V and W 00302 00303 //if (inner_prod(dir_distance, dir_distance) / v_in_v < 1e-10) //v0 and w0 are the same point 00304 // return std::make_pair(v0, w0); 00305 00306 coord_type v_in_dir_distance = viennagrid::inner_prod(dir_v, dir_distance); 00307 coord_type w_in_dir_distance = viennagrid::inner_prod(dir_w, dir_distance); 00308 00309 coord_type s = (v_in_w * w_in_dir_distance - w_in_w * v_in_dir_distance) / denominator; 00310 coord_type t = (v_in_v * w_in_dir_distance - v_in_w * v_in_dir_distance) / denominator; 00311 //std::cout << "s = " << s << std::endl; 00312 //std::cout << "t = " << s << std::endl; 00313 00314 //Check whether minimizing distance is within line segment, i.e. s,t \in [0,1] 00315 if ( (s < 0) || (s > 1) || (t < 0) || (t > 1) ) //Note: A more fine-grained check is possible for those looking for optimizations 00316 return point_pair_with_shortest_distance(closest_points_point_line(w0, v0, v1), 00317 closest_points_point_line(w1, v0, v1), 00318 closest_points_point_line(v0, w0, w1), 00319 closest_points_point_line(v1, w0, w1)); 00320 00321 // compute points on V(s) an W(t) for which distance is smallest: 00322 PointT min_dist_point_V = v0 + s * dir_v; 00323 PointT min_dist_point_W = w0 + t * dir_w; 00324 00325 return std::make_pair(min_dist_point_V, min_dist_point_W); 00326 } 00327 00328 00329 //convenience overload: simplex line and simplex line 00330 template <typename PointAccessorT, typename WrappedConfigT1, typename WrappedConfigT2> 00331 std::pair<typename PointAccessorT::value_type, typename PointAccessorT::value_type> 00332 closest_points_impl(PointAccessorT const accessor, 00333 viennagrid::element<simplex_tag<1>,WrappedConfigT1> const & line0, 00334 viennagrid::element<simplex_tag<1>,WrappedConfigT2> const & line1) 00335 { 00336 return closest_points_line_line(accessor(vertices(line0)[0]), 00337 accessor(vertices(line0)[1]), 00338 accessor(vertices(line1)[0]), 00339 accessor(vertices(line1)[1])); 00340 } 00341 00342 //convenience overload: hypercube line and simplex line 00343 template <typename PointAccessorT, typename WrappedConfigT1, typename WrappedConfigT2> 00344 std::pair<typename PointAccessorT::value_type, typename PointAccessorT::value_type> 00345 closest_points_impl(PointAccessorT const accessor, 00346 viennagrid::element<hypercube_tag<1>,WrappedConfigT1> const & line0, 00347 viennagrid::element<simplex_tag<1>,WrappedConfigT2> const & line1) 00348 { 00349 return closest_points_line_line(accessor(vertices(line0)[0]), 00350 accessor(vertices(line0)[1]), 00351 accessor(vertices(line1)[0]), 00352 accessor(vertices(line1)[1])); 00353 } 00354 00355 //convenience overload: simplex line and hypercube line 00356 template <typename PointAccessorT, typename WrappedConfigT1, typename WrappedConfigT2> 00357 std::pair<typename PointAccessorT::value_type, typename PointAccessorT::value_type> 00358 closest_points_impl(PointAccessorT const accessor, 00359 viennagrid::element<simplex_tag<1>,WrappedConfigT1> const & line0, 00360 viennagrid::element<hypercube_tag<1>,WrappedConfigT2> const & line1) 00361 { 00362 return closest_points_line_line(accessor(vertices(line0)[0]), 00363 accessor(vertices(line0)[1]), 00364 accessor(vertices(line1)[0]), 00365 accessor(vertices(line1)[1])); 00366 } 00367 00368 //convenience overload: hypercube line and hypercube line 00369 template <typename PointAccessorT, typename WrappedConfigT1, typename WrappedConfigT2> 00370 std::pair<typename PointAccessorT::value_type, typename PointAccessorT::value_type> 00371 closest_points_impl(PointAccessorT const accessor, 00372 viennagrid::element<hypercube_tag<1>,WrappedConfigT1> const & line0, 00373 viennagrid::element<hypercube_tag<1>,WrappedConfigT2> const & line1) 00374 { 00375 return closest_points_line_line(accessor(vertices(line0)[0]), 00376 accessor(vertices(line0)[1]), 00377 accessor(vertices(line1)[0]), 00378 accessor(vertices(line1)[1])); 00379 } 00380 00381 00382 // 00383 // Distance between point and triangle 00384 // 00385 00386 // Implementation: Supposed to work for arbitrary dimensions 00387 // Projects p onto the plane spanned by the triangle, then computes the shortest distance in the plane and uses Pythagoras for the full distance 00388 template <typename PointT> 00389 std::pair<PointT,PointT> 00390 closest_points_point_triangle(PointT const & p, 00391 PointT const & v0, 00392 PointT const & v1, 00393 PointT const & v2) //endpoints of second line 00394 { 00395 typedef typename viennagrid::result_of::coord< PointT >::type coord_type; 00396 00397 // write T(s) = v0 + s * (v1 - v0) + t * (v2 - v0), {s,t} \in [0,1], s+t < 1 for the triangle T 00398 // =: v0 + s * u0 + t * u1 00399 // 00400 // Projection p' of p is given by the solution of the system 00401 // 00402 // <u0, u0> * s + <u1, u0> * t = <u, u0> 00403 // <u0, u1> * s + <u1, u1> * t = <u, u1>, 00404 // 00405 // where u = p - v0. This is a 2x2-system that is directly inverted. 00406 // 00407 // Write 00408 // 00409 // A = ( <u0, u0> <u1, u0> ) = ( a b ) 00410 // ( <u0, u1> <u1, u1> ) ( c d ) 00411 // 00412 // and use b = c because of symmetry of the inner product. 00413 00414 PointT u0 = v1 - v0; 00415 PointT u1 = v2 - v0; 00416 PointT u = p - v0; 00417 00418 coord_type a = viennagrid::inner_prod(u0, u0); 00419 coord_type b = viennagrid::inner_prod(u0, u1); 00420 coord_type d = viennagrid::inner_prod(u1, u1); 00421 coord_type u_in_u0 = viennagrid::inner_prod(u, u0); 00422 coord_type u_in_u1 = viennagrid::inner_prod(u, u1); 00423 00424 coord_type denominator = a * d - b * b; 00425 00426 if (denominator < 1e-6 * a * d) //triangle is VERY thin: TODO: NUMERIC!!! 1e-6 00427 { 00428 std::cerr << "ViennaGrid: Warning: Strongly degenerated triangle detected: " << std::endl 00429 << "vertex 0: " << v0 << std::endl 00430 << "vertex 1: " << v1 << std::endl 00431 << "vertex 2: " << v2 << std::endl; 00432 } 00433 00434 coord_type s = ( d * u_in_u0 - b * u_in_u1) / denominator; 00435 coord_type t = (- b * u_in_u0 + a * u_in_u1) / denominator; 00436 00437 PointT p_prime = v0 + s * u0 + t * u1; //projection of p onto triangular plane 00438 00439 // nonzero distance is encountered only if p_prime is outside the triangle 00440 if ( (s < 0 || s > 1) 00441 || (t < 0 || t > 1) 00442 || (s + t > 1) ) //p_prime is outside the triangle 00443 { 00444 // safe mode: Compute distances to all edges. Can be optimized further by using information from s, t, etc. 00445 return point_pair_with_shortest_distance(std::make_pair(p, closest_points_point_line(p_prime, v0, v1).second), 00446 std::make_pair(p, closest_points_point_line(p_prime, v0, v2).second), 00447 std::make_pair(p, closest_points_point_line(p_prime, v1, v2).second)); 00448 } 00449 00450 return std::make_pair(p, p_prime); 00451 } 00452 00453 00454 //convenience overload: point 00455 template <typename PointAccessorT, typename CoordT, typename CoordinateSystemT, typename WrappedConfigT> 00456 std::pair<spatial_point<CoordT, CoordinateSystemT>, typename PointAccessorT::value_type> 00457 closest_points_impl(PointAccessorT const accessor, 00458 spatial_point<CoordT, CoordinateSystemT> const & p, 00459 viennagrid::element<simplex_tag<2>,WrappedConfigT> const & el) 00460 { 00461 return closest_points_point_triangle(p, 00462 accessor(vertices(el)[0]), 00463 accessor(vertices(el)[1]), 00464 accessor(vertices(el)[2])); 00465 } 00466 00467 //convenience overload: vertex 00468 template <typename PointAccessorT, typename WrappedConfigT1, typename WrappedConfigT2> 00469 std::pair<typename PointAccessorT::value_type,typename PointAccessorT::value_type> 00470 closest_points_impl(PointAccessorT const accessor, 00471 viennagrid::element<vertex_tag,WrappedConfigT1> const & v, 00472 viennagrid::element<simplex_tag<2>,WrappedConfigT2> const & el) 00473 { 00474 return closest_points_impl( accessor(v), el ); 00475 } 00476 00477 00478 00479 00480 00481 // 00482 // Distance between point and quadrilateral (using decomposition into two triangles) 00483 // 00484 00485 //convenience overload: point 00486 // 00487 // Keep reference orientation in mind: 00488 // 00489 // v[2] v[3] 00490 // * --------- * 00491 // | \ | 00492 // | \ | 00493 // * --------- * 00494 // v[0] v[1] 00495 // 00496 template <typename PointAccessorT, typename CoordT, typename CoordinateSystemT, typename WrappedConfigT> 00497 std::pair<spatial_point<CoordT, CoordinateSystemT>, typename PointAccessorT::value_type> 00498 closest_points_impl(PointAccessorT const accessor, 00499 spatial_point<CoordT, CoordinateSystemT> const & p, 00500 viennagrid::element<hypercube_tag<2>,WrappedConfigT> const & el) 00501 { 00502 return point_pair_with_shortest_distance(closest_points_point_triangle(p, 00503 accessor(vertices(el)[0]), 00504 accessor(vertices(el)[1]), 00505 accessor(vertices(el)[2])), 00506 closest_points_point_triangle(p, 00507 accessor(vertices(el)[1]), 00508 accessor(vertices(el)[3]), 00509 accessor(vertices(el)[2])) 00510 ); 00511 } 00512 00513 //convenience overload: vertex 00514 template <typename PointAccessorT, typename WrappedConfigT1, typename WrappedConfigT2> 00515 std::pair<typename PointAccessorT::value_type,typename PointAccessorT::value_type> 00516 closest_points_impl(PointAccessorT const accessor, 00517 viennagrid::element<vertex_tag,WrappedConfigT1> const & v, 00518 viennagrid::element<hypercube_tag<2>,WrappedConfigT2> const & el) 00519 { 00520 return closest_points_impl( accessor(v), el); 00521 } 00522 00523 00524 00525 // 00526 // Distance between point and tetrahedron 00527 // 00528 00529 // Implementation: Supposed to work for arbitrary geometric dimensions 00530 // Projects p onto the plane spanned by the triangle, then computes the shortest distance in the plane and uses Pythagoras for the full distance 00531 template <typename PointT> 00532 std::pair<PointT,PointT> 00533 closest_points_point_tetrahedron(PointT const & p, 00534 PointT const & v0, 00535 PointT const & v1, 00536 PointT const & v2, 00537 PointT const & v3) 00538 { 00539 typedef typename viennagrid::result_of::coord< PointT >::type coord_type; 00540 00541 // write T(s) = v0 + r * (v1 - v0) + s * (v2 - v0) + t * (v3 - v1), {r,s,t} \in [0,1], r+s+t < 1 for the tetrahedron T 00542 // =: v0 + r * u0 + s * u1 + t * u2 00543 // 00544 // Projection p' of p is given by the solution of the system 00545 // 00546 // <u0, u0> * r + <u1, u0> * s + <u2, u0> * t = <u, u0> 00547 // <u0, u1> * r + <u1, u1> * s + <u2, u1> * t = <u, u1> 00548 // <u0, u2> * r + <u1, u2> * s + <u2, u2> * t = <u, u2> 00549 // 00550 // where u = p - v0. This is a 3x3-system that is directly inverted using Cramer's rule. 00551 // 00552 00553 PointT u0 = v1 - v0; 00554 PointT u1 = v2 - v0; 00555 PointT u2 = v3 - v0; 00556 PointT u = p - v0; 00557 00558 coord_type a_00 = viennagrid::inner_prod(u0, u0); 00559 coord_type a_01 = viennagrid::inner_prod(u1, u0); 00560 coord_type a_02 = viennagrid::inner_prod(u2, u0); 00561 00562 coord_type a_10 = a_01; 00563 coord_type a_11 = viennagrid::inner_prod(u1, u1); 00564 coord_type a_12 = viennagrid::inner_prod(u2, u1); 00565 00566 coord_type a_20 = a_02; 00567 coord_type a_21 = a_12; 00568 coord_type a_22 = viennagrid::inner_prod(u2, u2); 00569 00570 coord_type u_in_u0 = viennagrid::inner_prod(u, u0); 00571 coord_type u_in_u1 = viennagrid::inner_prod(u, u1); 00572 coord_type u_in_u2 = viennagrid::inner_prod(u, u2); 00573 00574 coord_type det_A = a_00 * a_11 * a_22 + a_01 * a_12 * a_20 + a_02 * a_10 * a_21 00575 - a_20 * a_11 * a_02 - a_21 * a_12 * a_00 - a_22 * a_10 * a_01; 00576 00577 if (det_A < 1e-6 * std::sqrt(a_00 * a_11 * a_22)) //tetrahedron is VERY thin: 00578 { 00579 std::cerr << "ViennaGrid: Warning: Strongly degenerated tetrahedron detected: " << std::endl 00580 << "vertex 0: " << v0 << std::endl 00581 << "vertex 1: " << v1 << std::endl 00582 << "vertex 2: " << v2 << std::endl 00583 << "vertex 3: " << v3 << std::endl; 00584 } 00585 00586 // | <u, u0> a_01 a_02 | 00587 // r = | <u, u1> a_11 a_12 | / det_A and similarly for s and t 00588 // | <u, u2> a_21 a_22 | 00589 coord_type r = ( u_in_u0 * a_11 * a_22 + a_01 * a_12 * u_in_u2 + a_02 * u_in_u1 * a_21 00590 - u_in_u2 * a_11 * a_02 - a_21 * a_12 * u_in_u0 - a_22 * u_in_u1 * a_01 ) / det_A; 00591 coord_type s = ( a_00 * u_in_u1 * a_22 + u_in_u0 * a_12 * a_20 + a_02 * a_10 * u_in_u2 00592 - a_20 * u_in_u1 * a_02 - u_in_u2 * a_12 * a_00 - a_22 * a_10 * u_in_u0 ) / det_A; 00593 coord_type t = ( a_00 * a_11 * u_in_u2 + a_01 * u_in_u1 * a_20 + u_in_u0 * a_10 * a_21 00594 - a_20 * a_11 * u_in_u0 - a_21 * u_in_u1 * a_00 - u_in_u2 * a_10 * a_01 ) / det_A; 00595 00596 PointT p_prime = v0 + r * u0 + s * u1 + t * u2; //projection of p onto triangular plane 00597 00598 // nonzero distance is encountered only if p_prime is outside the triangle 00599 if ( (r < 0 || r > 1) 00600 || (s < 0 || s > 1) 00601 || (t < 0 || t > 1) 00602 || (r + s + t > 1) ) //p_prime is outside the triangle 00603 { 00604 // safe mode: Compute distances to all edges. Can be optimized further by using information from s, t, etc. 00605 return point_pair_with_shortest_distance(std::make_pair(p, closest_points_point_triangle(p_prime, v0, v1, v2).second), 00606 std::make_pair(p, closest_points_point_triangle(p_prime, v0, v1, v3).second), 00607 std::make_pair(p, closest_points_point_triangle(p_prime, v0, v2, v3).second), 00608 std::make_pair(p, closest_points_point_triangle(p_prime, v1, v2, v3).second)); 00609 } 00610 00611 return std::make_pair(p, p_prime); 00612 } 00613 00614 00615 //convenience overload: point 00616 template <typename PointAccessorT, typename CoordT, typename CoordinateSystemT, typename WrappedConfigT> 00617 std::pair<spatial_point<CoordT, CoordinateSystemT>,typename PointAccessorT::value_type> 00618 closest_points_impl(PointAccessorT const accessor, 00619 spatial_point<CoordT, CoordinateSystemT> const & p, 00620 viennagrid::element<simplex_tag<3>,WrappedConfigT> const & el) 00621 { 00622 return closest_points_point_tetrahedron(p, 00623 accessor(vertices(el)[0]), 00624 accessor(vertices(el)[1]), 00625 accessor(vertices(el)[2]), 00626 accessor(vertices(el)[3])); 00627 } 00628 00629 //convenience overload: vertex 00630 template <typename PointAccessorT, typename WrappedConfigT1, typename WrappedConfigT2> 00631 std::pair<typename PointAccessorT::value_type,typename PointAccessorT::value_type> 00632 closest_points_impl(PointAccessorT const accessor, 00633 viennagrid::element<vertex_tag,WrappedConfigT1> const & v, 00634 viennagrid::element<simplex_tag<3>,WrappedConfigT2> const & el) 00635 { 00636 return closest_points_impl( accessor(v), el); 00637 } 00638 00639 00640 00641 // 00642 // Distance between point and hexahedron (using decomposition into six tetrahedra) 00643 // 00644 00645 //convenience overload: point 00646 template <typename PointAccessorT, typename PointT, typename WrappedConfigT> 00647 std::pair<PointT,typename PointAccessorT::value_type> 00648 closest_points_impl(PointAccessorT const accessor, 00649 PointT const & p, 00650 viennagrid::element<hypercube_tag<3>,WrappedConfigT> const & el) 00651 { 00652 return point_pair_with_shortest_distance(closest_points_point_tetrahedron(p, 00653 accessor(vertices(el)[0]), 00654 accessor(vertices(el)[1]), 00655 accessor(vertices(el)[3]), 00656 accessor(vertices(el)[7])), 00657 closest_points_point_tetrahedron(p, 00658 accessor(vertices(el)[0]), 00659 accessor(vertices(el)[1]), 00660 accessor(vertices(el)[7]), 00661 accessor(vertices(el)[5])), 00662 closest_points_point_tetrahedron(p, 00663 accessor(vertices(el)[0]), 00664 accessor(vertices(el)[2]), 00665 accessor(vertices(el)[6]), 00666 accessor(vertices(el)[7])), 00667 closest_points_point_tetrahedron(p, 00668 accessor(vertices(el)[0]), 00669 accessor(vertices(el)[3]), 00670 accessor(vertices(el)[2]), 00671 accessor(vertices(el)[7])), 00672 closest_points_point_tetrahedron(p, 00673 accessor(vertices(el)[0]), 00674 accessor(vertices(el)[5]), 00675 accessor(vertices(el)[7]), 00676 accessor(vertices(el)[4])), 00677 closest_points_point_tetrahedron(p, 00678 accessor(vertices(el)[0]), 00679 accessor(vertices(el)[6]), 00680 accessor(vertices(el)[4]), 00681 accessor(vertices(el)[7])) 00682 ); 00683 } 00684 00685 //convenience overload: vertex 00686 template <typename PointAccessorT, typename WrappedConfigT1, typename WrappedConfigT2> 00687 std::pair<typename PointAccessorT::value_type,typename PointAccessorT::value_type> 00688 closest_points_impl(PointAccessorT const accessor, 00689 viennagrid::element<vertex_tag,WrappedConfigT1> const & v, 00690 viennagrid::element<hypercube_tag<3>,WrappedConfigT2> const & el) 00691 { 00692 return closest_points_impl( accessor(v), el); 00693 } 00694 00695 // 00697 // 00698 00699 // 00700 // Closest points between points (trivial) 00701 // 00702 00703 // Closest points between two points (overloads for vertices follow): 00704 template <typename PointT> 00705 std::pair<PointT, PointT> 00706 closest_points_on_boundary_impl(PointT const & p1, 00707 PointT const & p2) 00708 { 00709 return std::make_pair(p1, p2); 00710 } 00711 00712 template <typename PointAccessorT, typename CoordT1, typename CoordinateSystemT1, typename CoordT2, typename CoordinateSystemT2> 00713 std::pair< spatial_point<CoordT1, CoordinateSystemT1>, spatial_point<CoordT1, CoordinateSystemT1> > 00714 closest_points_on_boundary_impl(PointAccessorT const, 00715 spatial_point<CoordT1, CoordinateSystemT1> const & p1, 00716 spatial_point<CoordT2, CoordinateSystemT2> const & p2) 00717 { 00718 return closest_points_on_boundary_impl(p1, p2); 00719 } 00720 00721 template <typename PointAccessorT, typename PointT, typename WrappedConfigT> 00722 std::pair<PointT, typename viennagrid::result_of::point<viennagrid::element< vertex_tag,WrappedConfigT> >::type> 00723 closest_points_on_boundary_impl(PointAccessorT const accessor, 00724 PointT const & p1, 00725 viennagrid::element<vertex_tag,WrappedConfigT> const & v2) 00726 { 00727 return std::make_pair(p1, accessor(v2)); 00728 } 00729 00730 template <typename PointAccessorT, typename PointT, typename WrappedConfigT> 00731 std::pair<typename viennagrid::result_of::point<viennagrid::element< vertex_tag,WrappedConfigT> >::type, PointT> 00732 closest_points_on_boundary_impl(PointAccessorT const accessor, 00733 PointT const & p2, 00734 viennagrid::element<vertex_tag,WrappedConfigT> const & v1) 00735 { 00736 return std::make_pair(accessor(v1), p2); 00737 } 00738 00739 // Closest points between vertices: 00740 template <typename PointAccessorT, typename PointT, typename WrappedConfigT> 00741 std::pair<typename viennagrid::result_of::point<viennagrid::element< vertex_tag,WrappedConfigT> >::type, typename viennagrid::result_of::point<viennagrid::element< vertex_tag,WrappedConfigT> >::type> 00742 closest_points_on_boundary_impl(PointAccessorT const accessor, 00743 viennagrid::element<vertex_tag,WrappedConfigT> const & v1, 00744 viennagrid::element<vertex_tag,WrappedConfigT> const & v2) 00745 { 00746 return std::make_pair( accessor(v1), accessor(v2) ); 00747 } 00748 00750 00752 template <typename FacetTypeOrTag, typename PointAccessorT, typename PointT, typename SomethingT> 00753 std::pair<PointT, typename PointAccessorT::value_type> 00754 closest_points_on_boundary_point_to_any(PointAccessorT const accessor, 00755 PointT const & p, 00756 SomethingT const & cont) 00757 { 00758 typedef typename viennagrid::result_of::element_tag<FacetTypeOrTag>::type FacetTag; 00759 typedef std::pair<PointT, typename viennagrid::result_of::point<SomethingT>::type> PairType; 00760 00761 typedef typename viennagrid::result_of::const_element_range<SomethingT, FacetTag>::type FacetRange; 00762 typedef typename viennagrid::result_of::iterator<FacetRange>::type FacetIterator; 00763 00764 PairType closest_pair; 00765 double shortest_distance = std::numeric_limits<double>::max(); 00766 00767 FacetRange facets(cont); 00768 for (FacetIterator fit = facets.begin(); 00769 fit != facets.end(); 00770 ++fit) 00771 { 00772 if (!is_boundary(cont, *fit)) 00773 continue; 00774 00775 PairType pair = closest_points_impl(accessor, p, *fit); 00776 double cur_norm = norm_2(pair.first - pair.second); 00777 if (cur_norm < shortest_distance) 00778 { 00779 closest_pair = pair; 00780 shortest_distance = cur_norm; 00781 } 00782 } 00783 00784 return closest_pair; 00785 } 00786 00787 00788 00789 template <typename PointAccessorT, typename PointT, typename ElementTag, typename WrappedConfigT> 00790 std::pair<PointT, typename PointAccessorT::value_type> 00791 closest_points_on_boundary_impl(PointAccessorT const & mesh_obj, 00792 PointT const & v, 00793 viennagrid::element<ElementTag,WrappedConfigT> const & el) 00794 { 00795 return closest_points_on_boundary_point_to_any<typename ElementTag::facet_tag>(mesh_obj, v, el); 00796 } 00797 00798 template <typename PointAccessorT, typename WrappedConfigT, typename PointT> 00799 std::pair<PointT, typename viennagrid::result_of::point< viennagrid::mesh<WrappedConfigT> >::type> 00800 closest_points_on_boundary_impl(PointAccessorT const point_accessor, 00801 PointT const & p, 00802 mesh<WrappedConfigT> const & mesh_obj) 00803 { 00804 typedef typename viennagrid::result_of::cell_tag< mesh<WrappedConfigT> >::type CellTag; 00805 return closest_points_on_boundary_point_to_any<typename CellTag::facet_tag>( point_accessor, p, mesh_obj); 00806 } 00807 00808 template <typename PointAccessorT, typename SegmentationT, typename PointT> 00809 std::pair<PointT, typename viennagrid::result_of::point< segment_handle<SegmentationT> >::type> 00810 closest_points_on_boundary_impl(PointAccessorT const point_accessor, 00811 PointT const & p, 00812 segment_handle<SegmentationT> const & segment) 00813 { 00814 typedef typename viennagrid::result_of::cell_tag< segment_handle<SegmentationT> >::type CellTag; 00815 return closest_points_on_boundary_point_to_any<typename CellTag::facet_tag>( point_accessor, p, segment ); 00816 } 00817 00818 00819 00820 template <typename PointAccessorT, typename WrappedConfigT1, typename EA2, typename WrappedConfigT2> 00821 std::pair<typename PointAccessorT::value_type, typename PointAccessorT::value_type> 00822 closest_points_on_boundary_impl(PointAccessorT const accessor, 00823 viennagrid::element<vertex_tag,WrappedConfigT1> const & v, 00824 viennagrid::element<EA2,WrappedConfigT2> const & el) 00825 { 00826 return closest_points_on_boundary_impl( accessor(v), el); 00827 } 00828 00829 template <typename WrappedMeshConfigType, typename WrappedConfigT> 00830 std::pair< typename viennagrid::result_of::point< mesh<WrappedMeshConfigType> >::type, typename viennagrid::result_of::point< mesh<WrappedMeshConfigType> >::type > 00831 closest_points_on_boundary_impl(viennagrid::element<vertex_tag,WrappedConfigT> const & v, 00832 mesh<WrappedMeshConfigType> const & mesh_obj) 00833 { 00834 return closest_points_on_boundary_impl( default_point_accessor(mesh_obj)(v), mesh_obj); 00835 } 00836 00837 00839 00840 template <typename PointAccessorT, typename SomethingT1, typename SomethingT2> 00841 std::pair< typename PointAccessorT::value_type, typename PointAccessorT::value_type > 00842 closest_points_on_boundary_generic(PointAccessorT const accessor, 00843 SomethingT1 const & el1, 00844 SomethingT2 const & el2) 00845 { 00846 00847 00848 typedef typename viennagrid::result_of::facet_tag<SomethingT1>::type FacetTag1; 00849 typedef typename viennagrid::result_of::facet_tag<SomethingT2>::type FacetTag2; 00850 00851 00852 typedef std::pair< typename PointAccessorT::value_type, typename PointAccessorT::value_type > PairType; 00853 00854 typedef typename viennagrid::result_of::const_element_range<SomethingT1, FacetTag1>::type FacetRange1; 00855 typedef typename viennagrid::result_of::iterator<FacetRange1>::type FacetIterator1; 00856 00857 typedef typename viennagrid::result_of::const_element_range<SomethingT2, FacetTag2>::type FacetRange2; 00858 typedef typename viennagrid::result_of::iterator<FacetRange2>::type FacetIterator2; 00859 00860 PairType closest_pair; 00861 double shortest_distance = std::numeric_limits<double>::max(); 00862 00863 00864 FacetRange1 facets1(el1); 00865 FacetRange2 facets2(el2); 00866 00867 for (FacetIterator1 fit1 = facets1.begin(); 00868 fit1 != facets1.end(); 00869 ++fit1) 00870 { 00871 if (!is_boundary(el1, *fit1)) 00872 continue; 00873 00874 for (FacetIterator2 fit2 = facets2.begin(); 00875 fit2 != facets2.end(); 00876 ++fit2) 00877 { 00878 if (!is_boundary(el2, *fit2)) 00879 continue; 00880 00881 PairType p = closest_points_impl(accessor, *fit1, *fit2); 00882 double cur_norm = norm_2(p.first - p.second); 00883 if (cur_norm < shortest_distance) 00884 { 00885 closest_pair = p; 00886 shortest_distance = cur_norm; 00887 } 00888 } 00889 } 00890 00891 return closest_pair; 00892 } 00893 00894 00895 template <typename PointAccessorT, 00896 typename ElementTag1, typename WrappedConfigT1, 00897 typename ElementTag2, typename WrappedConfigT2> 00898 std::pair< typename PointAccessorT::value_type, typename PointAccessorT::value_type > 00899 closest_points_on_boundary_impl(PointAccessorT const accessor, 00900 viennagrid::element<ElementTag1,WrappedConfigT1> const & el1, 00901 viennagrid::element<ElementTag2,WrappedConfigT2> const & el2) 00902 { 00903 return closest_points_on_boundary_generic(accessor, el1, el2); 00904 } 00905 00906 template <typename PointAccessorT, 00907 typename WrappedMeshConfigType, 00908 typename ElementTag, typename WrappedConfigT> 00909 std::pair< typename PointAccessorT::value_type, typename PointAccessorT::value_type > 00910 closest_points_on_boundary_impl(PointAccessorT const accessor, 00911 viennagrid::element<ElementTag,WrappedConfigT> const & el1, 00912 mesh<WrappedMeshConfigType> const & mesh_obj) 00913 { 00914 return closest_points_on_boundary_generic(accessor, mesh_obj, el1); 00915 } 00916 00917 template <typename PointAccessorT, 00918 typename SegmentationT, 00919 typename ElementTag, typename WrappedConfigT> 00920 std::pair< typename PointAccessorT::value_type, typename PointAccessorT::value_type > 00921 closest_points_on_boundary_impl(PointAccessorT const accessor, 00922 viennagrid::element<ElementTag,WrappedConfigT> const & el1, 00923 segment_handle<SegmentationT> const & segment) 00924 { 00925 return closest_points_on_boundary_generic(accessor, segment, el1); 00926 } 00927 00928 00929 template <typename PointAccessorT, 00930 typename Segmentation1T, 00931 typename Segmentation2T> 00932 std::pair< typename PointAccessorT::value_type, typename PointAccessorT::value_type > 00933 closest_points_on_boundary_impl(PointAccessorT const accessor, 00934 segment_handle<Segmentation1T> const & segment1, 00935 segment_handle<Segmentation2T> const & segment2) 00936 { 00937 return closest_points_on_boundary_generic(accessor, segment1, segment2); 00938 } 00939 00940 00941 00943 00944 template <typename T> 00945 struct topological_id 00946 { 00947 //by default, typename is unknown, thus force error by not defining 'value' 00948 }; 00949 00950 template <typename CoordT, typename CoordinateSystemT> 00951 struct topological_id< spatial_point<CoordT, CoordinateSystemT> > 00952 { 00953 static const int value = 1; 00954 }; 00955 00956 00957 template <int dim, typename WrappedConfigT> 00958 struct topological_id< viennagrid::element<simplex_tag<dim>, WrappedConfigT> > 00959 { 00960 static const int value = 10000 + dim; //10.000 dimensions are certainly far from being ever instantiated 00961 }; 00962 00963 template <int dim, typename WrappedConfigT> 00964 struct topological_id< viennagrid::element<hypercube_tag<dim>, WrappedConfigT> > 00965 { 00966 static const int value = 20000 + dim; 00967 }; 00968 00969 template <typename WrappedConfigT> 00970 struct topological_id< viennagrid::mesh<WrappedConfigT> > 00971 { 00972 static const int value = 100000; 00973 }; 00974 00975 template <typename SegmentationT> 00976 struct topological_id< segment_handle<SegmentationT> > 00977 { 00978 static const int value = 100000 + 1; 00979 }; 00980 00981 00982 00983 template <typename T, typename U> 00984 struct topologically_sorted 00985 { 00986 static const bool value = (topological_id<T>::value <= topological_id<U>::value) ? true : false; 00987 }; 00988 00989 00990 template <typename T, 00991 typename U, 00992 bool correct_order = topologically_sorted<T, U>::value > 00993 struct ascending_topological_order 00994 { 00995 static T const & first(T const & t, U const &) { return t; } 00996 static U const & second(T const &, U const & u) { return u; } 00997 }; 00998 00999 template <typename T, 01000 typename U> 01001 struct ascending_topological_order<T, U, false> 01002 { 01003 static U const & first(T const &, U const & u) { return u; } 01004 static T const & second(T const & t, U const &) { return t; } 01005 }; 01006 01007 } //namespace detail 01008 01009 // 01010 // The public interface functions 01011 // 01013 template <typename PointAccessorT, typename SomethingT1, typename SomethingT2> 01014 std::pair< typename PointAccessorT::value_type, typename PointAccessorT::value_type > 01015 closest_points(PointAccessorT const accessor, 01016 SomethingT1 const & el1, 01017 SomethingT2 const & el2) 01018 { 01019 return detail::closest_points_impl(accessor, 01020 detail::ascending_topological_order<SomethingT1, SomethingT2>::first(el1, el2), 01021 detail::ascending_topological_order<SomethingT1, SomethingT2>::second(el1, el2)); 01022 } 01023 01025 template <typename SomethingT1, typename SomethingT2> 01026 std::pair< typename viennagrid::result_of::point<SomethingT1>::type, typename viennagrid::result_of::point<SomethingT1>::type > 01027 closest_points(SomethingT1 const & el1, 01028 SomethingT2 const & el2) 01029 { 01030 return closest_points_impl(default_point_accessor(el1), 01031 detail::ascending_topological_order<SomethingT1, SomethingT2>::first(el1, el2), 01032 detail::ascending_topological_order<SomethingT1, SomethingT2>::second(el1, el2)); 01033 } 01034 01035 01037 template <typename PointAccessorT, typename SomethingT1, typename SomethingT2> 01038 std::pair< typename PointAccessorT::value_type, typename PointAccessorT::value_type > 01039 closest_points_on_boundary(PointAccessorT const accessor, 01040 SomethingT1 const & el1, 01041 SomethingT2 const & el2) 01042 { 01043 return detail::closest_points_on_boundary_impl(accessor, 01044 detail::ascending_topological_order<SomethingT1, SomethingT2>::first(el1, el2), 01045 detail::ascending_topological_order<SomethingT1, SomethingT2>::second(el1, el2)); 01046 } 01047 01049 template <typename SomethingT1, typename SomethingT2> 01050 std::pair< typename viennagrid::result_of::point<SomethingT1>::type, typename viennagrid::result_of::point<SomethingT1>::type > 01051 closest_points_on_boundary(SomethingT1 const & el1, 01052 SomethingT2 const & el2) 01053 { 01054 return detail::closest_points_on_boundary_impl(default_point_accessor(el1), 01055 detail::ascending_topological_order<SomethingT1, SomethingT2>::first(el1, el2), 01056 detail::ascending_topological_order<SomethingT1, SomethingT2>::second(el1, el2)); 01057 } 01058 01059 01060 } //namespace viennagrid 01061 #endif
1.7.6.1