Bug triangulations

[nichomueller]

Found a bug when computing the view of an appended triangulation (SubCellTriangulation + BodyFittedTriangulation).

using Gridap
using Gridap.Arrays
using Gridap.Geometry
using GridapEmbedded

D = 2

model = CartesianDiscreteModel((0,1,0,1),(10,10))

geo = disk(0.1,x0=Point(0.5,0.5))
cutgeo = cut(model,geo)

Ω = Triangulation(cutgeo,PHYSICAL_OUT)
Ω_view = view(Ω,[num_cells(Ω)])

glue = get_glue(Ω_view,Val(D))
tface_to_mface_map = glue.tface_to_mface_map

c = array_cache(tface_to_mface_map)
getindex!(c,tface_to_mface_map,1) # crashes

# source of problem
parent_glue = get_glue(Ω_view.parent,Val(D))
ids = Ω_view.cell_to_parent_cell
tface_to_mface_map = lazy_map(Reindex(parent_glue.tface_to_mface_map),ids) # wrong eltype

[JordiManyer]

Just as a follow-up: The issue has to do with a type instability on the cellmaps coming from the two triangulations that are being appended.

[nichomueller]

Here is a possible fix:

function Base.view(trian::Geometry.AppendedTriangulation,ids::AbstractArray)
  Ti = eltype(ids)
  ids1 = Ti[]
  ids2 = Ti[]
  n1 = num_cells(trian.a)
  for i in ids
    if i <= n1
      push!(ids1,i)
    else
      push!(ids2,i-n1)
    end
  end
  trian1 = view(trian.a,ids1)
  trian2 = view(trian.b,ids2)
  num_cells(trian1) == 0 ? trian2 : lazy_append(trian1,trian2)
end