Assertion Error when domain is aligned with background mesh

[diliphridya]

@fverdugo @santiagobadia
When solving a heat equation using AgFEM (aggregation only in space) , CG(1) discretization in space and DG(0) discretization in time, if the domain is aligned to the Background mesh then an error, AssertionError: n_ldofs == size(mat, 1) is observed in line op = AffineFEOperator(U,V,t_Ω,t_Γ,t_s) in the following code.

using Gridap
using GridapEmbedded
using Gridap.ReferenceFEs
using Gridap: ∇
using Gridap.Visualization
using GridapEmbedded.LevelSetCutters

# Setting directory for storing output files
d = "./"

# Manufactured Solution
u(x) = x[1]^2*x[2]
const k1 = TensorValue(1.0,0.0,0.0,0.0)
const k2 = VectorValue(0.0,1.0)
f(x) = x[1]^2 - 2*x[2]

# Background Cartesian mesh
domain = (0,1,0,1)
n_x = 10
n_t = 10
partition = (n_x,n_t)
bgmodel = CartesianDiscreteModel(domain,partition)
fi = joinpath(d,"bgmodel")
writevtk(bgmodel,fi)

# Domain
geo = quadrilateral(x0=Point(0,0),d1=VectorValue(1,0),d2=VectorValue(0,1))
cutgeo = cut(bgmodel,geo)
model = DiscreteModel(cutgeo)
fi = joinpath(d,"model")
writevtk(model,fi)

# AgFEM strategy
strategy = AggregateSpaceCutCells()
aggregates = aggregatespace(strategy,cutgeo)

# Trial & Test Space
T = Float64
order = (1,0)
reffe = LagrangianRefFE(T,QUAD,order)
conf = CDConformity((CONT,DISC))
face_own_dofs = get_face_own_dofs(reffe,conf)
V0 = TestFESpace(reffe=reffe, model=model, conformity=conf,dof_space=:physical)
V = AgFEMSpace(V0,aggregates)
U = TrialFESpace(V)

# Triangulation
trian_Ω = Triangulation(cutgeo)
trian_Γ = EmbeddedBoundary(cutgeo)
strian = SkeletonTriangulation(model,reffe,face_own_dofs)
n_Γ = get_normal_vector(trian_Γ)
ns = get_normal_vector(strian)
fi = joinpath(d,"strian")
writevtk(strian,fi,cellfields=["normal"=>ns])

# Quadrature
degree = 2 * max(order...)
quad_Ω = CellQuadrature(trian_Ω,degree)
quad_Γ = CellQuadrature(trian_Γ,degree)
squad = CellQuadrature(strian,degree)

# Weak formulation
a_Ω(u,v) = ∇(v)⊙(k1⋅∇(u)) + v*(k2⋅∇(u))
b_Ω(v) = v*f
t_Ω = AffineFETerm(a_Ω,b_Ω,trian_Ω,quad_Ω)

a_s(u,v) = jump(u)*(v.outward)*(-1.0)
t_s = LinearFETerm(a_s,strian,squad)

const γd = 10.0
h = 1/n_x
n_γ = k1⋅n_Γ
a_Γ(u,v) = (γd/h)*v*u  - v*(n_γ⋅(k1⋅∇(u))) - (n_γ⋅(k1⋅∇(v)))*u
l_Γ(v) = (γd/h)*v*u - (n_γ⋅(k1⋅∇(v)))*u
t_Γ = AffineFETerm(a_Γ,l_Γ,trian_Γ,quad_Γ)

# Solving
op = AffineFEOperator(U,V,t_Ω,t_Γ,t_s)
uh = solve(op)
uh_Ω = restrict(uh,trian_Ω)
fi = joinpath(d,"results")
writevtk(trian_Ω,fi,cellfields=["uh"=>uh_Ω])

# Error
e = u-uh_Ω
fi = joinpath(d,"error")
writevtk(trian_Ω,fi,cellfields=["e" => e])

# L2 Error
l2(u) = u*u
el2 = sqrt(sum(integrate(l2(e),trian_Ω,quad_Ω)))
tol = 1.0e-10
@assert el2 < tol

trian = Triangulation(bgmodel)
colors = color_aggregates(aggregates,bgmodel)
writevtk(trian,"trian",celldata=["aggregate"=>aggregates,"color"=>colors],cellfields=["uh"=>uh])

[fverdugo]

@diliphridya Thanks for reporting. I'll take a look.