Postdoc at Umeå University, Sweden

GSoC project update time! So after the first evaluation, I continued
working on the project to extend the interpolation algorithm to
Raviart-Thomas elements. In this post I will be summarizing how I went
about doing that. Before I begin, I want to point you to couple of
things. The main GSoC blog
page has a list that
summarizes my project progress. I figured I should have it in the main
page for easy access. You can also find some **Gridap.jl** code
that I wrote during the project. Second are the two other blog posts:

containing details of a couple of side-projects that I have been working on. One is showcasing Julia’s multiple dispatch applied to a problem that I worked on during my PhD. The second one is on convergence of mixed DG methods for biharmonic equation. Do take a look if you’re interested.

In the previous blog
post here, I
talked about how to build the `CellField`

objects for Lagrange
elements. Recalling the code block,

```
function _cell_vals(fₕ, V₂)
phys_point = get_cell_points(get_fe_dof_basis(V₂)).cell_phys_point
fₕ_phys_coords(x) = evaluate(fₕ, x)
phys_point_fx = lazy_map(fₕ_phys_coords, phys_point)
gₕ = CellField(V₂, phys_point_fx)
end
```

This works fine if `V₂`

is a Lagrange `FESpace`

, since the nodal dofs
there are simply point values at the nodes. Raviart Thomas functions
are vector-valued functions described using moment based dofs like

where \(\mathbf{v}_h \in RT_{0,h}\). Now the above code does not work
because the argument `phys_point_fx`

inside `CellField`

**are not
the** degrees of freedom anymore. They are point-wise values of `fₕ`

evaluated at the new mesh points. So an extra step to convert the
point values into dofs is necessary. Luckily we can use one of the
methods of `evaluate!`

that is defined in Gridap.

So a more generalized method would be to

```
# Interpolatable is a new object that will be discussed in detail in upcoming posts.
function FESpaces._cell_vals(V::SingleFieldFESpace, f::Interpolatable)
fe_basis = get_fe_dof_basis(V)
fh = f.uh
trian = get_triangulation(V)
cell_maps = get_cell_map(trian)
bs = get_data(fe_basis)
cache = return_cache(testitem(bs), fh)
if(DomainStyle(fe_basis) == ReferenceDomain())
cell_vals = lazy_map(i -> evaluate!(cache, _to_physical_domain(bs[i], cell_maps[i]), fh), 1:num_cells(trian))
else
cell_vals = lazy_map(i -> evaluate!(cache, bs[i], fh), 1:num_cells(trian))
end
end
```

get the basis functions of the new `FESpace`

and use the method

```
evaluate!(cache, b::MomentBasedDofBasis, field::Field)
```

defined in Gridap. However, inside `_cell_vals`

, the routine `b`

needs
to be in the physical domain. So if `b`

exists in the reference
domain, a function to convert them to the physical domain is
essential. To this end, we define these two methods for
`MomentDofBasis`

and `LagrangeDofBasis`

:

```
function _to_physical_domain(b::MomentBasedDofBasis, cell_map::Field)
nodes = b.nodes
f_moments = b.face_moments
f_phys_nodes = cell_map(nodes)
MomentBasedDofBasis(f_phys_nodes, f_moments, b.face_nodes)
end
function _to_physical_domain(b::LagrangianDofBasis, cell_map::Field)
nodes = b.nodes
f_phys_nodes = cell_map(nodes)
LagrangianDofBasis(f_phys_nodes, b.dof_to_node, b.dof_to_comp, b.node_and_comp_to_dof)
end
```

An interesting tidbit here is the occurrence of multiple dispatch:

```
cache = return_cache(testitem(bs), fh)
```

returns the KDTree and the necessary data structures needed to
evaluate `fh`

on arbitrary points. This `cache`

then works directly
with `evaluate!`

without any extra definition, since
`evaluate(fh::Field, arbitrary_point::Point)`

is defined as a new
feature. The

```
evaluate!(cache, b::MomentBasedDofBasis, field::Field) # Uses evaluate(fh::Field, arbitrary_point::Point)
```

uses the former to evaluate the function `field`

on
`b.nodes`

. The only thing I need to do was to ensure `b.nodes`

is
defined in the right domain. Neat! Also neat is what my mentor Eric
Neiva came up with

```
gₕ = evaluate(x -> fh(x), V₂)
```

which works always, but not optimized since it builds the KDTree for
each point in the triangulation of `V₂`

. We are getting close to a
good solution which can be added to Gridap soon. Check out the PR
Here.