CellField structs and FE interpolation about gridap.jl

[fverdugo]

Yes, this is one of the (at least 2 ingredients) I am aware of that are missing in order to implement the FE spaces and FE functions.

We have:

Type                 Iteration provides at each cell
==========   ===================================
CellValue         Number, FieldValue, etc.
CellArray         Array
CellField          Field (Missing)

It would be nice to be able to iterate on any CellField directly (i.e., without evaluating it) in the same way we iterate on CellValues or CellArrays. That is, iterating over a CellField would generate a Field for each cell. This Field has to be evaluable at a given array of points.

Then, a new concrete implementation of CellField will be CellFieldFromAnalyticalFieldWithGeomap. Iterating over an instance of CellFieldFromAnalyticalFieldWithGeomap will produce a Field for each cell, a Field that is evaluated a the reference space of the corresponding cell.

What do you think?

EDIT:
The requested type already exists. It is called CellFieldFromCompose. The only part missing it to make the CellFields iterable.

[fverdugo]

The second ingredient is the construction of a new CellArray struct that given, the local to global dof info, the free dof values, and the Dirichlet dof values, provides the corresponding Vector of dofs values at each cell. This is trivial. I'll do it now. For the case of no Dirichlet Boundary conditions it is already implemented in CellVectorFromLocalToGlobal.

[santiagobadia]

OK, I will create a draft of CellField and its concrete CellFieldFromAnalyticalFieldWithGeomap

Please, take a look at CellVectorByComposition that I have added to Wrappers.jl in order to express the composition of cell vectors. I have used to work with a NFace-based DOF info, and later compose it with the Cell-NFace info to have a Cell-based VEF info that is needed in the assembly process.

[fverdugo]

It is not needed to create a concrete type. There is already one. It is called CellFieldFromCompose (the name can be changed if needed). I will also think on how to make all CellField types iterable. I think, it should not be very difficult.

[santiagobadia]

Not sure I understand

There is already one. It is called CellFieldFromCompose (the name can be changed if needed).

You say above

Then, a new concrete implementation of CellField will be CellFieldFromAnalyticalFieldWithGeomap. Iterating over an instance of CellFieldFromAnalyticalFieldWithGeomap will produce a Field for each cell, a Field that is evaluated a the reference space of the corresponding cell.

It seems CellFieldFromCompose is a composition of a Function and CellField, not sure it is exactly what we need...

[fverdugo]

I was thinking at the beginning to make all CellField types iterable. In this scenario, CellFieldFromCompose would be iterable as well and could be used for your needs.

Anyway, I have realized that this is not the way to go, since making all CellField types iterable would require some work. Perhaps, the way to go (at least for the moment) is to make a new abstract type called IterCellField or IndexCellFielld (since the discretization used to build the FE space will be indexable). Then, implement the concrete struct CellFieldFromAnalyticalFieldWithGeomap which is iterable, or iterable and indexable.

In the future, we can think if it pays off to make all CellField types iterable as I was thinking at the beginning...

[santiagobadia]

I have put the following sketch of the code we need in CellFieldTests.jl

"""
Cell-wise field created from a `Field`
"""
struct CellFieldFromField{D,T,F<:Field{D,T}} <: CellField{D,T}
  field::F
end

function evaluate(self::CellFieldFromField{D,T},points::CellPoints{D}) where D
  CellFieldValuesFromField(self.field,points)
end

function gradient(self::CellFieldFromField)
  gradfield = gradient(self.field)
  CellFieldFromField(gradfield)
end

"""
Result of the lazy evaluation of a `CellFieldFromField` on a cell-wise array
of points `CellPoints{D}`
"""
struct CellFieldValuesFromField{D,T,F,C} # <: IndexCellArray{T,1,C}
  # @santiagobadia : Not sure from where to extend
  field::F
  cellpoints::C
  cv::CachedVector{T,Vector{T}}
end

function CellFieldValuesFromField(
  field::Field{D,T},cellpoints::CellPoints{D}) where {D,T}
  F = typeof(field)
  C = typeof(cellpoints)
  a = Vector{T}(undef,celllength(cellpoints))
  cv = CachedArray(a)
  CellFieldValuesFromField{D,T,F,C}(field,cellpoints,cv)
end

# @santiagobadia : I am probably missing interface methods

@propagate_inbounds function getindex(self::CellFieldValuesFromField,cell::Int)
  points = self.cellpoints[cell]
  setsize!(self.cv,(length(points),))
  evaluate!(self.field,points,cv)
end

inputcellarray(self::CellFieldValuesFromField) = self.cellpoints

computesize(self::CellFieldValuesFromField, asize) = asize[1]

function computevals!(self::CellFieldValuesFromField, a, v)
  evaluate!(self.field,a,v)
end

[fverdugo]

As far as I understand, this is not the type you need. I think we are facing a fundamental limitation of the current design of CellField that requires some work. Tomorrow I’ll think on this and propose something

[santiagobadia]

I would say that the CellFieldFromField will be needed. I am not sure about CellFieldValuesFromField but it is reasonable to have it for the lazy evaluation of the first one. The first one is the counterpart of CellBasisWithSingleInterpolation and the second one the same for values.

In fact, there is something that I don't like much. I believe that Basis and Field should be understood as a particular cases of Function{D,T,N} with T <: FieldValue and write their code only once, e.g., all the structs in the first paragraph would be only one.

