Reference documentation for deal.II version 9.4.1
Searching...
No Matches
function_signed_distance.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2019 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
18
19#include <algorithm>
20
22
23namespace Functions
24{
25 namespace SignedDistance
26 {
27 template <int dim>
28 Sphere<dim>::Sphere(const Point<dim> &center, const double radius)
29 : center(center)
31 {
33 }
34
35
36
37 template <int dim>
38 double
40 const unsigned int component) const
41 {
42 AssertIndexRange(component, this->n_components);
43 (void)component;
44
46 }
47
48
49
50 template <int dim>
53 const unsigned int component) const
54 {
55 AssertIndexRange(component, this->n_components);
56 (void)component;
57
58 const Tensor<1, dim> center_to_point = point - center;
59 const Tensor<1, dim> grad = center_to_point / center_to_point.norm();
61 }
62
63
64
65 template <int dim>
68 const unsigned int component) const
69 {
70 AssertIndexRange(component, this->n_components);
71 (void)component;
72
73 const Tensor<1, dim> center_to_point = point - center;
74 const double distance = center_to_point.norm();
75
76 const SymmetricTensor<2, dim> hess =
77 unit_symmetric_tensor<dim>() / distance -
78 symmetrize(outer_product(center_to_point, center_to_point)) /
79 std::pow(distance, 3);
80
81 return hess;
82 }
83
84
85
86 template <int dim>
87 Plane<dim>::Plane(const Point<dim> &point, const Tensor<1, dim> &normal)
88 : point_in_plane(point)
89 , normal(normal)
90 {
91 Assert(normal.norm() > 0, ExcMessage("Plane normal must not be 0."))
92 }
93
94
95
96 template <int dim>
97 double
99 const unsigned int component) const
100 {
101 AssertIndexRange(component, this->n_components);
102 (void)component;
103
104 return normal * (point - point_in_plane);
105 }
106
107
108
109 template <int dim>
111 Plane<dim>::gradient(const Point<dim> &, const unsigned int component) const
112 {
113 AssertIndexRange(component, this->n_components);
114 (void)component;
115
116 return normal;
117 }
118
119
120
121 template <int dim>
123 Plane<dim>::hessian(const Point<dim> &, const unsigned int component) const
124 {
125 AssertIndexRange(component, this->n_components);
126 (void)component;
127
129 }
130
131
132
133 template <int dim>
136 const double tolerance,
137 const unsigned int max_iter)
138 : center(center)
140 , tolerance(tolerance)
141 , max_iter(max_iter)
142 {
143 for (unsigned int d = 0; d < dim; ++d)
145 }
146
147
148
149 template <int dim>
150 double
152 const unsigned int component) const
153 {
154 AssertIndexRange(component, this->n_components);
155 (void)component;
156
157 if (dim == 1)
159 else if (dim == 2)
160 return compute_signed_distance_ellipse(point);
161 else
162 Assert(false, ExcNotImplemented());
163
164 return 0.0;
165 }
166
167
168
169 template <int dim>
172 const unsigned int component) const
173 {
174 AssertIndexRange(component, this->n_components);
175 (void)component;
176
178 if (dim == 1)
179 grad = point - center;
180 else if (dim == 2)
181 {
182 const Point<dim> point_in_centered_coordinate_system =
183 Point<dim>(compute_closest_point_ellipse(point) - center);
185 point_in_centered_coordinate_system);
186 }
187 else
189
192 else
194 }
195
196
197
198 template <int dim>
199 double
201 {
202 double val = 0.0;
203 for (unsigned int d = 0; d < dim; ++d)
204 val += std::pow((point[d] - center[d]) / radii[d], 2);
205 return val - 1.0;
206 }
207
208
209
210 template <int dim>
213 {
214 AssertDimension(dim, 2);
215
216 /*
217 * Function to compute the closest point on an ellipse (adopted from
218 * https://wet-robots.ghost.io/simple-method-for-distance-to-ellipse/ and
220 *
221 * Since the ellipse is symmetric to the two major axes through its
222 * center, the point is moved so the center coincides with the origin and
223 * into the first quadrant.
224 * 1. Choose a point on the ellipse (x), here x = a*cos(pi/4) and y =
225 * b*sin(pi/4).
226 * 2. Find second point on the ellipse, that has the same distance.
227 * 3. Find midpoint on the ellipse (must be closer).
228 * 4. Repeat 2.-4. until convergence.
229 */
230 // get equivalent point in first quadrant of centered ellipse
231 const double px = std::abs(point[0] - center[0]);
232 const double py = std::abs(point[1] - center[1]);
233 const double sign_px = std::copysign(1.0, point[0] - center[0]);
234 const double sign_py = std::copysign(1.0, point[1] - center[1]);
235 // get semi axes radii
236 const double &a = radii[0];
237 const double &b = radii[1];
238 // initial guess (t = angle from x-axis)
239 double t = numbers::PI_4;
240 double x = a * std::cos(t);
241 double y = b * std::sin(t);
242
243 unsigned int iter = 0;
244 double delta_t;
245 do
246 {
247 // compute the ellipse evolute (center of curvature) for the current t
248 const double ex = (a * a - b * b) * std::pow(std::cos(t), 3) / a;
249 const double ey = (b * b - a * a) * std::pow(std::sin(t), 3) / b;
250 // compute distances from current point on ellipse to its evolute
251 const double rx = x - ex;
252 const double ry = y - ey;
253 // compute distances from point to the current evolute
254 const double qx = px - ex;
255 const double qy = py - ey;
256 // compute the curvature radius at the current point on the ellipse
257 const double r = std::hypot(rx, ry);
258 // compute the distance from evolute to the point
259 const double q = std::hypot(qx, qy);
260 // compute step size on ellipse
261 const double delta_c = r * std::asin((rx * qy - ry * qx) / (r * q));
262 // compute approximate angle step
263 delta_t = delta_c / std::sqrt(a * a + b * b - x * x - y * y);
264 t += delta_t;
265 // make sure the angle stays in first quadrant
266 t = std::min(numbers::PI_2, std::max(0.0, t));
267 x = a * std::cos(t);
268 y = b * std::sin(t);
269 iter++;
270 }
271 while (std::abs(delta_t) > tolerance && iter < max_iter);
272 AssertIndexRange(iter, max_iter);
273
276
277 return center + Point<dim>(sign_px * x, sign_py * y);
278 }
279
280
281
282 template <int dim>
285 const Point<dim> &) const
286 {
288 return Tensor<1, dim, double>();
289 }
290
291
292
293 template <>
296 const Point<2> &point) const
297 {
298 const auto &a = radii[0];
299 const auto &b = radii[1];
300 const auto &x = point[0];
301 const auto &y = point[1];
302 return Tensor<1, 2, double>({b * x / a, a * y / b});
303 }
304
305
306
307 template <int dim>
308 double
310 {
312 return 0;
313 }
314
315
316
317 template <>
318 double
320 {
321 // point corresponds to center
322 if (point.distance(center) < tolerance)
324
325 const Point<2> &closest_point = compute_closest_point_ellipse(point);
326
327 const double distance =
328 std::hypot(closest_point[0] - point[0], closest_point[1] - point[1]);
329
330 return evaluate_ellipsoid(point) < 0.0 ? -distance : distance;
331 }
332 } // namespace SignedDistance
333} // namespace Functions
334
335#include "function_signed_distance.inst"
336
double value(const Point< dim > &point, const unsigned int component=0) const override
Tensor< 1, dim > gradient(const Point< dim > &, const unsigned int component=0) const override
double compute_signed_distance_ellipse(const Point< dim > &point) const
Point< dim > compute_closest_point_ellipse(const Point< dim > &point) const
Ellipsoid(const Point< dim > &center, const std::array< double, dim > &radii, const double tolerance=1e-14, const unsigned int max_iter=10)
double evaluate_ellipsoid(const Point< dim > &point) const
const std::array< double, dim > radii
Tensor< 1, dim, double > compute_analyical_normal_vector_on_ellipse(const Point< dim > &point) const
Plane(const Point< dim > &point, const Tensor< 1, dim > &normal)
double value(const Point< dim > &point, const unsigned int component=0) const override
SymmetricTensor< 2, dim > hessian(const Point< dim > &, const unsigned int component=0) const override
Tensor< 1, dim > gradient(const Point< dim > &, const unsigned int component=0) const override
Sphere(const Point< dim > &center=Point< dim >(), const double radius=1)
SymmetricTensor< 2, dim > hessian(const Point< dim > &point, const unsigned int component=0) const override
double value(const Point< dim > &point, const unsigned int component=0) const override
Tensor< 1, dim > gradient(const Point< dim > &point, const unsigned int component=0) const override
Definition: point.h:111
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
Definition: tensor.h:503
numbers::NumberTraits< Number >::real_type norm() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
Point< 3 > center
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertIsFinite(number)
Definition: exceptions.h:1758
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
static constexpr double PI_2
Definition: numbers.h:238
static constexpr double PI_4
Definition: numbers.h:243
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)