Reference documentation for deal.II version GIT relicensing-437-g81ec864850 2024-04-19 07:30:02+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2021 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
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)
30 , radius(radius)
31 {
32 Assert(radius > 0, ExcMessage("Radius must be positive."));
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
45 return point.distance(center) - radius;
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();
60 return grad;
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 Utilities::fixed_power<3>(distance);
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>
135 const std::array<double, dim> &radii,
136 const double tolerance,
137 const unsigned int max_iter)
138 : center(center)
139 , radii(radii)
140 , tolerance(tolerance)
141 , max_iter(max_iter)
142 {
143 for (unsigned int d = 0; d < dim; ++d)
144 Assert(radii[d] > 0, ExcMessage("All radii must be positive."));
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)
158 return point.distance(center) - radii[0];
159 else if (dim == 2)
160 return compute_signed_distance_ellipse(point);
161 else
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
177 Tensor<1, dim> grad;
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);
184 grad = compute_analyical_normal_vector_on_ellipse(
185 point_in_centered_coordinate_system);
186 }
187 else
189
190 if (grad.norm() > 1e-12)
191 return grad / grad.norm();
192 else
193 return grad;
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 += Utilities::fixed_power<2>((point[d] - center[d]) / radii[d]);
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
219 * https://github.com/0xfaded/ellipse_demo):
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 =
249 (a * a - b * b) * Utilities::fixed_power<3>(std::cos(t)) / a;
250 const double ey =
251 (b * b - a * a) * Utilities::fixed_power<3>(std::sin(t)) / b;
252 // compute distances from current point on ellipse to its evolute
253 const double rx = x - ex;
254 const double ry = y - ey;
255 // compute distances from point to the current evolute
256 const double qx = px - ex;
257 const double qy = py - ey;
258 // compute the curvature radius at the current point on the ellipse
259 const double r = std::hypot(rx, ry);
260 // compute the distance from evolute to the point
261 const double q = std::hypot(qx, qy);
262 // compute step size on ellipse
263 const double delta_c = r * std::asin((rx * qy - ry * qx) / (r * q));
264 // compute approximate angle step
265 delta_t = delta_c / std::sqrt(a * a + b * b - x * x - y * y);
266 t += delta_t;
267 // make sure the angle stays in first quadrant
268 t = std::min(numbers::PI_2, std::max(0.0, t));
269 x = a * std::cos(t);
270 y = b * std::sin(t);
271 ++iter;
272 }
273 while (std::abs(delta_t) > tolerance && iter < max_iter);
274 AssertIndexRange(iter, max_iter);
275
278
279 return center + Point<dim>(sign_px * x, sign_py * y);
280 }
281
282
283
284 template <int dim>
292
293
294
295 template <>
298 const Point<2> &point) const
299 {
300 const auto &a = radii[0];
301 const auto &b = radii[1];
302 const auto &x = point[0];
303 const auto &y = point[1];
304 return Tensor<1, 2, double>({b * x / a, a * y / b});
305 }
306
307
308
309 template <int dim>
310 double
312 {
314 return 0;
315 }
316
317
318
319 template <>
320 double
322 {
323 // point corresponds to center
324 if (point.distance(center) < tolerance)
325 return *std::min_element(radii.begin(), radii.end()) * -1.;
326
327 const Point<2> &closest_point = compute_closest_point_ellipse(point);
328
329 const double distance =
330 std::hypot(closest_point[0] - point[0], closest_point[1] - point[1]);
331
332 return evaluate_ellipsoid(point) < 0.0 ? -distance : distance;
333 }
334
335
336
337 template <int dim>
339 const Point<dim> &top_right)
340 : bounding_box({bottom_left, top_right})
341 {}
342
343
344
345 template <int dim>
347 : bounding_box(bounding_box)
348 {}
349
350
351
352 template <int dim>
353 double
355 const unsigned int component) const
356 {
357 AssertDimension(component, 0);
358 (void)component;
359
360 return bounding_box.signed_distance(p);
361 }
362
363
364
365 template <int dim>
367 const double radius,
368 const double notch_width,
369 const double notch_height)
370 : sphere(center, radius)
371 , notch(// bottom left
372 (dim == 1) ?
373 Point<dim>(center[0] - 0.5 * notch_width) :
374 (dim == 2) ?
375 Point<dim>(center[0] - 0.5 * notch_width,
376 std::numeric_limits<double>::lowest()) : /* notch is open in negative y-direction*/
377 Point<dim>(center[0] - 0.5 * notch_width,
378 std::numeric_limits<
379 double>::lowest() /* notch is open in negative y-direction*/,
380 std::numeric_limits<double>::lowest()), /* notch is open in negative z-direction*/
381 // top right
382 (dim == 1) ?
383 Point<dim>(center[0] + 0.5 * notch_width) :
384 (dim == 2) ?
385 Point<dim>(center[0] + 0.5 * notch_width,
386 center[1] + notch_height - radius) :
387 Point<dim>(center[0] + 0.5 * notch_width,
388 std::numeric_limits<
389 double>::max() /* notch is open in y-direction*/,
390 center[2] + notch_height - radius))
391 {
392 Assert(
393 notch_width <= 2 * radius,
395 "The width of the notch must be less than the circle diameter."));
396 Assert(
397 notch_height <= 2 * radius,
399 "The height of the notch must be less than the circle diameter."));
400 }
401
402
403
404 template <int dim>
405 double
407 const unsigned int component) const
408 {
409 (void)component;
410 AssertDimension(component, 0);
411
412 // calculate the set difference between the level set functions of the
413 // sphere and the notch
414 return std::max(sphere.value(p), -notch.value(p));
415 }
416 } // namespace SignedDistance
417} // namespace Functions
418
419#include "function_signed_distance.inst"
420
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:111
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
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)
#define DEAL_II_NOT_IMPLEMENTED()
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)