search

gridap / GridapEmbedded.jl

Embedded finite element methods in Julia
Other
42 stars 14 forks source link

Bug triangulations #98

Open nichomueller opened 1 week ago

nichomueller commented 1 week ago

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 commented 2 days ago

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 commented 2 days ago

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