Add signed distance function for Zalesak's disk #15033

mschreter · 2023-04-06T13:08:12Z

This PR proposes to add a signed distance function for Zalesak's disk (Zalesak, S. T. "Fully multidimensional flux-corrected transport algorithms for fluids." Journal of computational physics 31.3 (1979): 335-362.), which is a common benchmark in the level set community.

bangerth

Since you already have the reference ready, would you mind adding an entry to doc/doxygen/references.bib and then in the documentation of the class cite the paper via @cite ...?

bangerth · 2023-04-06T17:49:28Z

include/deal.II/base/function_signed_distance.h

+    /**
+     * Signed-distance level set function of a Zalesak's disk.
+     *
+     * This function is zero on the disk, negative "inside" and positive


Suggested change

* This function is zero on the disk, negative "inside" and positive

* This function is zero on the surface of the disk, negative "inside" and positive

bangerth · 2023-04-06T17:51:41Z

source/base/function_signed_distance.cc

+                        std::numeric_limits<
+                          double>::lowest() /* notch is open in y-direction*/,


Nice job -- very few people know about the difference between std::numeric_limits::lowest() and std::numeric_limit::min() and either incorrectly use the latter, or give up and use -std::numeric_limits::max() :-)

bangerth · 2023-04-06T17:52:57Z

source/base/function_signed_distance.cc

+      AssertThrow(
+        notch_width <= 2 * radius,
+        ExcMessage(
+          "The width of the notch must be less than the circle diameter. Abort..."));
+      AssertThrow(
+        notch_height <= 2 * radius,
+        ExcMessage(
+          "The height of the notch must be less than the circle diameter. Abort..."));


Let's just use regular Assert here and drop the last half-sentence (our error messages typically describe what is wrong -- a state -- rather than what will happen -- an action):

Suggested change

AssertThrow(

notch_width <= 2 * radius,

ExcMessage(

"The width of the notch must be less than the circle diameter. Abort..."));

AssertThrow(

notch_height <= 2 * radius,

ExcMessage(

"The height of the notch must be less than the circle diameter. Abort..."));

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."));

peterrum

A few comments. I'll continue later.

Could we move the picture to https://github.com/dealii/old-website/tree/main/images? We don't need to fill the code repo with pictures.

peterrum · 2023-04-06T15:27:26Z

include/deal.II/base/function_signed_distance.h

+      Functions::SignedDistance::Sphere<dim> sphere;
+      BoundingBox<dim>                       notch;
+      double                                 radius;
+      Point<dim>                             center_sphere;


peterrum · 2023-04-06T15:27:53Z

include/deal.II/base/function_signed_distance.h

+     * Contour surfaces of the signed distance function of a 3D Zalesak's disk
+     * are illustrated below:
+     *
+     * @image html signed_distance_zalesak_disk.png


I guess the right figure is the same for 2D?

What happens in 1D?

I am not sure, it depends on what Sphere<1> does.

source/base/function_signed_distance.cc

bangerth · 2023-04-06T20:07:58Z

@mschreter If you send me pictures, I'm happy to add them to the repository.

simonsticko · 2023-04-06T20:17:54Z

source/base/function_signed_distance.cc

