Interpolable for non-conforming (octree) meshes

[ericneiva]

Hello @JordiManyer and @amartinhuertas,

For non-conforming meshes, I think we should find a way in the function below to provide a list of vertex_to_cells that includes the non-conforming connectivities:

Gridap.jl/src/CellData/Interpolation.jl

Lines 99 to 118 in 62586b6

	function _point_to_cell_cache(searchmethod::KDTreeSearch,trian::Triangulation)
	model = get_active_model(trian)
	topo = get_grid_topology(model)
	if num_nodes(model) == num_vertices(model)
	# Non-periodic case
	vertex_coordinates = Geometry.get_vertex_coordinates(topo)
	vertex_to_cells = get_faces(topo, 0, num_cell_dims(trian))
	else
	# Periodic case
	vertex_coordinates = collect1d(Geometry.get_node_coordinates(model))
	cell_to_vertices = Table(get_cell_node_ids(model))
	vertex_to_cells = Arrays.inverse_table(cell_to_vertices, num_nodes(model))
	end
	kdtree = KDTree(map(nc -> SVector(Tuple(nc)), vertex_coordinates))
	cell_to_ctype = get_cell_type(trian)
	ctype_to_polytope = get_polytopes(trian)
	cell_map = get_cell_map(trian)
	table_cache = array_cache(vertex_to_cells)
	return searchmethod, kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map, table_cache
	end

For OctreeDistributedDiscreteModel, one possibility could be to use the NonConformingGridTopolgy that I am developing here.

Otherwise, I need to increase the num_nearest_vertices to the minimum value that ensures there is always (at least one) conforming node in the list of nearest vertices. This is not ideal.

Here is a MWE that illustrates the issue (running it with 1 proc).

using Gridap
using GridapDistributed
using GridapP4est
using PartitionedArrays
using MPI

using GridapEmbedded.AlgoimUtils

if !MPI.Initialized()
  MPI.Init()
end

# mpiexec -n 1 julia +1.11 -O1 --project test/dev/gridap_issue_1266.jl

with_mpi() do distribute

  # Initial Octree Distributed Discrete Model
  ranks        = distribute(LinearIndices((MPI.Comm_size(MPI.COMM_WORLD),)))
  coarse_model = CartesianDiscreteModel((0,1,0,1),(1,1))
  num_levels_initial_refinement = 2
  dmodel = OctreeDistributedDiscreteModel(ranks,
                                          coarse_model,
                                          num_levels_initial_refinement)

  # Coarsen first four cells
  fmodel_refine_coarsen_flags = 
    map(ranks,partition(get_cell_gids(dmodel.dmodel))) do rank,indices
      flags = zeros(Int,length(indices))
      flags[1:4] .=  coarsen_flag
      flags
  end
  fmodel,_ = Gridap.Adaptivity.adapt(dmodel,fmodel_refine_coarsen_flags);
  # writevtk(fmodel,"fmodel")

  # Interpolate a fun on a point of the non-conforming grid
  Ω = Triangulation(fmodel)
  order = 1
  reffe_s = ReferenceFE(lagrangian,Float64,order)
  Qₕ = TestFESpace(Ω,reffe_s,conformity=:H1)
  φ = zero(Qₕ)

  sm = Gridap.CellData.KDTreeSearch(num_nearest_vertices=5)
  iφ = map(local_views(φ)) do liφ
    Gridap.CellData.Interpolable(liφ,searchmethod=sm)
  end |> GridapDistributed.DistributedInterpolable
  iφ(Point(0.25,0.25)) # OK

  φ(Point(0.25,0.25)) # KO
  # Using default num_nearest_vertices = 1 fails because the only candidate
  # vertex is a hanging node and does not have in its list of cells the coarse cell 
  # that contains the point. _Detail:_ To be precise, the point that is not found is 
  # `mean(testitem(get_cell_coordinates(trian)))` of the DistributedInterpolable constructor.

  true
end

# MPI.Finalize()

I also attach a picture of the mesh in the unit square.

FYI @nannaberre

We can talk about this next time we meet. If you have any questions, lmk. Thanks.

[ericneiva]

On a first assessment of this, we agreed that the way to go is to make vertex_to_cells or the topology and optional argument and overload _point_to_cell_cache.