Inaccurate cuts with disk

[oriolcg]

I'm getting inaccurate cuts when intersecting with a disk geometry, see attached figure. You can see that only the intersection with vertical edges are well captured, while the intersection with horizontal edges seems to be following another geometry.

Did anyone faced similar issues?

[ericneiva]

Hi, Oriol! I have not faced similar issues, because I have always used uniform Cartesian meshes. But as soon as I prescribe an anisotropic Cartesian mesh, I reproduce the behaviour you have identified.

using Gridap
using GridapEmbedded

const R = 0.42

geom = disk(R,x0=Point(0.5,0.5))
partition = (9,3)
domain = (0,1,0,1)

bgmodel = CartesianDiscreteModel(domain,partition)
cutgeo = cut(bgmodel,geom)
Ω = Triangulation(cutgeo)
writevtk(Ω,"trian")

[fverdugo]

Hi @oriolcg,

I am not sure that this is a bug. Marching cubes has some geometrical error (it considers piecewise linear approximations), and perhaps you are seen this error.

can you try with a levelset with a piecewise linear 0-isovalue ? like quadrilateral ?

GridapEmbedded.jl/src/LevelSetCutters/AnalyticalGeometries.jl

Line 135 in c552a14

function quadrilateral(;x0=Point(0,0),d1=VectorValue(1,0),d2=VectorValue(0,1),name="quadrilateral")

[oriolcg]

Thanks for the answer @fverdugo, I'll try with quadrilateral and let you know.

How is the Marching cubes implemented? I would expect that the intersection between element edges and an analytical geometry would be exact. It's true that Marching cubes does a linear approximation, but that's the reconstructed element and it shouldn't affect the accuracy of the intersection point.

[oriolcg]

Indeed, quadrilateral seems to work well. But I'm still not convinced about the Marching cubes argument. I'll try to investigate a bit more.

[oriolcg]

I think I see where is the problem, and it is not related to the Marching cubes algorithm.

I thought that the intersections were obtained from intersecting the background grid edges with the analytical geometry, but they aren't. They are obtained from the levelset function interpolated from the levelset values at the background grid points. For the mesh that I was testing this approach results in big errors because the levelset gives the closest point distance (orthogonal to the circle), while one would be interested on the distance following the edge direction.

I presume that this issue will appear for all geometries that are not linear. The errors can be large for background grids in which at least one of the directions does not capture well the geometry (even if all the other directions are fine enough).

Does it make sense?

A solution for this problem would be to find the actual intersection between each edge and the analytical geometry. But this is quite involved and may be more expensive.