add function get_return_type_and_codim

src/cellfunction/cellfunction.jl

@@ -638,10 +638,26 @@ LazyOperators.materialize(a::AbstractArray{<:Number}, ::AbstractLazy) = a
 LazyOperators.materialize(a::Number, ::AbstractLazy) = a
 LazyOperators.materialize(a::LinearAlgebra.UniformScaling, ::AbstractLazy) = a
 
-# TODO : find a better place and a better granularity
-function get_return_type(f, domain::AbstractDomain)
-    f1 = materialize(f, first(DomainIterator(domain)))
-    elementPoint = get_dummy_element_point(domain)
-    value = materialize(f1, elementPoint)
-    return eltype(value)
+"""
+    get_return_type_and_codim(f::AbstractLazy, elementInfo::AbstractDomainIndex)
+    get_return_type_and_codim(f::AbstractLazy, domain::AbstractDomain)
+    get_return_type_and_codim(f::AbstractLazy, mesh::AbstractMesh)
+
+Evaluate the returned type and the codimension of `materialize(f,x)` where
+`x` is a `ElementPoint` from `elementInfo` (or `domain`/`mesh`).
+The returned codimension is always a Tuple, even for a scalar.
+"""
+function get_return_type_and_codim(f::AbstractLazy, elementInfo::AbstractDomainIndex)
+    fₑ = materialize(f, elementInfo)
+    elementPoint = get_dummy_element_point(elementInfo)
+    value = materialize(fₑ, elementPoint)
+    N = value isa Number ? (1,) : size(value)
+    T = eltype(value)
+    return T, N
+end
+function get_return_type_and_codim(f::AbstractLazy, domain::AbstractDomain)
+    get_return_type_and_codim(f, first(DomainIterator(domain)))
+end
+function get_return_type_and_codim(f::AbstractLazy, mesh::AbstractMesh)
+    get_return_type_and_codim(f, CellInfo(mesh, 1))
 end

src/cellfunction/eval_point.jl

@@ -221,18 +221,25 @@ ElementPoint(x, elementInfo::CellInfo, ds) = CellPoint(x, elementInfo, ds)
 ElementPoint(x, elementInfo::FaceInfo, ds) = FacePoint(x, elementInfo, ds)
 
 """
+    get_dummy_element_point(elementInfo::AbstractDomainIndex)
     get_dummy_element_point(domain::AbstractDomain)
 
 Return a `CellPoint` (or a `FacePoint` depending on the
-type of `domain`) associated with the first element on the
-domain and whose coordinates correspond to the origin of
-its `ReferenceDomain`.
-This utility function can be used to easily materialize
-a `CellFunction` and know the type of the result for example.
+type of `elementInfo`) whose coordinates are equals
+to the center of the reference shape of the element.
+
+For a `domain` argument, the point is built from its firt element.
+
+# Devs notes:
+These utility functions can be used to easily materialize
+a `CellFunction` and get the type of the result for example.
 """
+function get_dummy_element_point(elementInfo::AbstractDomainIndex)
+    x = center(shape(get_element_type(elementInfo)))
+    ElementPoint(x, elementInfo, ReferenceDomain())
+end
+
 function get_dummy_element_point(domain::AbstractDomain)
-    mesh = get_mesh(domain)
     elementInfo = first(DomainIterator(domain))
-    xnode1 = get_coords(get_nodes(mesh, 1))
-    ElementPoint(zero(xnode1), elementInfo, ReferenceDomain())
+    get_dummy_element_point(elementInfo)
 end

src/feoperator/projection.jl

@@ -10,7 +10,7 @@ all the data obtained from the neighbor cells of this node. We could use
 the surface area (among other possible choices).
 """
 function var_on_vertices(f::AbstractLazy, mesh::Mesh)
-    N, T = _codim_and_type(f, mesh)
+    T, N = get_return_type_and_codim(f, mesh)
     @assert length(N) <= 1 "N = $(length(N)) > 1 not supported yet"
     values = zeros(T, nnodes(mesh), N[1])
     _var_on_vertices!(values, f, mesh)
@@ -51,29 +51,6 @@ function _var_on_vertices!(values, f::AbstractLazy, mesh::Mesh)
     end
 end
 
-"""
-    _codim_and_type(f::AbstractLazy, mesh::Mesh)
-
-Evaluate the codimension of `f` and returned type. The returned codimension
-is always a Tuple of codimension(s), even for a scalar.
-"""
-function _codim_and_type(f::AbstractLazy, mesh::Mesh)
-    # Get info about first cell of the mesh
-    cInfo = CellInfo(mesh, 1)
-
-    # Materialize FE function on CellInfo
-    _f = materialize(f, cInfo)
-
-    # Evaluate the function on the center of the cell
-    cPoint = CellPoint(center(shape(celltype(cInfo))), cInfo, ReferenceDomain())
-    value = materialize(_f, cPoint)
-
-    # Codim and type
-    N = value isa Number ? (1,) : size(value)
-    T = eltype(value)
-    return N, T
-end
-
 """
     var_on_centers(f::AbstractLazy, mesh::AbstractMesh)
 
@@ -82,7 +59,7 @@ Interpolate solution on mesh centers.
 The result is a (ncells, ncomps) matrix if ncomps > 1, or a (ncells) vector otherwise.
 """
 function var_on_centers(f::AbstractLazy, mesh::AbstractMesh)
-    N, T = _codim_and_type(f, mesh)
+    T, N = get_return_type_and_codim(f, mesh)
     @assert length(N) <= 1 "N = $(length(N)) > 1 not supported yet"
     values = zeros(T, ncells(mesh), N[1])
     _var_on_centers!(values, f, mesh)

src/feoperator/projection_newapi.jl

@@ -50,7 +50,8 @@ function projection_l2!(u::AbstractSingleFieldFEFunction, f, dΩ::Measure; mass
         A = mass
     end
     l(v) = ∫(f ⋅ v)dΩ
-    b = assemble_linear(l, V; T = get_return_type(f, get_domain(dΩ)))
+    T, = get_return_type_and_codim(f, get_domain(dΩ))
+    b = assemble_linear(l, V; T = T)
     x = A \ b
     set_dof_values!(u, x)
     return nothing

src/mesh/domain.jl

@@ -233,6 +233,13 @@ function BoundaryFaceDomain(mesh::AbstractMesh, args...; kwargs...)
     BoundaryFaceDomain(parent(mesh), args...; kwargs...)
 end
 
+"""
+    abstract type AbstractDomainIndex
+
+# Devs notes
+All subtypes should implement the following functions:
+- get_element_type(::AbstractDomainIndex)
+"""
 abstract type AbstractDomainIndex end
 
 abstract type AbstractCellInfo <: AbstractDomainIndex end
@@ -250,6 +257,7 @@ end
 @inline celltype(c::CellInfo) = c.ctype
 @inline nodes(c::CellInfo) = c.nodes
 @inline get_nodes_index(c::CellInfo) = c.c2n
+get_element_type(c::CellInfo) = celltype(c)
 
 """ Legacy constructor for CellInfo : no information about node indices """
 CellInfo(icell, ctype, nodes) = CellInfo(icell, ctype, nodes, nothing)
@@ -285,6 +293,9 @@ end
 @inline celltype(c::CellSide) = c.ctype
 @inline nodes(c::CellSide) = c.nodes
 @inline cell2nodes(c::CellSide) = c.c2n
+get_element_type(c::CellSide) = celltype(c)
+
+abstract type AbstractFaceInfo <: AbstractDomainIndex end
 
 """
     FaceInfo{CN<:CellInfo,CP<:CellInfo,FT,FN,F2N,I}
@@ -302,7 +313,7 @@ are duplicate from the negative ones.
 is stored explicitely in `FaceInfo` even if it could have been
 computed by collecting info from the side of the negative or positive cells.
 """
-struct FaceInfo{CN <: CellInfo, CP <: CellInfo, FT, FN, F2N, I}
+struct FaceInfo{CN <: CellInfo, CP <: CellInfo, FT, FN, F2N, I} <: AbstractFaceInfo
     cellinfo_n::CN
     cellinfo_p::CP
     cellside_n::Int
@@ -379,6 +390,7 @@ get_cellinfo_p(faceInfo::FaceInfo) = faceInfo.cellinfo_p
 @inline get_nodes_index(faceInfo::FaceInfo) = faceInfo.f2n
 get_cell_side_n(faceInfo::FaceInfo) = faceInfo.cellside_n
 get_cell_side_p(faceInfo::FaceInfo) = faceInfo.cellside_p
+get_element_type(c::FaceInfo) = facetype(c)
 
 """
 Return the opposite side of the `FaceInfo` : cellside "n" because cellside "p"

src/writers/vtk.jl

@@ -386,11 +386,11 @@ function write_vtk_lagrange(
     nd = ndofs(dhl_export)
 
     # Get ncomps and type of each `point` variable
-    dimtype = map(var -> _codim_and_type(var, mesh), values(vars_point))
+    type_dim = map(var -> get_return_type_and_codim(var, mesh), values(vars_point))
 
     # VTK stuff
     coords_vtk = zeros(spacedim(mesh), nd)
-    values_vtk = map(_dimtype -> zeros(last(_dimtype), first(_dimtype)..., nd), dimtype)
+    values_vtk = map(((_t, _d),) -> zeros(_t, _d..., nd), type_dim)
     cells_vtk = MeshCell[]
     sizehint!(cells_vtk, ncells(mesh))
     nodeweigth_vtk = zeros(nd)

test/interpolation/test_projection.jl

@@ -85,13 +85,13 @@
 
     @testset "misc" begin
         mesh = one_cell_mesh(:quad)
-        N, T = Bcube._codim_and_type(PhysicalFunction(x -> x[1]), mesh)
+        T, N = Bcube.get_return_type_and_codim(PhysicalFunction(x -> x[1]), mesh)
         @test N == (1,)
         @test T == Float64
-        N, T = Bcube._codim_and_type(PhysicalFunction(x -> 1), mesh)
+        T, N = Bcube.get_return_type_and_codim(PhysicalFunction(x -> 1), mesh)
         @test N == (1,)
         @test T == Int
-        N, T = Bcube._codim_and_type(PhysicalFunction(x -> x), mesh)
+        T, N = Bcube.get_return_type_and_codim(PhysicalFunction(x -> x), mesh)
         @test N == (2,)
         @test T == Float64
     end