General Remapping Operator Application

[jb-mackay]

Remapping operators are applied to fields as sparse matrix-vector multiplication.

From @simonbyrne: For application to general fields (vector or tuple-valued fields) we will need to write our own sparse matrix multiply (e.g. here)

field1 # on space1 (source space)

# allocate field on space2 (target space)
field2 = similar(Fields.coordinate_field(space2), eltype(field1))

# 1. if we just have a scalar field:
mul!(vec(parent(field2)), M, vec(parent(field1)))

# 2. for more general case (vectors or tuple-valued fields) we have to write our own sparse matrix multiply
#  x =  M * y
nzv = nonzeros(M) # vector of non-zero values
rv = rowvals(M) # vector of row indices
for col = 1:ncol
   a = y[col] # (y is a collection of vectors)
   #  for FV you can just do y[1,1,1,1,e]
    # <= need to convert a to a global vector bases (can reuse machinery from DSS)
    # requires getting the local_geometry object in space1
   for j = 1:nzrange(M, col) # the number of non-zero rows in the column of M
      aa = nzv[j]
      # <= convert to the vector basis of space2 (requires local_geometry object)
      x[rv[j]] = x[rv[j]] ⊞ nzv[j] ⊠ aa # <= need to use recursive operation for tuples, structs, vector components
   end
end

[jb-mackay]

For Cartesian spaces, this is a simple broadcast. On the sphere we will need the above that is aware of vector bases.