Open
Description
Find below a MWE which reproduces a (very) nasty BUG we have found in the implementation of change_domain among SkeletonTriangulations. The BUG exposes itself whenever one wants to change_domain among SkeletonTriangulations for which the plus and minus side is swapped. change_domain generates a CellField of the proper Triangultaion, BUT, it does not swap the plus and minus sides of the corresponding arrays.
using Gridap
using Test
domain = (0.0, 1.0, 0.0, 1.0)
partition = (1,2)
model = CartesianDiscreteModel(domain,partition)
order=1
reffe_u = ReferenceFE(lagrangian,Float64,order)
U = FESpace(model,reffe_u,conformity=:L2)
v = get_fe_basis(U)
num_facets = num_faces(model,1)
face_to_mask = Vector{Bool}(undef,num_facets)
face_to_mask .= false
face_to_mask[2] = true
left_cell_around = 1
plus = BoundaryTriangulation(model,face_to_mask,left_cell_around)
right_cell_around = 2
minus = BoundaryTriangulation(model,face_to_mask,right_cell_around)
Λ1 = SkeletonTriangulation(plus,minus)
Λ2 = SkeletonTriangulation(minus,plus)
degree = 0
dΛ1 = Measure(Λ1,degree)
dΛ2 = Measure(Λ2,degree)
x1=get_cell_points(dΛ1.quad)
x2=get_cell_points(dΛ2.quad)
jv=jump(v)
jv1 = Gridap.CellData.change_domain(jv,
get_triangulation(x1),
Gridap.CellData.DomainStyle(jv))
jv2 = Gridap.CellData.change_domain(jv,
get_triangulation(x2),
Gridap.CellData.DomainStyle(jv))
jv1_to_jv2 = Gridap.CellData.change_domain(jv1,
get_triangulation(x2),
Gridap.CellData.DomainStyle(jv))
jv1_to_jv2_at_x2=jv1_to_jv2(x2)
jv2_at_x2=jv2(x2)
@test_broken norm(jv1_to_jv2_at_x2[1][1]-jv2_at_x2[1][1]) ≈ 0.0 # MUST BE TRUE AND ITS FALSE
@test_broken norm(jv1_to_jv2_at_x2[1][2]-jv2_at_x2[1][2]) ≈ 0.0 # MUST BE TRUE AND ITS TRUE
@test_broken !(norm(jv1_to_jv2_at_x2[1][1]+jv2_at_x2[1][2]) ≈ 0.0) # MUST BE FALSE AND ITS TRUE
@test_broken !(norm(jv1_to_jv2_at_x2[1][2]+jv2_at_x2[1][1]) ≈ 0.0) # MUST BE FALSE AND ITS TRUE