Discrepancy in boundary normal for embedded geometries

[diliphridya]

@santiagobadia @fverdugo
The normal vector obtained from the following code is shown in the figure.

using Gridap
using GridapEmbedded
using Gridap.ReferenceFEs
u(x) = x[2]
f(x) = -Δ(u)(x)
ud(x) = u(x)

L = 1.0
domain = (0.0,L,0.0,L)
n = 3
partition = (n,n)
bgmodel = CartesianDiscreteModel(domain,partition)

geo = quadrilateral(x0=Point(0.1,0.1),d1=VectorValue(0.8,0.0),d2=VectorValue(0.0,0.8))
cutgeo = cut(bgmodel,geo)

trian_Γ = EmbeddedBoundary(cutgeo)
n_Γ = get_normal_vector(trian_Γ)
writevtk(trian_Γ,"trian_G",cellfields=["n"=>n_Γ])

Here, the two edges of the quadrilateral are not seen.

On the other hand , if the following geometry is considered,

R = 0.5
L = 0.8*(2*R)
p1 = Point(0.0,0.0)
p2 = p1 + VectorValue(L,0.0)

geo1 = disk(R,x0=p1)
geo2 = disk(R,x0=p2)
geo3 = setdiff(geo1,geo2)

a = 1.01
pmin = p1-a*R
pmax = p1+a*R

bgmodel = CartesianDiscreteModel(pmin,pmax,partition)
cutgeo = cut(bgmodel,geo3)

Then, the normal vector is obtained is :

[fverdugo]

Hi @diliphridya

It seems that there is some problem in the definition of the quadrilateral function. From this line

geo = quadrilateral(x0=Point(0.1,0.1),d1=VectorValue(0.8,0.0),d2=VectorValue(0.0,0.8))

I would expect a square with size = 0.8, but I get this:

In green you see the background mesh and in pink the integration mesh for the generated square. Clearly the pink shape is not a square. Moreover, you can see that part of the boundary (yellow lines) is not contained in the background mesh.

This is the code I have used to visualize all the geometries:

using Gridap
using GridapEmbedded
using Gridap.ReferenceFEs
u(x) = x[2]
f(x) = -Δ(u)(x)
ud(x) = u(x)

L = 1.0
domain = (0.0,L,0.0,L)
n = 3
partition = (n,n)
bgmodel = CartesianDiscreteModel(domain,partition)

bgtrian = Triangulation(bgmodel)
writevtk(bgtrian,"bgtrian")

geo = quadrilateral(x0=Point(0.1,0.1),d1=VectorValue(0.8,0.0),d2=VectorValue(0.0,0.8))
cutgeo = cut(bgmodel,geo)

trian_Ω = Triangulation(cutgeo)
writevtk(trian_Ω,"trian_Ω")

trian_Γ = EmbeddedBoundary(cutgeo)
n_Γ = get_normal_vector(trian_Γ)
writevtk(trian_Γ,"trian_G",cellfields=["n"=>n_Γ])

[fverdugo]

BTW, by using the square function instead of quadrilateral, every thing works OK:

using Gridap
using GridapEmbedded
using Gridap.ReferenceFEs
u(x) = x[2]
f(x) = -Δ(u)(x)
ud(x) = u(x)

L = 1.0
domain = (0.0,L,0.0,L)
n = 3
partition = (n,n)
bgmodel = CartesianDiscreteModel(domain,partition)

bgtrian = Triangulation(bgmodel)
writevtk(bgtrian,"bgtrian")

#geo = quadrilateral(x0=Point(0.1,0.1),d1=VectorValue(0.8,0.0),d2=VectorValue(0.0,0.8))
geo = square(L=0.8,x0=Point(0.5,0.5))
cutgeo = cut(bgmodel,geo)

trian_Ω = Triangulation(cutgeo)
writevtk(trian_Ω,"trian_Ω")

trian_Γ = EmbeddedBoundary(cutgeo)
n_Γ = get_normal_vector(trian_Γ)
writevtk(trian_Γ,"trian_G",cellfields=["n"=>n_Γ])

[fverdugo]

@diliphridya I recommend you to investigate what is wrong with quadrilateral.

[diliphridya]

Thanks @fverdugo . Sure, I will have a look at the quadrilateral function.

[diliphridya]

There was small error in signs in the quadrilateral function. Now, with the fix I can see the same result as from square.

[fverdugo]

Perfect! Please, Open a PR with the fix.

[diliphridya]

Sure.