7
$\begingroup$

I'm trying to understand a Monte-Carlo Laplace/Poisson PDE solver: http://www.cs.cmu.edu/~kmcrane/Projects/MonteCarloGeometryProcessing/paper.pdf

This method inspired by random-walks and ray-tracing has an advantage over FEM mesh methods in that it requires little to no mesh pre-processing. I do not expect great performance from Mathematica but that's okay - I'm just trying to play with concepts in the paper.

I want to perform a walk on spheres in Mathematica shown in this diagram in the paper:

walk on sphere

Questions:

  1. Given a mesh such as mesh = DiscretizeGraphics@ExampleData[{"Geometry3D", "Cow"}] how do I get a random point $x_0$ on the interior of the mesh $\Omega$ ? I've tried RandomPoint but it only gives points on the surface.
  2. How do I find the closest point on the mesh surface $\partial\Omega$ to $x_0$ so I can then determine the radius $r_0$ of the first sphere?

Once I've got $x_0, r_0$ I can do RandomPoint on the sphere to generate $x_1$ and Nest this process until the sphere radius is below some threshold.

Improve this question
$\endgroup$
7
  • $\begingroup$ Replace DiscretizeGraphics by BoundaryDiscretizeGraphics and use RegionNearest. RegionNearest[mesh] generates a RegionNearestFunction that can (and should) be reused. $\endgroup$
    – Henrik Schumacher
    May 26, 2020 at 20:31
  • $\begingroup$ If I do that I get "The first argument BoundaryDiscretizeGraphics[] is expected to be parameter-free" when calling RandomPoint $\endgroup$
    – flinty
    May 26, 2020 at 20:34
  • 2
    $\begingroup$ Ah, the cow is not watertight. You may try ExampleData[{"Geometry3D", "Triceratops"}, "BoundaryMeshRegion"] instead (it is already a BoundaryMeshRegion has `RegionDimension equal to 3). $\endgroup$
    – Henrik Schumacher
    May 26, 2020 at 20:39
  • $\begingroup$ With the above and rnf=RegionNearest[mesh], the point x0=RandomPoint[mesh] is correctly on the interior, but the nearest point I expected to lie on the surface rnf[x0] is actually the same point x0 $\endgroup$
    – flinty
    May 26, 2020 at 21:02
  • 2
    $\begingroup$ @flinty If Henrik doesn't write it up as an answer, then perhaps you could. As you probably know, self-answers are encouraged! $\endgroup$
    – MarcoB
    May 26, 2020 at 21:45

1 Answer 1

Reset to default
9
$\begingroup$

Thanks to Henrik Schumacher I got this to work with NestWhileList and visualized the walk on spheres. The other outputs of nextPoint besides RandomPoint are for visualization and the termination of the NestWhileList when the radius is small enough:

mesh = ExampleData[{"Geometry3D", "Triceratops"}, "BoundaryMeshRegion"];
rnf = RegionNearest@RegionBoundary@mesh;

nextPoint[p_] := Block[{r = EuclideanDistance[rnf[p], p]},
 {RandomPoint[Sphere[p, r]], p, r}]

walk = NestWhileList[
   nextPoint[#[[1]]] &, {RandomPoint[mesh], {}, ∞}, #[[3]] > 10^-3 &];

Graphics3D[{Opacity[.0], mesh,
  Opacity[1], Red, Thick, Line[walk[[All, 1]]],
  Blue, Opacity[0.04], 
  Sphere[#[[2]], #[[3]]] & /@ Rest[walk]}]

walk on spheres triceratops mesh

And with simple modifications it works for 2D, here demonstrated with a random polygon:

reg = RandomPolygon[12];
rnf = RegionNearest@RegionBoundary@reg;

nextPoint[p_] := 
 Block[{r = EuclideanDistance[rnf[p], p]}, {RandomPoint[Circle[p, r]], p, r}]

walk = NestWhileList[
   nextPoint[#[[1]]] &, {RandomPoint[reg], {}, ∞}, #[[3]] > 10^-5 &];

Graphics[{Opacity[.1], reg, Opacity[1], Red, Thick, 
  Line[walk[[All, 1]]], Blue, Opacity[0.04], 
  Disk[#[[2]], #[[3]]] & /@ Rest[walk]}]

walk on circles random polygon 2d

Improve this answer
$\endgroup$
1
  • 3
    $\begingroup$ Love the paper, thanks for bringing it up in a Mathematica context! Would be awesome to see the whole poisson/laplace-solver working in Mathematica in case you are implementing it. $\endgroup$
    – Thies Heidecke
    May 26, 2020 at 22:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.