Reference documentation for deal.II version GIT 9042b9283b 2023-12-02 14:50: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\}}\)
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 
18 #include <deal.II/base/utilities.h>
19 
20 #include <algorithm>
21 
23 
24 namespace 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>
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 
129  return SymmetricTensor<2, dim>();
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
189  AssertThrow(false, ExcNotImplemented());
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>
212  Point<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 
277  AssertIsFinite(x);
278  AssertIsFinite(y);
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  {
290  AssertThrow(false, ExcNotImplemented());
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  {
314  AssertThrow(false, ExcNotImplemented());
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,
395  ExcMessage(
396  "The width of the notch must be less than the circle diameter."));
397  Assert(
398  notch_height <= 2 * radius,
399  ExcMessage(
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
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
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:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
Point< 3 > center
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcNotImplemented()
#define AssertIsFinite(number)
Definition: exceptions.h:1915
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1732
Expression asin(const Expression &x)
Expression copysign(const Expression &value, const Expression &sign)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:192
Point< dim > closest_point(const TopoDS_Shape &in_shape, const Point< dim > &origin, const double tolerance=1e-7)
Definition: utilities.cc:789
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
static constexpr double PI_2
Definition: numbers.h:265
static constexpr double PI_4
Definition: numbers.h:270
constexpr DEAL_II_HOST SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
constexpr DEAL_II_HOST SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)