distance.hh 7.37 KB
 Ned Coltman committed Jul 29, 2020 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 ``````/***************************************************************************** * See the file COPYING for full copying permissions. * * * * This program is free software: you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation, either version 3 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program. If not, see . * *****************************************************************************/ /*! * \file * \ingroup Geometry * \brief Helper functions for distance queries */ #ifndef DUMUX_GEOMETRY_DISTANCE_HH #define DUMUX_GEOMETRY_DISTANCE_HH #include #include namespace Dumux { /*! * \ingroup Geometry * \brief Compute the average distance from a point to a geometry by integration */ template inline typename Geometry::ctype averageDistancePointGeometry(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry, std::size_t integrationOrder = 2) { typename Geometry::ctype avgDist = 0.0; const auto& quad = Dune::QuadratureRules::rule(geometry.type(), integrationOrder); for (const auto& qp : quad) avgDist += (geometry.global(qp.position())-p).two_norm()*qp.weight()*geometry.integrationElement(qp.position()); return avgDist/geometry.volume(); } /*! * \ingroup Geometry * \brief Compute the distance from a point to a line through the points a and b */ template inline typename Point::value_type distancePointLine(const Point& p, const Point& a, const Point& b) { const auto ab = b - a; const auto t = (p - a)*ab/ab.two_norm2(); auto proj = a; proj.axpy(t, ab); return (proj - p).two_norm(); } /*! * \ingroup Geometry * \brief Compute the distance from a point to a line given by a geometry with two corners * \note We currently lack the a representation of a line geometry. This convenience function * assumes a segment geometry (with two corners) is passed which represents a line geometry. */ template inline typename Geometry::ctype distancePointLine(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry) { static_assert(Geometry::mydimension == 1, "Geometry has to be a line"); const auto& a = geometry.corner(0); const auto& b = geometry.corner(1); return distancePointLine(p, a, b); } /*! * \ingroup Geometry * \brief Compute the distance from a point to the segment connecting the points a and b */ template inline typename Point::value_type distancePointSegment(const Point& p, const Point& a, const Point& b) { const auto ab = b - a; `````` Timo Koch committed Jul 22, 2021 87 88 `````` const auto ap = p - a; const auto t = ap*ab; `````` Ned Coltman committed Jul 29, 2020 89 90 `````` if (t <= 0.0) `````` Timo Koch committed Jul 22, 2021 91 `````` return ap.two_norm(); `````` Ned Coltman committed Jul 29, 2020 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 `````` const auto lengthSq = ab.two_norm2(); if (t >= lengthSq) return (b - p).two_norm(); auto proj = a; proj.axpy(t/lengthSq, ab); return (proj - p).two_norm(); } /*! * \ingroup Geometry * \brief Compute the distance from a point to a given segment geometry */ template inline typename Geometry::ctype distancePointSegment(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry) { static_assert(Geometry::mydimension == 1, "Geometry has to be a segment"); const auto& a = geometry.corner(0); const auto& b = geometry.corner(1); return distancePointSegment(p, a, b); } `````` Timo Koch committed Oct 27, 2020 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 ``````/*! * \ingroup Geometry * \brief Compute the average distance from a segment to a geometry by integration */ template inline typename Geometry::ctype averageDistanceSegmentGeometry(const typename Geometry::GlobalCoordinate& a, const typename Geometry::GlobalCoordinate& b, const Geometry& geometry, std::size_t integrationOrder = 2) { typename Geometry::ctype avgDist = 0.0; const auto& quad = Dune::QuadratureRules::rule(geometry.type(), integrationOrder); for (const auto& qp : quad) avgDist += distancePointSegment(geometry.global(qp.position()), a, b)*qp.weight()*geometry.integrationElement(qp.position()); return avgDist/geometry.volume(); } `````` Ned Coltman committed Jul 29, 2020 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 ``````/*! * \ingroup Geometry * \brief Compute the shortest distance between two points */ template inline ctype distance(const Dune::FieldVector& a, const Dune::FieldVector& b) { return (a-b).two_norm(); } namespace Detail { // helper struct to compute distance between two geometries, specialized below template struct GeometryDistance { static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions"); static auto distance(const Geo1& geo1, const Geo2& geo2) { DUNE_THROW(Dune::NotImplemented, "Geometry distance computation not implemented for dimworld = " << dimWorld << ", dim1 = " << dim1 << ", dim2 = " << dim2); } }; // distance point-point template struct GeometryDistance { static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions"); static auto distance(const Geo1& geo1, const Geo2& geo2) { return Dumux::distance(geo1.corner(0), geo2.corner(0)); } }; // distance segment-point template struct GeometryDistance { static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions"); static auto distance(const Geo1& geo1, const Geo2& geo2) { return distancePointSegment(geo2.corner(0), geo1); } }; // distance point-segment template struct GeometryDistance { static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions"); static auto distance(const Geo1& geo1, const Geo2& geo2) { return distancePointSegment(geo1.corner(0), geo2); } }; } // end namespace Detail /*! * \ingroup Geometry * \brief Compute the shortest distance between two geometries */ template inline auto distance(const Geo1& geo1, const Geo2& geo2) { return Detail::GeometryDistance::distance(geo1, geo2); } } // end namespace Dumux #endif``````