On the other hand, I think we need type inheritance at the assembled type, i.e., the CellFunction level, and the component level, e.g., the component of CellFunction{S,M,T,N} must be a function with range S,M and domain T,N. This is the problem I am facing now for fields, the CellFieldmust provide a Field. Again, if we define FunctionArray{S,M,T,N}, the case I am interested in is a particular case, FunctionArray{Point{D},1,T<:FieldValue,1}.

Finally, when you say I don't need the previous object, I believe you are thinking about something I need and not included in the previous sketch. I still need a "FE-function" that when restricted to a cell will provide me a Field. It is what we need for interpolation, after composition with another "Fe-function" for the geomap as composition of cell-wise fields. When defined this object, the interpolation will be easy. Again, the composition method can be defined at CellFunction level.

[santiagobadia]

For instance, to have this CellFunction instead of CellBasisWithSingleInterpolation and CellFieldFromField above (just define constants).

"""
Concrete implementation of CellFunction for the case of the same function on all cells
"""
struct ConstantCellFunction{S,M,T,N}
  A::Function{S,M,T,N}  # Or put an acceptable name here
end

function CellBasisFromSingleInterpolation(f::ConstantCellFunction{S,M,T,N}) where {S,M,T,N}
  ConstantCellFunction{S,M,T,N}(f)
end

[fverdugo]

I would say that the CellFieldFromField will be needed.

Yes, sure. It will be useful in some context and it makes sense to add it. My point was that this is not the precise type needed for the FE interpolation.

In fact, there is something that I don't like much. I believe that Basis and Field should be understood as a particular cases of Function{D,T,N} with T <: FieldValue and write their code only once, e.g., all the structs in the first paragraph would be only one.

I am also aware of this. We need an umbrella type for Field and Basis. Note that we cannot named it Function since it is already defined in julia....

I still need a "FE-function" that when restricted to a cell will provide me a Field. It is what we need for interpolation, after composition with another "Fe-function" for the geomap as composition of cell-wise fields. When defined this object, the interpolation will be easy. Again, the composition method can be defined at CellFunction level.

Yes sure. This is what I was trying to explain with the table:

Type                 Iteration provides at each cell
==========   ===================================
CellValue         Number, FieldValue, etc.
CellArray         Array
CellField          Field (Missing)

Now we can iterate in a CellValue (and provides a e.g., FieldValue at the cell level). We can iterate a CellArray (and provides an Array at the cell level). however, now we cannot iterate a CellFunction/CellField/CellBasis. We can only evaluate it on a CellArray. This is the missing part. We need to iterate CellFunction/CellField/CellBasis objects, which furnish a Function?/Field/Basis at the cell level.

My only concern is whether we need to iterate all CellFunction/CellField/CellBasis objects or only a subset of them... Making all of them iterable is the elegant way to go. Making iterable only a subset of them is a shortcut that can be useful for the moment.

[santiagobadia]

OK, I think we agree in the diagnostic.
I am not sure it has much sense to waste the time looking for shortcuts.
It is probably better to start thinking about the right solution.
I will create a list of things (I think) are needed.

[fverdugo]

OK. Totally agree by doing things right from the beginning. I am now tiding up the geometry-related part.

[santiagobadia]

• Create a new object Map{S,M,T,N} that given an M-dim array of S provides an N-dim array of T
• Define its constants for Basis and Field
• Create an abstract type iterable IterCellMap{S,M,T,N} that returns in iterate and object of type Map.
• Create an abstract type indexable IndexCellMap{S,M,T,N,V<:Map{S,M,T,N}}<:AbstractVector{V} that returns in getindex and object of type V and implements the iterate methods.
• const CellMap{S,M,T,N} = Union{IterCellMap{S,M,T,N},indexCellMap{S,M,T,N}} where {S,M,T,N}
• Replace CellBasisWithSingleInterpolation with a struct ConstantCellMap{S,M,T,N,A<:Map{S,M,T,N}}<:IterCellMap{S,M,T,N,A}, which will be useful to create, e.g., a CellField with an AnalyticalField.
• Create an struct IterCellMapValues{S,M,T,N,A<:CellMap{S,M,T,N},B<:CellArray{S,M}}<:IterCellArray{T,N} where the result of CellMap is stored (lazy evaluation). Perhaps also for IndexCellMapValues{S,M,T,N,A<:IndexCellMap{S,M,T,N},B<:IndexCellArray{...}<:IndexCellArray{T,N,...} The CellBasisValues objects and the not made yet CellFieldValues not needed, particular cases
• Create a CellMapCompose by composition of maps
• Make the CellFieldFromExtend (i.e., a FE function) to fit in the previous setting

For the moment, I would create new objects, with new names, and next, replace CellFunction with CellMap (or the other way around). But I would like to use the same name in the aggregated object and object and now Function already in Julia.

[santiagobadia]

@fverdugo please take a look at the list and edit it if you want

[fverdugo]

@santiagobadia I have edited a bit the list. The extra type parameters will be needed to have type stable structs... Remove the bold face markup once you agreed.

[fverdugo]

Essentially, at the end CellMap follows the same desing as CellValue and CellArray, which I find a good thing.

[santiagobadia]

OK, the only thing I don't get is why the indexable implements the iterate method. I thought that with getindex and the fact that <:AbstractArray the iterate is not needed, it is already provided by Julia. I will read in detain the interfaces for abstract iterators and arrays to understand this point.

[fverdugo]

Yes of course. Its a typo. At the beginning, IndexCellMap was not extending AbstractVector but IterCellMap. So, it was required to implement the iterate interface...Now implementing getindex et al will be sufficient.