3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
geometry/intersectspointgeometry.hh
Go to the documentation of this file.
1/*****************************************************************************
2 * See the file COPYING for full copying permissions. *
3 * *
4 * This program is free software: you can redistribute it and/or modify *
5 * it under the terms of the GNU General Public License as published by *
6 * the Free Software Foundation, either version 3 of the License, or *
7 * (at your option) any later version. *
8 * *
9 * This program is distributed in the hope that it will be useful, *
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
12 * GNU General Public License for more details. *
13 * *
14 * You should have received a copy of the GNU General Public License *
15 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
16 *****************************************************************************/
22#ifndef DUMUX_GEOMETRY_INTERSECTS_POINT_GEOMETRY_HH
23#define DUMUX_GEOMETRY_INTERSECTS_POINT_GEOMETRY_HH
24
25#include <cmath>
26#include <dune/common/exceptions.hh>
27#include <dune/common/fvector.hh>
28
30
31namespace Dumux {
32
37template <class ctype, int dimworld, class Geometry, typename std::enable_if_t<(Geometry::mydimension == 3), int> = 0>
38bool intersectsPointGeometry(const Dune::FieldVector<ctype, dimworld>& point, const Geometry& g)
39{
40 // get the g type
41 const auto type = g.type();
42
43 // if it's a tetrahedron we can check directly
44 if (type.isTetrahedron())
45 return intersectsPointSimplex(point, g.corner(0), g.corner(1), g.corner(2), g.corner(3));
46
47 // split hexahedrons into five tetrahedrons
48 else if (type.isHexahedron())
49 {
50 if (intersectsPointSimplex(point, g.corner(0), g.corner(1), g.corner(3), g.corner(5))) return true;
51 if (intersectsPointSimplex(point, g.corner(0), g.corner(5), g.corner(6), g.corner(4))) return true;
52 if (intersectsPointSimplex(point, g.corner(5), g.corner(3), g.corner(6), g.corner(7))) return true;
53 if (intersectsPointSimplex(point, g.corner(0), g.corner(3), g.corner(2), g.corner(6))) return true;
54 if (intersectsPointSimplex(point, g.corner(5), g.corner(3), g.corner(0), g.corner(6))) return true;
55 return false;
56 }
57
58 // split pyramids into two tetrahedrons
59 else if (type.isPyramid())
60 {
61 if (intersectsPointSimplex(point, g.corner(0), g.corner(1), g.corner(2), g.corner(4))) return true;
62 if (intersectsPointSimplex(point, g.corner(1), g.corner(3), g.corner(2), g.corner(4))) return true;
63 return false;
64 }
65
66 // split prisms into three tetrahedrons
67 else if (type.isPrism())
68 {
69 if (intersectsPointSimplex(point, g.corner(0), g.corner(1), g.corner(2), g.corner(4))) return true;
70 if (intersectsPointSimplex(point, g.corner(3), g.corner(0), g.corner(2), g.corner(4))) return true;
71 if (intersectsPointSimplex(point, g.corner(2), g.corner(5), g.corner(3), g.corner(4))) return true;
72 return false;
73 }
74
75 else
76 DUNE_THROW(Dune::NotImplemented,
77 "Intersection for point and geometry type "
78 << type << " in " << dimworld << "-dimensional world.");
79}
80
85template <class ctype, int dimworld, class Geometry, typename std::enable_if_t<(Geometry::mydimension == 2), int> = 0>
86bool intersectsPointGeometry(const Dune::FieldVector<ctype, dimworld>& point, const Geometry& g)
87{
88 // get the g type
89 const auto type = g.type();
90
91 // if it's a triangle we can check directly
92 if (type.isTriangle())
93 return intersectsPointSimplex(point, g.corner(0), g.corner(1), g.corner(2));
94
95 // split quadrilaterals into two triangles
96 else if (type.isQuadrilateral())
97 {
98 if (intersectsPointSimplex(point, g.corner(0), g.corner(1), g.corner(3))) return true;
99 if (intersectsPointSimplex(point, g.corner(0), g.corner(3), g.corner(2))) return true;
100 return false;
101 }
102
103 else
104 DUNE_THROW(Dune::NotImplemented,
105 "Intersection for point and geometry type "
106 << type << " in " << dimworld << "-dimensional world.");
107}
108
113template <class ctype, int dimworld, class Geometry, typename std::enable_if_t<(Geometry::mydimension == 1), int> = 0>
114bool intersectsPointGeometry(const Dune::FieldVector<ctype, dimworld>& point, const Geometry& g)
115{
116 return intersectsPointSimplex(point, g.corner(0), g.corner(1));
117}
118
119} // end namespace Dumux
120
121#endif
bool intersectsPointGeometry(const Dune::FieldVector< ctype, dimworld > &point, const Geometry &g)
Find out whether a point is inside a three-dimensional geometry.
Definition: geometry/intersectspointgeometry.hh:38
bool intersectsPointSimplex(const Dune::FieldVector< ctype, dimworld > &point, const Dune::FieldVector< ctype, dimworld > &p0, const Dune::FieldVector< ctype, dimworld > &p1, const Dune::FieldVector< ctype, dimworld > &p2, const Dune::FieldVector< ctype, dimworld > &p3)
Find out whether a point is inside a tetrahedron (p0, p1, p2, p3) (dimworld is 3)
Definition: geometry/intersectspointsimplex.hh:36
Definition: adapt.hh:29
Detect if a point intersects a simplex (including boundary)