Error integrating on embedded interface

[oriolcg]

trafficstars

Hi @fverdugo ,

The following code reproduces the error that I get when integrating over embedded interfaces:

using Gridap
using GridapEmbedded

𝒯 = CartesianDiscreteModel((0,1,0,1),(8,8)) |> simplexify

# Embedded geometry
R = 0.25
C = Point(0.5,0.5)
solid_geo = disk(R,x0=C,name="solid")
fluid_geo = !(solid_geo,name="fluid")
𝒯c = cut(𝒯,union(fluid_geo,solid_geo))

# Triangulations
Ω = Interior(𝒯)
Ωc = Interior(𝒯c)
Ωs_act = Triangulation(𝒯c,ACTIVE,solid_geo)
Ωf_act = Triangulation(𝒯c,ACTIVE,fluid_geo)
Ωs = Triangulation(𝒯c,PHYSICAL,solid_geo)
Ωf = Triangulation(𝒯c,PHYSICAL,fluid_geo)
Γi = EmbeddedBoundary(𝒯c,fluid_geo,solid_geo)
dΓi = Measure(Γi,2)

# AgFEM
strategy = AggregateAllCutCells()
agg_f = aggregate(strategy,𝒯c,fluid_geo)
agg_s = aggregate(strategy,𝒯c,solid_geo)

model_fluid = get_active_model(Ωf_act)
model_solid = get_active_model(Ωs_act)
order = 2
FEᵥ_f_std = FiniteElements(
  PhysicalDomain(),
  model_fluid,
  lagrangian,
  VectorValue{2,Float64},
  order
)
FEᵥ_s_std = FiniteElements(
  PhysicalDomain(),
  model_solid,
  lagrangian,
  VectorValue{2,Float64},
  order
)
FEᵥ_f_ser = FiniteElements(
  PhysicalDomain(),
  model_fluid,
  lagrangian,
  VectorValue{2,Float64},
  order,
  space=:S,
  conformity=:L2
)
FEᵥ_s_ser = FiniteElements(
  PhysicalDomain(),
  model_solid,
  lagrangian,
  VectorValue{2,Float64},
  order,
  space=:S,
  conformity=:L2
)

Vᵥ_f_std = FESpace(Ωf_act,FEᵥ_f_std)
Vᵥ_s_std = FESpace(Ωs_act,FEᵥ_s_std)
Vᵥ_f_ser = FESpace(Ωf_act,FEᵥ_f_ser)
Vᵥ_s_ser = FESpace(Ωs_act,FEᵥ_s_ser)
Vᵥ_f = AgFEMSpace(Vᵥ_f_std,agg_f,Vᵥ_f_ser)
Vᵥ_s = AgFEMSpace(Vᵥ_s_std,agg_s,Vᵥ_s_ser)
Uᵥ_f = TrialFESpace(Vᵥ_f)
Uᵥ_s = TrialFESpace(Vᵥ_s)

uf = FEFunction(Uᵥ_f,rand(num_free_dofs(Uᵥ_f)))
us = FEFunction(Uᵥ_s,rand(num_free_dofs(Uᵥ_s)))

val = ∑(∫(uf-us)dΓi)

[fverdugo]

@oriolcg

the code below works. Two remarks:

• Use interpolations in the reference space. We have encountered some other issues in the past when combining AgFEM + interpolations in the physical space.
• The serendipity extension is only needed for the fluid space.

using Gridap
using GridapEmbedded

𝒯 = CartesianDiscreteModel((0,1,0,1),(8,8)) |> simplexify

# Embedded geometry
R = 0.25
C = Point(0.5,0.5)
solid_geo = disk(R,x0=C,name="solid")
fluid_geo = !(solid_geo,name="fluid")
𝒯c = cut(𝒯,union(fluid_geo,solid_geo))

# Triangulations
Ω = Interior(𝒯)
Ωc = Interior(𝒯c)
Ωs_act = Triangulation(𝒯c,ACTIVE,solid_geo)
Ωf_act = Triangulation(𝒯c,ACTIVE,fluid_geo)
Ωs = Triangulation(𝒯c,PHYSICAL,solid_geo)
Ωf = Triangulation(𝒯c,PHYSICAL,fluid_geo)
Γi = EmbeddedBoundary(𝒯c,fluid_geo,solid_geo)
dΓi = Measure(Γi,2)

# AgFEM
strategy = AggregateAllCutCells()
agg_f = aggregate(strategy,𝒯c,fluid_geo)
agg_s = aggregate(strategy,𝒯c,solid_geo)

Tv = VectorValue{2,Float64}
k = 2
reffe_std = ReferenceFE(lagrangian,Tv,k)
reffe_ser = ReferenceFE(lagrangian,Tv,k,space=:S)

Vf_std = FESpace(Ωf_act,reffe_std)
Vs_std = FESpace(Ωs_act,reffe_std)
Vf_ser = FESpace(Ωf_act,reffe_ser)

Vf = AgFEMSpace(Vf_std,agg_f,Vf_ser)
Vs = AgFEMSpace(Vs_std,agg_s)

Uf = TrialFESpace(Vf)
Us = TrialFESpace(Vs)

uf = FEFunction(Uf,rand(num_free_dofs(Uf)))
us = FEFunction(Us,rand(num_free_dofs(Us)))

val = ∑(∫(uf-us)dΓi)