Reference documentation for deal.II version GIT db2fd67796 2022-12-07 23:35: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 - 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
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 
19 #include <algorithm>
20 
22 
23 namespace 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  std::pow(distance, 3);
80 
81  return hess;
82  }
83 
84 
85 
86  template <int dim>
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 
128  return SymmetricTensor<2, dim>();
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
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 
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
188  AssertThrow(false, ExcNotImplemented());
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 += std::pow((point[d] - center[d]) / radii[d], 2);
205  return val - 1.0;
206  }
207 
208 
209 
210  template <int dim>
211  Point<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 = (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 
274  AssertIsFinite(x);
275  AssertIsFinite(y);
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  {
287  AssertThrow(false, ExcNotImplemented());
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  {
311  AssertThrow(false, ExcNotImplemented());
312  return 0;
313  }
314 
315 
316 
317  template <>
318  double
320  {
321  // point corresponds to center
322  if (point.distance(center) < tolerance)
323  return *std::min_element(radii.begin(), radii.end()) * -1.;
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
numbers::NumberTraits< Number >::real_type norm() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
Point< 3 > center
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcNotImplemented()
#define AssertIsFinite(number)
Definition: exceptions.h:1786
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1695
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1611
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:189
Point< dim > closest_point(const TopoDS_Shape &in_shape, const Point< dim > &origin, const double tolerance=1e-7)
Definition: utilities.cc:786
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:253
static constexpr double PI_4
Definition: numbers.h:258
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)