Line |
Branch |
Exec |
Source |
1 |
|
|
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- |
2 |
|
|
// vi: set et ts=4 sw=4 sts=4: |
3 |
|
|
// |
4 |
|
|
// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder |
5 |
|
|
// SPDX-License-Identifier: GPL-3.0-or-later |
6 |
|
|
// |
7 |
|
|
/*! |
8 |
|
|
* \file |
9 |
|
|
* \ingroup Geometry |
10 |
|
|
* \brief A function to compute the circumsphere of a given geometry |
11 |
|
|
* \note The circumscribed sphere (or circumsphere) is a sphere touching all vertices of the geometry |
12 |
|
|
* Therefore it does not necessarily exist for all geometries. |
13 |
|
|
*/ |
14 |
|
|
#ifndef DUMUX_GEOMETRY_CIRCUMSPHERE_HH |
15 |
|
|
#define DUMUX_GEOMETRY_CIRCUMSPHERE_HH |
16 |
|
|
|
17 |
|
|
#include <type_traits> |
18 |
|
|
#include <dumux/common/math.hh> |
19 |
|
|
#include <dumux/geometry/sphere.hh> |
20 |
|
|
|
21 |
|
|
namespace Dumux { |
22 |
|
|
|
23 |
|
|
/*! |
24 |
|
|
* \ingroup Geometry |
25 |
|
|
* \brief Computes the circumsphere of a triangle |
26 |
|
|
* \note See https://www.ics.uci.edu/~eppstein/junkyard/circumcenter.html |
27 |
|
|
*/ |
28 |
|
|
template<class Point> |
29 |
|
|
static inline Sphere<typename Point::value_type, Point::dimension> |
30 |
|
240 |
circumSphereOfTriangle(const Point& a, const Point& b, const Point& c) |
31 |
|
|
{ |
32 |
|
240 |
const auto ac = c - a; |
33 |
|
240 |
const auto ab = b - a; |
34 |
|
240 |
const auto n = crossProduct(ab, ac); |
35 |
|
1440 |
const auto distCenterToA = (crossProduct(n, ab)*ac.two_norm2() + crossProduct(ac, n)*ab.two_norm2()) / (2.0*n.two_norm2()); |
36 |
|
|
|
37 |
|
240 |
return { a + distCenterToA, distCenterToA.two_norm() }; |
38 |
|
|
} |
39 |
|
|
|
40 |
|
|
/*! |
41 |
|
|
* \ingroup Geometry |
42 |
|
|
* \brief Computes the circumsphere of a triangle |
43 |
|
|
*/ |
44 |
|
|
template<class Geometry, typename std::enable_if_t<Geometry::mydimension == 2, int> = 0> |
45 |
|
|
static inline Sphere<typename Geometry::ctype, Geometry::coorddimension> |
46 |
|
240 |
circumSphereOfTriangle(const Geometry& triangle) |
47 |
|
|
{ |
48 |
|
480 |
assert(triangle.corners() == 3 && "Geometry is not a triangle."); |
49 |
|
600 |
return circumSphereOfTriangle(triangle.corner(0), triangle.corner(1), triangle.corner(2)); |
50 |
|
|
} |
51 |
|
|
|
52 |
|
|
} // end namespace Dumux |
53 |
|
|
|
54 |
|
|
#endif |
55 |
|
|
|