year={2016},
publisher={SIAM}
}
+
+@article{zalesak1979fully,
+ title={Fully multidimensional flux-corrected transport algorithms for fluids},
+ author={Zalesak, Steven T},
+ journal={Journal of Computational Physics},
+ volume={31},
+ number={3},
+ pages={335--362},
+ year={1979},
+ publisher={Elsevier}
+}
--- /dev/null
+New: Added Functions::SignedDistance::ZalesakDisk()
+to compute the signed distance function of Zalesak's disk.
+<br>
+(Magdalena Schreter, Simon Sticko, Peter Munch, 2023/04/06)
* @param top_right Top right point of the rectangle.
*/
Rectangle(const Point<dim> &bottom_left, const Point<dim> &top_right);
+ /**
+ * Constructor, takes a bounding box.
+ *
+ * @param bounding_box Bounding box.
+ */
+ Rectangle(const BoundingBox<dim> bounding_box);
/**
* Calculates the signed distance from a given point @p p to the rectangle.
private:
const BoundingBox<dim> bounding_box;
};
+
+
+ /**
+ * Signed-distance level set function of Zalesak's disk proposed in
+ * @cite zalesak1979fully.
+ *
+ * It is calculated by the set difference $\psi(x) = \max(\psi_{S}(x),
+ * -\psi_{N}(x))$ of the level set functions of a sphere $\psi_{S}$ and a
+ * rectangle $\psi_{N}$. This function is zero on the surface of the disk,
+ * negative "inside" and positive in the rest of $\mathbb{R}^{dim}$.
+ *
+ * Contour surfaces of the signed distance function of a 3D Zalesak's disk
+ * are illustrated below:
+ *
+ * @image html https://www.dealii.org/images/grids/signed_distance_zalesak_disk.png
+ *
+ * @ingroup functions
+ */
+ template <int dim>
+ class ZalesakDisk : public Function<dim>
+ {
+ public:
+ /**
+ * Constructor, takes the radius and center of the disk together
+ * with the width and the height of the notch.
+ *
+ * @param radius Radius of the disk.
+ * @param center Center of the disk.
+ * @param notch_width Width of the notch of the disk.
+ * @param notch_height Height of the notch of the disk.
+ */
+ ZalesakDisk(const Point<dim> ¢er,
+ const double radius,
+ const double notch_width,
+ const double notch_height);
+
+ /**
+ * Calculates the signed distance from a given point @p p to the disk.
+ * The signed distance is negative for points inside the disk, zero
+ * for points on the disk and positive for points outside the
+ * disk.
+ */
+ double
+ value(const Point<dim> & p,
+ const unsigned int component = 0) const override;
+
+ private:
+ /**
+ * Disk described by a sphere.
+ */
+ const Functions::SignedDistance::Sphere<dim> sphere;
+
+ /**
+ * Notch described by a rectangle.
+ */
+ const Functions::SignedDistance::Rectangle<dim> notch;
+ };
} // namespace SignedDistance
} // namespace Functions
+ template <int dim>
+ Rectangle<dim>::Rectangle(const BoundingBox<dim> bounding_box)
+ : bounding_box(bounding_box)
+ {}
+
+
+
template <int dim>
double
Rectangle<dim>::value(const Point<dim> & p,
return bounding_box.signed_distance(p);
}
+
+
+
+ template <int dim>
+ ZalesakDisk<dim>::ZalesakDisk(const Point<dim> ¢er,
+ const double radius,
+ const double notch_width,
+ const double notch_height)
+ : sphere(center, radius)
+ , notch(// bottom left
+ (dim == 1) ?
+ Point<dim>(center[0] - 0.5 * notch_width) :
+ (dim == 2) ?
+ Point<dim>(center[0] - 0.5 * notch_width,
+ std::numeric_limits<double>::lowest()) : /* notch is open in negative y-direction*/
+ Point<dim>(center[0] - 0.5 * notch_width,
+ std::numeric_limits<
+ double>::lowest() /* notch is open in negative y-direction*/,
+ std::numeric_limits<double>::lowest()), /* notch is open in negative z-direction*/
+ // top right
+ (dim == 1) ?
+ Point<dim>(center[0] + 0.5 * notch_width) :
+ (dim == 2) ?
+ Point<dim>(center[0] + 0.5 * notch_width,
+ center[1] + notch_height - radius) :
+ Point<dim>(center[0] + 0.5 * notch_width,
+ std::numeric_limits<
+ double>::max() /* notch is open in y-direction*/,
+ center[2] + notch_height - radius))
+ {
+ Assert(
+ notch_width <= 2 * radius,
+ ExcMessage(
+ "The width of the notch must be less than the circle diameter."));
+ Assert(
+ notch_height <= 2 * radius,
+ ExcMessage(
+ "The height of the notch must be less than the circle diameter."));
+ }
+
+
+
+ template <int dim>
+ double
+ ZalesakDisk<dim>::value(const Point<dim> & p,
+ const unsigned int component) const
+ {
+ (void)component;
+ AssertDimension(component, 0);
+
+ // calculate the set difference between the level set functions of the
+ // sphere and the notch
+ return std::max(sphere.value(p), -notch.value(p));
+ }
} // namespace SignedDistance
} // namespace Functions
template class Functions::SignedDistance::Plane<deal_II_dimension>;
template class Functions::SignedDistance::Ellipsoid<deal_II_dimension>;
template class Functions::SignedDistance::Rectangle<deal_II_dimension>;
+ template class Functions::SignedDistance::ZalesakDisk<deal_II_dimension>;
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2007 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+/**
+ * Test the classes Functions::SignedDistance::ZalesakDisk by creating an
+ * object of each class and calling the value functions at a point where we know
+ * what the function values should be.
+ */
+
+#include <deal.II/base/function_level_set.h>
+
+#include "../tests.h"
+
+
+namespace
+{
+ template <int dim>
+ void
+ print_values_at_point(const Function<dim> &function, const Point<dim> &point)
+ {
+ deallog << "point = " << point << std::endl;
+ deallog << "value = " << function.value(point) << std::endl;
+ }
+
+
+
+ template <int dim>
+ void
+ test_disk_level_set()
+ {
+ deallog << "test_zalesak_disk_level_set" << std::endl;
+
+ const Point<dim> center;
+ const double radius = 1;
+
+ const Functions::SignedDistance::ZalesakDisk<dim> level_set(center,
+ radius,
+ radius * 0.5,
+ 0.75 * radius);
+
+ {
+ Point<dim> p;
+
+ deallog << "center" << std::endl;
+ print_values_at_point(level_set, p);
+ deallog << "inside" << std::endl;
+ p[0] += 0.1;
+ print_values_at_point(level_set, p);
+ p[0] = 0.3;
+ print_values_at_point(level_set, p);
+ }
+ deallog << "notch" << std::endl;
+ Point<dim> p;
+ p[dim - 1] = -0.45;
+ p[std::max(0, dim - 2)] = -0.1;
+ print_values_at_point(level_set, p);
+ deallog << "outside" << std::endl;
+ p[dim - 1] = -2;
+ print_values_at_point(level_set, p);
+ p[dim - 1] = -1.2;
+ p[0] = -0.3;
+ p[std::max(0, dim - 2)] = -0.3;
+ print_values_at_point(level_set, p);
+ p[dim - 1] = -2;
+ p[0] = 0;
+ p[std::max(0, dim - 2)] = 0;
+ print_values_at_point(level_set, p);
+ p[dim - 1] = 2;
+ print_values_at_point(level_set, p);
+ p[0] = -2;
+ print_values_at_point(level_set, p);
+ p[0] = 2;
+ print_values_at_point(level_set, p);
+ }
+
+ template <int dim>
+ void
+ run_test()
+ {
+ deallog << "dim = " << dim << std::endl;
+ deallog << std::endl;
+ test_disk_level_set<dim>();
+ deallog << std::endl;
+ }
+} // namespace
+
+
+
+int
+main()
+{
+ initlog();
+ run_test<1>();
+ run_test<2>();
+ run_test<3>();
+}
--- /dev/null
+
+DEAL::dim = 1
+DEAL::
+DEAL::test_zalesak_disk_level_set
+DEAL::center
+DEAL::point = 0.00000
+DEAL::value = 0.250000
+DEAL::inside
+DEAL::point = 0.100000
+DEAL::value = 0.150000
+DEAL::point = 0.300000
+DEAL::value = -0.0500000
+DEAL::notch
+DEAL::point = -0.100000
+DEAL::value = 0.150000
+DEAL::outside
+DEAL::point = -2.00000
+DEAL::value = 1.00000
+DEAL::point = -0.300000
+DEAL::value = -0.0500000
+DEAL::point = 0.00000
+DEAL::value = 0.250000
+DEAL::point = 2.00000
+DEAL::value = 1.00000
+DEAL::point = -2.00000
+DEAL::value = 1.00000
+DEAL::point = 2.00000
+DEAL::value = 1.00000
+DEAL::
+DEAL::dim = 2
+DEAL::
+DEAL::test_zalesak_disk_level_set
+DEAL::center
+DEAL::point = 0.00000 0.00000
+DEAL::value = -0.250000
+DEAL::inside
+DEAL::point = 0.100000 0.00000
+DEAL::value = -0.250000
+DEAL::point = 0.300000 0.00000
+DEAL::value = -0.254951
+DEAL::notch
+DEAL::point = -0.100000 -0.450000
+DEAL::value = 0.150000
+DEAL::outside
+DEAL::point = -0.100000 -2.00000
+DEAL::value = 1.00250
+DEAL::point = -0.300000 -1.20000
+DEAL::value = 0.236932
+DEAL::point = 0.00000 -2.00000
+DEAL::value = 1.00000
+DEAL::point = 0.00000 2.00000
+DEAL::value = 1.00000
+DEAL::point = -2.00000 2.00000
+DEAL::value = 1.82843
+DEAL::point = 2.00000 2.00000
+DEAL::value = 1.82843
+DEAL::
+DEAL::dim = 3
+DEAL::
+DEAL::test_zalesak_disk_level_set
+DEAL::center
+DEAL::point = 0.00000 0.00000 0.00000
+DEAL::value = -0.250000
+DEAL::inside
+DEAL::point = 0.100000 0.00000 0.00000
+DEAL::value = -0.250000
+DEAL::point = 0.300000 0.00000 0.00000
+DEAL::value = -0.254951
+DEAL::notch
+DEAL::point = 0.00000 -0.100000 -0.450000
+DEAL::value = 0.200000
+DEAL::outside
+DEAL::point = 0.00000 -0.100000 -2.00000
+DEAL::value = 1.00250
+DEAL::point = -0.300000 -0.300000 -1.20000
+DEAL::value = 0.272792
+DEAL::point = 0.00000 0.00000 -2.00000
+DEAL::value = 1.00000
+DEAL::point = 0.00000 0.00000 2.00000
+DEAL::value = 1.00000
+DEAL::point = -2.00000 0.00000 2.00000
+DEAL::value = 1.82843
+DEAL::point = 2.00000 0.00000 2.00000
+DEAL::value = 1.82843
+DEAL::
projects / dealii.git / commitdiff