Fix: Grid(ReferenceFE{d}, model::MappedDiscreteModel) returning cell grid for all dimensions instead of face grid

Start in src/Geometry/MappedDiscreteModels.jl at Grid(::Type{ReferenceFE{d}}, model::MappedDiscreteModel), then inspect the existing DiscreteModel fallback for lower-dimensional requests. Run the BoundaryTriangulation(model) minimal example from the issue; done means the mapped model returns a face grid for lower dimensions while retaining the cell grid at the model dimension.

Description

In src/Geometry/MappedDiscreteModels.jl, the method

Grid(::Type{ReferenceFE{d}}, model::MappedDiscreteModel)

is implemented as:

function Grid(::Type{ReferenceFE{d}}, model::MappedDiscreteModel) where {d}
  get_grid(model)
end

This implementation matches all values of d and always returns the full cell grid of the model.

However, when d < Dc (for example when constructing boundary or skeleton triangulations), the function should return the corresponding lower-dimensional grid instead of the cell grid.

MappedDiscreteModel is commonly used for mapped geometries, and operations such as

BoundaryTriangulation(model)
SkeletonTriangulation(model)

internally request lower-dimensional grids through calls like:

Grid(ReferenceFE{D-1}, model)

With the current implementation, the method still returns the cell grid, which can lead to incorrect behavior or runtime errors when building boundary or skeleton triangulations.

---

Expected behavior

Grid(::Type{ReferenceFE{d}}, model::MappedDiscreteModel) should only return the cell grid when d equals the topological dimension of the model.

For lower dimensions, the existing fallback defined for DiscreteModel should handle the request.

This can be achieved by constraining the method signature:

function Grid(::Type{ReferenceFE{d}}, model::MappedDiscreteModel{d}) where {d}
  get_grid(model)
end

so that the method only applies when d matches the model dimension.

---

Minimal example

using Gridap
using Gridap.Geometry

model0 = CartesianDiscreteModel((0,1,0,1), (4,4))
phi(x) = VectorValue(x[1] + 0.1*x[1]*x[2], x[2])
model = MappedDiscreteModel(model0, phi)

Γ = BoundaryTriangulation(model)

This call should construct the boundary triangulation correctly using the face grid.