Postdoc at Umeå University, Sweden

It’s time for a GSoC project update! May 18 to June 8 in the Google Summer of Code timeline was the community bonding period, where the student spends about a month knowing more about the organization and the community. In the previous post, I spoke about the general application process, converging on a topic for the summer of code. In this short post, I will be talking about what goes on during the bonding period. Note that this can vary between the various organizations in GSoC, for example, differing in the tech that they use for communication etc.

I will also try to briefly discuss the technical details of the work that I’ve been working on so far. I don’t intend to explain the entire thing in this single blog post, as it would make it long and hard to follow. I will however try to make a separate page documenting the routines, if possible.

To get started, I needed to know how to contribute to a project. This step is almost uniform throughout the open-source community. There are several articles online that describe the various terminologies I use throughout this post. An excellent place to start with is the official Github’s documentation page

Gridap uses Github for hosting their
projects. The first step was to create a fork of the official
repository. A *fork* is a copy of the remote repository in your account,
which can then be *cloned* to make a local copy on your computer. Any
contributions/changes that you want to make to the code will be done
locally and pushed to the fork. The next step is to submit a pull
request to the original repository from your fork once you commit the
local changes. The admins then review the pull request and merge the
changes if everything looks good.

Aside from communicating through Slack, my mentors advised me to create a task list in the issue tracker of the fork, similar to the one here. This list is to keep track of where I am with the implementation. Github flavored markdown offers support to create a task list.

The first ingredient towards building the fast finite element
interpolator in Gridap is to evaluate the function at arbitrary points
in the domain. The algorithm works by finding the nearest
neighbour
of the point in a KD-Tree, built using the mesh
vertices. Subsequently, we construct an inverse map to get the
corresponding cell/element location and then evaluate the
function using Gridap’s `evaluate`

method. To make things simple, it
is easy to think of points as vertices, and the vertices make up the
cell/element. The goal is to make `evaluate`

work on any point inside
the cell.

This algorithm was written by Erik
Schnetter for Gridap a couple of months
ago. I will go through the algorithm briefly here. Feel free to skip
to the conclusions if you are not familiar with FEM. The first step is
to build a mesh of an arbitrary dimension `D`

. The code snippet below
works for any value of `D`

, subject to computational power of your
computer of course.

```
D = 3
domain = repeat([0, 1], D) # domain = [0,1]^D
partition = repeat([3], D) # 3 cells
model = CartesianDiscreteModel(domain, partition)
```

The next step is to create a finite element space on the mesh.

```
reffe = ReferenceFE(lagrangian, Float64, 2)
V = FESpace(model, reffe)
```

This creates a Lagrangian **F**inite **E**lement **S**pace of
order 2. The reference finite element space is where most of the
computation occurs and this makes the implementation efficient. Now
let us define an affine function on the finite element space. We chose
this for the test, since affine functions are represented exactly on a
second order FESpace.

```
coeff0 = rand(Float64)
coeffs = rand(SVector{D,Float64})
f(x) = coeffs ⋅ SVector(Tuple(x)) + coeff0
fh = interpolate_everywhere(f, V)
```

The `interpolate_everywhere`

method in Gridap interpolates the
function on the degrees of freedom of the finite element space. The
new `evaluate`

method then works if we specify an **arbitrary
point/array of points** in
the domain, say something like

```
evaluate(fh, Point(√2-1, √3-1))
evaluate(fh, [Point(√2-1, √3-1), Point(0.121,0.123)])
```

The signature of the new method would be

```
evaluate!(f::CellField, xs::AbstractVector{<:Point}) ## xs can be a vector of points
```

By default, the `evaluate`

method works for inputs of the form

```
evaluate(get_data(fh), Point(√2-1, √3-1))
```

but returns the piecewise finite element function evaluated on all the
cells/elements. A k-dimensional tree is then used to perform an efficient
search of the cell where the point is located in. This is done mainly
using three new methods present in the `Gridap.CellData`

module

```
_point_to_cell_cache(trian::Triangulation)
_point_to_cell!(cache, x::Point)
make_inverse_table(i2j::AbstractVector{<:Integer},nj::Int)
```

The `_point_to_cell_cache`

method takes in the triangulation (of
Gridap type `Triangulation`

) and returns the `KDTree`

of the mesh
vertices along with other caches required to compute the location of
the point in the triangulation. The method `_point_to_cell!`

takes in the
cache returned by `_point_to_cell_cache`

along with the point itself,
and returns the cell closest to the point. The method `make_inverse_table`

then computes the global and local positions of the point in the
cell, with local being the position of the point in the input array
(remember `evaluate`

now works for an array of points). This is useful
if there are multiple points inside the same cell. Once the
appropriate cell-positions `cell_to_points`

are found, the old
`evaluate`

function can then be used,

```
cell_to_xs = lazy_map(Broadcasting(Reindex(xs)),cell_to_points)
evaluate(get_data(fh), cell_to_xs)
```

where `cell_to_xs`

is the table containing the final positions of the
input points `xs`

in the cell array. This returns the function value
on the queried points `xs`

.

One of the key functions in Gridap is the `lazy_map`

, which is
described in detail in the developer tutorials

In short, `lazy_map`

is one of, if not the most important function in
Gridap. It returns a data structure known as `LazyArray`

which are
arrays that are conceptually defined at the indices. This means that
they are not actually computed as a whole, but **rather contain the
information** on how it is built. The actual value of the array is
computed **only** when the indices are called, and hence the term
**Lazy**. This results in more performant code as opposed to
**greedy** evaluations which result in memory intensive and slower
code.

This work was started a couple of months ago in a pull request by Erik
Schnetter, who wrote
most of the algorithm. I then volunteered to finish up Erik’s
implementation, mostly cleaning up the code, extending the method to
`CellPoint`

structures, and completing the unit tests. The pull request
details are below if you would like to have a look at them.

- Pull Request by Erik Schnetter
- My Pull Requests to Erik’s branch, #1 and #2

With this implementation, the geometrical part is complete, and my new task list is here. The objective of the bonding period was to take time and get introduced to the codebase, read the contributing guidelines, submit and review pull requests etc. Once all the tests have passed, the last step is to add the list of new features in the changelog.

- See this: keep a changelog.

Over the past month, I learned a lot about how contributing to open-source organization works, and also had a chance to interact with the developers of Gridap. I also thank my mentors for guiding me towards submitting my first ever successful pull request! The official coding period starts on June 8. Super stoked for GSoC this year!

tags: