Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
function_signed_distance.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2019 - 2023 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
9// Public License as published by the Free Software Foundation; either
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
19
20#include <algorithm>
21
23
24namespace Functions
25{
26 namespace SignedDistance
27 {
28 template <int dim>
29 Sphere<dim>::Sphere(const Point<dim> &center, const double radius)
30 : center(center)
31 , radius(radius)
32 {
33 Assert(radius > 0, ExcMessage("Radius must be positive."))
34 }
35
36
37
38 template <int dim>
39 double
41 const unsigned int component) const
42 {
43 AssertIndexRange(component, this->n_components);
44 (void)component;
45
46 return point.distance(center) - radius;
47 }
48
49
50
51 template <int dim>
54 const unsigned int component) const
55 {
56 AssertIndexRange(component, this->n_components);
57 (void)component;
58
59 const Tensor<1, dim> center_to_point = point - center;
60 const Tensor<1, dim> grad = center_to_point / center_to_point.norm();
61 return grad;
62 }
63
64
65
66 template <int dim>
69 const unsigned int component) const
70 {
71 AssertIndexRange(component, this->n_components);
72 (void)component;
73
74 const Tensor<1, dim> center_to_point = point - center;
75 const double distance = center_to_point.norm();
76
77 const SymmetricTensor<2, dim> hess =
78 unit_symmetric_tensor<dim>() / distance -
79 symmetrize(outer_product(center_to_point, center_to_point)) /
80 Utilities::fixed_power<3>(distance);
81
82 return hess;
83 }
84
85
86
87 template <int dim>
88 Plane<dim>::Plane(const Point<dim> &point, const Tensor<1, dim> &normal)
89 : point_in_plane(point)
90 , normal(normal)
91 {
92 Assert(normal.norm() > 0, ExcMessage("Plane normal must not be 0."))
93 }
94
95
96
97 template <int dim>
98 double
100 const unsigned int component) const
101 {
102 AssertIndexRange(component, this->n_components);
103 (void)component;
104
105 return normal * (point - point_in_plane);
106 }
107
108
109
110 template <int dim>
112 Plane<dim>::gradient(const Point<dim> &, const unsigned int component) const
113 {
114 AssertIndexRange(component, this->n_components);
115 (void)component;
116
117 return normal;
118 }
119
120
121
122 template <int dim>
124 Plane<dim>::hessian(const Point<dim> &, const unsigned int component) const
125 {
126 AssertIndexRange(component, this->n_components);
127 (void)component;
128
130 }
131
132
133
134 template <int dim>
136 const std::array<double, dim> &radii,
137 const double tolerance,
138 const unsigned int max_iter)
139 : center(center)
140 , radii(radii)
141 , tolerance(tolerance)
142 , max_iter(max_iter)
143 {
144 for (unsigned int d = 0; d < dim; ++d)
145 Assert(radii[d] > 0, ExcMessage("All radii must be positive."))
146 }
147
148
149
150 template <int dim>
151 double
153 const unsigned int component) const
154 {
155 AssertIndexRange(component, this->n_components);
156 (void)component;
157
158 if (dim == 1)
159 return point.distance(center) - radii[0];
160 else if (dim == 2)
161 return compute_signed_distance_ellipse(point);
162 else
163 Assert(false, ExcNotImplemented());
164
165 return 0.0;
166 }
167
168
169
170 template <int dim>
173 const unsigned int component) const
174 {
175 AssertIndexRange(component, this->n_components);
176 (void)component;
177
178 Tensor<1, dim> grad;
179 if (dim == 1)
180 grad = point - center;
181 else if (dim == 2)
182 {
183 const Point<dim> point_in_centered_coordinate_system =
184 Point<dim>(compute_closest_point_ellipse(point) - center);
185 grad = compute_analyical_normal_vector_on_ellipse(
186 point_in_centered_coordinate_system);
187 }
188 else
190
191 if (grad.norm() > 1e-12)
192 return grad / grad.norm();
193 else
194 return grad;
195 }
196
197
198
199 template <int dim>
200 double
202 {
203 double val = 0.0;
204 for (unsigned int d = 0; d < dim; ++d)
205 val += Utilities::fixed_power<2>((point[d] - center[d]) / radii[d]);
206 return val - 1.0;
207 }
208
209
210
211 template <int dim>
214 {
215 AssertDimension(dim, 2);
216
217 /*
218 * Function to compute the closest point on an ellipse (adopted from
219 * https://wet-robots.ghost.io/simple-method-for-distance-to-ellipse/ and
220 * https://github.com/0xfaded/ellipse_demo):
221 *
222 * Since the ellipse is symmetric to the two major axes through its
223 * center, the point is moved so the center coincides with the origin and
224 * into the first quadrant.
225 * 1. Choose a point on the ellipse (x), here x = a*cos(pi/4) and y =
226 * b*sin(pi/4).
227 * 2. Find second point on the ellipse, that has the same distance.
228 * 3. Find midpoint on the ellipse (must be closer).
229 * 4. Repeat 2.-4. until convergence.
230 */
231 // get equivalent point in first quadrant of centered ellipse
232 const double px = std::abs(point[0] - center[0]);
233 const double py = std::abs(point[1] - center[1]);
234 const double sign_px = std::copysign(1.0, point[0] - center[0]);
235 const double sign_py = std::copysign(1.0, point[1] - center[1]);
236 // get semi axes radii
237 const double &a = radii[0];
238 const double &b = radii[1];
239 // initial guess (t = angle from x-axis)
240 double t = numbers::PI_4;
241 double x = a * std::cos(t);
242 double y = b * std::sin(t);
243
244 unsigned int iter = 0;
245 double delta_t;
246 do
247 {
248 // compute the ellipse evolute (center of curvature) for the current t
249 const double ex =
250 (a * a - b * b) * Utilities::fixed_power<3>(std::cos(t)) / a;
251 const double ey =
252 (b * b - a * a) * Utilities::fixed_power<3>(std::sin(t)) / b;
253 // compute distances from current point on ellipse to its evolute
254 const double rx = x - ex;
255 const double ry = y - ey;
256 // compute distances from point to the current evolute
257 const double qx = px - ex;
258 const double qy = py - ey;
259 // compute the curvature radius at the current point on the ellipse
260 const double r = std::hypot(rx, ry);
261 // compute the distance from evolute to the point
262 const double q = std::hypot(qx, qy);
263 // compute step size on ellipse
264 const double delta_c = r * std::asin((rx * qy - ry * qx) / (r * q));
265 // compute approximate angle step
266 delta_t = delta_c / std::sqrt(a * a + b * b - x * x - y * y);
267 t += delta_t;
268 // make sure the angle stays in first quadrant
269 t = std::min(numbers::PI_2, std::max(0.0, t));
270 x = a * std::cos(t);
271 y = b * std::sin(t);
272 iter++;
273 }
274 while (std::abs(delta_t) > tolerance && iter < max_iter);
275 AssertIndexRange(iter, max_iter);
276
279
280 return center + Point<dim>(sign_px * x, sign_py * y);
281 }
282
283
284
285 template <int dim>
288 const Point<dim> &) const
289 {
291 return Tensor<1, dim, double>();
292 }
293
294
295
296 template <>
299 const Point<2> &point) const
300 {
301 const auto &a = radii[0];
302 const auto &b = radii[1];
303 const auto &x = point[0];
304 const auto &y = point[1];
305 return Tensor<1, 2, double>({b * x / a, a * y / b});
306 }
307
308
309
310 template <int dim>
311 double
313 {
315 return 0;
316 }
317
318
319
320 template <>
321 double
323 {
324 // point corresponds to center
325 if (point.distance(center) < tolerance)
326 return *std::min_element(radii.begin(), radii.end()) * -1.;
327
328 const Point<2> &closest_point = compute_closest_point_ellipse(point);
329
330 const double distance =
331 std::hypot(closest_point[0] - point[0], closest_point[1] - point[1]);
332
333 return evaluate_ellipsoid(point) < 0.0 ? -distance : distance;
334 }
335
336
337
338 template <int dim>
340 const Point<dim> &top_right)
341 : bounding_box({bottom_left, top_right})
342 {}
343
344
345
346 template <int dim>
348 : bounding_box(bounding_box)
349 {}
350
351
352
353 template <int dim>
354 double
356 const unsigned int component) const
357 {
358 AssertDimension(component, 0);
359 (void)component;
360
361 return bounding_box.signed_distance(p);
362 }
363
364
365
366 template <int dim>
368 const double radius,
369 const double notch_width,
370 const double notch_height)
371 : sphere(center, radius)
372 , notch(// bottom left
373 (dim == 1) ?
374 Point<dim>(center[0] - 0.5 * notch_width) :
375 (dim == 2) ?
376 Point<dim>(center[0] - 0.5 * notch_width,
377 std::numeric_limits<double>::lowest()) : /* notch is open in negative y-direction*/
378 Point<dim>(center[0] - 0.5 * notch_width,
379 std::numeric_limits<
380 double>::lowest() /* notch is open in negative y-direction*/,
381 std::numeric_limits<double>::lowest()), /* notch is open in negative z-direction*/
382 // top right
383 (dim == 1) ?
384 Point<dim>(center[0] + 0.5 * notch_width) :
385 (dim == 2) ?
386 Point<dim>(center[0] + 0.5 * notch_width,
387 center[1] + notch_height - radius) :
388 Point<dim>(center[0] + 0.5 * notch_width,
389 std::numeric_limits<
390 double>::max() /* notch is open in y-direction*/,
391 center[2] + notch_height - radius))
392 {
393 Assert(
394 notch_width <= 2 * radius,
396 "The width of the notch must be less than the circle diameter."));
397 Assert(
398 notch_height <= 2 * radius,
400 "The height of the notch must be less than the circle diameter."));
401 }
402
403
404
405 template <int dim>
406 double
408 const unsigned int component) const
409 {
410 (void)component;
411 AssertDimension(component, 0);
412
413 // calculate the set difference between the level set functions of the
414 // sphere and the notch
415 return std::max(sphere.value(p), -notch.value(p));
416 }
417 } // namespace SignedDistance
418} // namespace Functions
419
420#include "function_signed_distance.inst"
421
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
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
double value(const Point< dim > &p, const unsigned int component=0) const override
Rectangle(const Point< dim > &bottom_left, const Point< dim > &top_right)
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
ZalesakDisk(const Point< dim > &center, const double radius, const double notch_width, const double notch_height)
double value(const Point< dim > &p, const unsigned int component=0) const override
Definition point.h:112
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
numbers::NumberTraits< Number >::real_type norm() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
Point< 3 > center
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIsFinite(number)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
static constexpr double PI_2
Definition numbers.h:264
static constexpr double PI_4
Definition numbers.h:269
STL namespace.
::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 > abs(const ::VectorizedArray< Number, width > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)