Incorrect Assertion on boundary faces integral

[Paulms]

I want to integrate on a subdomain of the boundary. However I get the following error message:

ERROR: AssertionError:

You are trying to build a CellField from an array of length 9
on a Triangulation with 18 cells. The length of the given array
and the number of cells should match.

It seems to work if the subdomain consists only on one side of the square, but not if it has more than one.

In the following code, I tried to isolate the problem:

using Gridap

# build mesh
domain = (0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition) |> simplexify

# define FE Spaces
Vh = FESpace(model, ReferenceFE(lagrangian,VectorValue{2,Float64},1),
        conformity=:H1)
Qh = FESpace(model, ReferenceFE(lagrangian,Float64,1),
        conformity=:L2, constraint=:zeromean)

Uh = TrialFESpace(Vh)
Ph = TrialFESpace(Qh)

Yh = MultiFieldFESpace([Vh, Qh])
Xh = MultiFieldFESpace([Uh, Ph])

# Discrete domain, quadratures and normals
# 5 = bottom, 7 = left, 8 = right, 6 = top
Ω = Triangulation(model)
ΓSlip = BoundaryTriangulation(model,tags=["tag_5"])
ΓDirichlet = BoundaryTriangulation(model,tags=["tag_6","tag_7","tag_8"])

hFD = integrate(1,CellQuadrature(ΓDirichlet,0))
hFS = integrate(1,CellQuadrature(ΓSlip,0))

dΩ = Measure(Ω,2)
dΓS = Measure(ΓSlip,2)
dΓD = Measure(ΓDirichlet,2)

n_ΓD = get_normal_vector(ΓDirichlet)
n_ΓS = get_normal_vector(ΓSlip)

f = VectorValue(0,0)

# Weak Form
length(hFD) == num_cells(ΓDirichlet) # 9 elements for both sides
a_Ω(u,v) = ∫(ε(v)⊙ε(u))dΩ
# This line errors:
a_ΓD(u,v) =  ∫( 1.0./hFD*(u⋅v))dΓD
l_Ω(v) = ∫( v⋅f )dΩ

B_S((u,p),(v,q)) = a_Ω(u,v) + a_ΓD(u,v)
L_S((v,q)) = l_Ω(v)

op = AffineFEOperator(B_S,L_S,Xh,Yh)

# But this one works fine:
length(hFS) == num_cells(ΓSlip)   # 3 elements for both sides
a_ΓS(u,v) =  ∫( 1.0./hFS*((u⋅n_ΓS)*(v⋅n_ΓS)))dΓS
B_S((u,p),(v,q)) = a_Ω(u,v) + a_ΓS(u,v)
L_S((v,q)) = l_Ω(v)

op = AffineFEOperator(B_S,L_S,Xh,Yh)

Environment:

(@v1.10) pkg> st
Status `~/.julia/environments/v1.10/Project.toml`
  [a93c6f00] DataFrames v1.6.1
  [56d4f2e9] Gridap v0.17.23
  [ade2ca70] Dates
  [56ddb016] Logging

[JordiManyer]

I believe this is happening because before integrating on dΓS it's trying to put everything on the same triangulation (and it might be choosing Ω). The problem is you are not telling Gridap where hFS is defined (beacuse it's an array). What you need is a CellField, which is an array associated to a certain triangulation.

Try creating your edge measures like this:

hFS = CellField(integrate(1,CellQuadrature(ΓSlip,0)),ΓSlip,ReferenceDomain())

That should give Gridap enough information to do the proper change of triangulation.

[Paulms]

Thanks a lot! Your solution fixed the issue perfectly!

[Paulms]

Closing the issue