Embedded Model with Arbitrary Level-Set Function

[MichaelLoibl]

Dear all,
I am trying to set up an embedded problem using GridapEmbedded.jl. My level-set functions are in general arbitrary functions and are not derived from standard shapes. In this context, I would need a generic description of a level-set function. I only saw embedded examples where basic geometries like disks, cylinders, spheres or cubes and boolean operations of them are supported.

Is it possible to have a generic function in GridapEmbedded?

If yes, could somebody provide me an example? E.g. for a triangle (even though that would be maybe supported by standard geometries as well) where the level-set functions are:

f_1(x,y) = - (2^(1/2)*x)/2 - (2^(1/2)*(y - 9/10))/2
f_2(x,y) = x - 9/10
f_3(x,y) = y - 9/10

Any other example is also welcomed.

Best regards,
Michael Loibl

[JordiManyer]

https://github.com/gridap/GridapEmbedded.jl/blob/5d31c824215c0bb3d2ec18d949a665ce956a1be0/test/LevelSetCuttersTests/DiscreteGeometriesTests.jl#L49

[MichaelLoibl]

Thank you. I was not looking at the test folder, but only in the examples folder and the tutorial. Thus, I oversaw it.

However, when I run it, it actually fails on that step with the arbitrary function:

julia "d:\Code\GridapEmbedded\test\LevelSetCuttersTests\DiscreteGeometriesTests.jl"
ERROR: LoadError: MethodError: no method matching GridapEmbedded.LevelSetCutters.DiscreteGeometry(::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{Gridap.CellData.ReferenceDomain}}, ::Gridap.Geometry.CartesianDiscreteModel{2, Float64, typeof(identity)})
The type `GridapEmbedded.LevelSetCutters.DiscreteGeometry` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  GridapEmbedded.LevelSetCutters.DiscreteGeometry(!Matched::GridapEmbedded.CSG.Node, !Matched::Array{Gridap.TensorValues.VectorValue{D, T}, 1}) where {D, T}
   @ GridapEmbedded C:\Users\Loibl\.julia\packages\GridapEmbedded\i6MoZ\src\LevelSetCutters\DiscreteGeometries.jl:3
  GridapEmbedded.LevelSetCutters.DiscreteGeometry(!Matched::AbstractVector, !Matched::AbstractVector; name)
   @ GridapEmbedded C:\Users\Loibl\.julia\packages\GridapEmbedded\i6MoZ\src\LevelSetCutters\DiscreteGeometries.jl:68

Stacktrace:
 [1] top-level scope
   @ d:\Code\GridapEmbedded\test\LevelSetCuttersTests\DiscreteGeometriesTests.jl:50
in expression starting at d:\Code\GridapEmbedded\test\LevelSetCuttersTests\DiscreteGeometriesTests.jl:1

I also run the test:

Test Summary:          | Pass  Error  Total     Time
LevelSetCutters        |  103      1    104  1m12.7s
  AnalyticalGeometries |   95            95    59.8s
  DiscreteGeometries   |    4      1      5     7.1s
  LookupTables         |    4             4     2.0s
  CutTriangulations    |                  0     1.1s
  LevelSetCutters      |                  0     2.7s

The DiscreteGeometriesTests is the only one failing for me.

Any idea what is going wrong?

[JordiManyer]

After looking at Git Blame, I believe those were added recently and I haven't done the release yet. You can point to the master branch for now, I want to hold off the release a little longer.

From your project, run

pkg > add GridapEmbedded#master

[MichaelLoibl]

Thank you.

My actual goal is to have a comparison for integration of cut elements. I think that your code could help me hereby. It would be perfect to compute the integral of an arbitrary function over the domain. It would be even better to get the corresponding quadrature points then. However, when I try to compute the area of the cut domain , I have problems to actually evaluate the integral.

This is my code:

module example_circle_1

using AbstractTrees
using Gridap
using Gridap.ReferenceFEs
using Gridap.Arrays
using GridapEmbedded
using GridapEmbedded.CSG
using GridapEmbedded.LevelSetCutters

#using a cellfield
domain = (0,1,0,1)
n = 10
bgmodel = CartesianDiscreteModel(domain,(n,n))

reffe = ReferenceFE(lagrangian,Float64,1)
Ω_bg = Triangulation(bgmodel)
V_bg = FESpace(Ω_bg,reffe)
φh = interpolate(x->(x[1] - 1/2)^2 + (x[2] - 1/2)^2 - 1/25,V_bg)
geo = DiscreteGeometry(φh,bgmodel)

# Cut the background model
cutgeo = cut(bgmodel,geo)

Ω_act = Triangulation(cutgeo,ACTIVE)

writevtk(Ω_bg,"bg_trian")
writevtk(Ω_act,"act_trian")

Ω = Triangulation(cutgeo,PHYSICAL)
writevtk(Ω,"phys_trian")

order = 1
degree = 2*order
dΩ = Measure(Ω,degree)

area = ∫(1.0)*dΩ

end # module

What would be the correct way to evaluate the integral (and maybe also get the quadrature points)?

Best,
Michael

[JordiManyer]

That works for me just fine. I don't know what error you are talking about. Make sure you have updated libraries (have a look at the versions).

The value of the area is given by collecting the DomainContribution. Just use area = sum(∫(1.0)*dΩ). The quadrature points (in the reference space) are given by pts = Gridap.CellData.get_data(get_cell_point(dΩ)).

The library is based on Gridap, so I woudl suggest having a look at the tutorials.