+      // case: inside disk
+      else if (sphere.value(p) <= 0)
+        {
+          return std::max(sphere.value(p), -notch.signed_distance(p));


I am super-confused about the implementation of the value-function. What you do on this line is a setminus operation with level set functions:

(table from http://doi.wiley.com/10.1002/nme.4823).
But then it looks to me like the whole value-function could be just this line. Why do you need the rest?

Thanks a lot, I will have a look at it asap.

In our case, the setminus operation results in the correct signed distance function for points lying inside the disk but not necessarily for points outside the disk. This is due to the following reasons: notch actually represents the signed distance function to the bounding box of the notch, i.e., the signed distance function to a closed hyperrectangle. Thus, for points lying outside the disk and in the vicinity of the corner points of the notch (2D: zones 6,7,10 and 11 in the drawing below taken from Mousavi, R. (2014). Level set method for simulating the dynamics of the fluid-fluid interfaces: Application of a Discontinuous Galerkin Method)

we need special treatment to get the correct signed distance function. To this end, I introduced the cases in the remaining lines

dealii/source/base/function_signed_distance.cc

Lines 459 to 489 in d332f08

// point is within bounds of the notch faces

if (*std::min_element(signed_distances_notch.begin(),

signed_distances_notch.end()) > 0)

{

const double r =

std::sqrt(std::pow(radius, 2) -

std::pow(notch.side_length(dim - 1) - radius, 2));

// inside notch

if ((dim == 2 && notch.point_inside(p)) ||

(dim == 3 && p.distance(center_sphere) <= radius))

{

return std::min(

distance_front_notch_edge,

std::min(distance_bottom_notch_edge,

*std::min_element(signed_distances_notch.begin(),

signed_distances_notch.end())));

}

// case: below top face of notch

else if (dim == 2 ||

(dim == 3 && (std::abs(p[1] - center_sphere[1])) < r))

{

return std::min(distance_bottom_notch_edge,

signed_distances_notch[idx_top]);

}

// case: corner, if the point is outside sphere but within notch

// bounds

else

return std::min(distance_bottom_notch_edge,

distance_front_notch_edge);

}

and

dealii/source/base/function_signed_distance.cc

Lines 495 to 514 in d332f08

else

{

// consider only points below top face of notch

if (signed_distances_notch[idx_top] > 0)

{

// compute predictor for closest point to sphere.

const auto temp = p + (center_sphere - p) /

(center_sphere.distance(p)) *

sphere.value(p);

// case: notch edge, if the predicted point would be inside the

// notch.

if (notch.signed_distance(temp, 0) <= 0)

return std::min(distance_bottom_notch_edge,

distance_front_notch_edge);

}

// case: sphere

return sphere.value(p);

}

.

I have to admit that the name notch may be a bit misleading here, and maybe that's what led to your confusion.
Since I went through the code lines again and found that they are not very easy to follow unless you are directly involved in the implementation, I will try to come up with a more concise and readable implementation as soon as possible.

I'm afraid I still don't understand why you need the cases. If you just make the notch, so that it extends outside the disk:

I don't think there is any problem with setminus. If you compare the implementation you have now with only using setminus it looks like this:

And I really believe the setminus will give the correct sign-distance function, because its gradient will have norm equal to 1, if the two functions you construct it from have gradients with norm equal to 1:

Convinced. I really could not believe it was that easy, but I am glad it is. Thank you @simonsticko for your efforts. I added your name to the changelog file, I hope that is okay with you as well.

mschreter · 2023-04-13T22:19:49Z

I've tried to restructure the value() function a bit and added some more explanations. We now distinguish 3 cases in 1D, 4 cases in 2D and 6 cases in 3D. That's the best I could come up with for now. If any of you has further suggestions for improvements, I will be happy to incorporate them. Otherwise I am done from my side.

Thanks to the discussion with @simonsticko the value() function is now a one-liner, which is really cool. Sorry for overcomplicating things at first.

peterrum · 2023-04-14T21:17:26Z

include/deal.II/base/function_signed_distance.h

+      /**
+       * Nothch described by a rectangle.
+       */
+      const Functions::SignedDistance::Rectangle<dim> notch;


I have pushed a few minor simplifications. In particular, I have converted the notch to a Function. I think this is useful once we want to extract the level-set merging operations to helper functions.

@mschreter Could you rebase? I think the PR should be ready now! @simonsticko Could you take a second look?

simonsticko · 2023-04-15T13:47:13Z

include/deal.II/base/function_signed_distance.h

+
+
+    /**
+     * Signed-distance level set function of a Zalesak's disk proposed in


I believe this "a" should go away, i.e. "Signed-distance level set function of Zalesak's disk proposed in ...".

Otherwise I think this looks great now. 💯

kronbichler

Looks good to me. Thank you for your insistence @simonsticko, I think this is a very simple code now. Regarding the added picture, we have #15036, so I guess we should have a decision there soon?

peterrum · 2023-04-15T18:34:11Z

Regarding the added picture, we have #15036, so I guess we should have a decision there soon?

For the time being, let's move it to the webpage. I have do that.

@bangerth

The comments by @bangerth have been addressed, so I'm merging now.

mschreter force-pushed the add_signed_distance_funcion_zalesak_disk branch from 596c12a to 7e70130 Compare April 6, 2023 13:14

bangerth previously requested changes Apr 6, 2023

View reviewed changes

peterrum reviewed Apr 6, 2023

View reviewed changes

peterrum mentioned this pull request Apr 6, 2023

Where to put pictures of the documentation? #15036

Closed

simonsticko reviewed Apr 6, 2023

View reviewed changes

drwells added the Enhancement label Apr 7, 2023

peterrum reviewed Apr 14, 2023

View reviewed changes

peterrum added the ready to test label Apr 14, 2023

mschreter force-pushed the add_signed_distance_funcion_zalesak_disk branch from 5d10a22 to 64c6dda Compare April 14, 2023 21:39

simonsticko reviewed Apr 15, 2023

View reviewed changes

kronbichler approved these changes Apr 15, 2023

View reviewed changes

mschreter force-pushed the add_signed_distance_funcion_zalesak_disk branch 2 times, most recently from e9a7e46 to 5a8bbfe Compare April 15, 2023 16:17

add signed distance function for Zalesak's disk

d8eaa2e

peterrum force-pushed the add_signed_distance_funcion_zalesak_disk branch from 5a8bbfe to d8eaa2e Compare April 15, 2023 18:33

kronbichler merged commit 464e8e0 into dealii:master Apr 17, 2023
14 checks passed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add signed distance function for Zalesak's disk #15033

Add signed distance function for Zalesak's disk #15033

mschreter commented Apr 6, 2023 •

edited by peterrum

bangerth left a comment

bangerth Apr 6, 2023

bangerth Apr 6, 2023

bangerth Apr 6, 2023

peterrum left a comment

peterrum Apr 6, 2023

peterrum Apr 6, 2023

peterrum Apr 6, 2023

mschreter Apr 6, 2023

mschreter Apr 13, 2023

bangerth commented Apr 6, 2023

simonsticko Apr 6, 2023

mschreter Apr 6, 2023

mschreter Apr 11, 2023

simonsticko Apr 14, 2023

mschreter Apr 14, 2023 •

edited

mschreter commented Apr 13, 2023 •

edited

peterrum Apr 14, 2023

simonsticko Apr 15, 2023

kronbichler left a comment

peterrum commented Apr 15, 2023

	* This function is zero on the disk, negative "inside" and positive
	* This function is zero on the surface of the disk, negative "inside" and positive

		std::numeric_limits<
		double>::lowest() /* notch is open in y-direction*/,

	// point is within bounds of the notch faces
	if (*std::min_element(signed_distances_notch.begin(),
	signed_distances_notch.end()) > 0)
	{
	const double r =
	std::sqrt(std::pow(radius, 2) -
	std::pow(notch.side_length(dim - 1) - radius, 2));

	// inside notch
	if ((dim == 2 && notch.point_inside(p)) \|\|
	(dim == 3 && p.distance(center_sphere) <= radius))
	{
	return std::min(
	distance_front_notch_edge,
	std::min(distance_bottom_notch_edge,
	*std::min_element(signed_distances_notch.begin(),
	signed_distances_notch.end())));
	}
	// case: below top face of notch
	else if (dim == 2 \|\|
	(dim == 3 && (std::abs(p[1] - center_sphere[1])) < r))
	{
	return std::min(distance_bottom_notch_edge,
	signed_distances_notch[idx_top]);
	}
	// case: corner, if the point is outside sphere but within notch
	// bounds
	else
	return std::min(distance_bottom_notch_edge,
	distance_front_notch_edge);
	}

	else
	{
	// consider only points below top face of notch
	if (signed_distances_notch[idx_top] > 0)
	{
	// compute predictor for closest point to sphere.
	const auto temp = p + (center_sphere - p) /
	(center_sphere.distance(p)) *
	sphere.value(p);

	// case: notch edge, if the predicted point would be inside the
	// notch.
	if (notch.signed_distance(temp, 0) <= 0)
	return std::min(distance_bottom_notch_edge,
	distance_front_notch_edge);
	}

	// case: sphere
	return sphere.value(p);
	}



		/**
		* Signed-distance level set function of a Zalesak's disk proposed in

Add signed distance function for Zalesak's disk #15033

Add signed distance function for Zalesak's disk #15033

Conversation

mschreter commented Apr 6, 2023 • edited by peterrum

bangerth left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

peterrum left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bangerth commented Apr 6, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mschreter Apr 14, 2023 • edited

Choose a reason for hiding this comment

mschreter commented Apr 13, 2023 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kronbichler left a comment

Choose a reason for hiding this comment

peterrum commented Apr 15, 2023

mschreter commented Apr 6, 2023 •

edited by peterrum

mschreter Apr 14, 2023 •

edited

mschreter commented Apr 13, 2023 •

edited