Interpolation for Raviart Thomas Basis Functions

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:

• Generating Pretty Unicode Plots
• Mixed DG methods for Biharmonic Equation

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.

Constructing RT CellFields

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

\[\begin{align} dof_k(\mathbf{v}_h) = \frac{1}{meas(e_k)} \int_{e_k} \mathbf{v}_h\cdot n_k \,ds\\ \end{align}\]

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.