Postdoc at Umeå University, Sweden

It’s time for another GSoC project update. During the last couple of weeks, I have been working on building the method for performing interpolations between different meshes. I will be explaining this task list in detail and finish up by describing where this will eventually lead to. I will also briefly talk about some of the challenges I faced as I was going through the list. There are six tasks in the list, and I have listed them down below:

- Build a cartesian mesh along with the FE Space etc.
- Interpolate and exact solution in the FE Space.
- Create a new mesh using the old mesh. Use the
`map`

option to perturb the original mesh. Example meshes here in Fig 6.1. - Evaluate the old function in the new node positions.
- Create an
`FEFunction`

struct using the nodal values in the previous step. - Check the result of the new
`FEFunction`

.

Two things before we start. I am maintaining another repository containing some of the code I use to test my ideas. To explain the tasks in the list, I will be considering this code here:

Also, I have split this topic into two blog posts to better explain the tasks individually. I will be discussing the first three tasks in this post. For the other three tasks, head on here. Here we go!

The first step is to build a cartesian mesh along with the finite element space. This is a simple process which can be done by running

```
# Build the cartesian model
domain = (0,1,0,1)
partition = (10,10)
model = CartesianDiscreteModel(domain, partition)
```

to build the model and running

```
# Build a second order reference finite element space
reffe = ReferenceFE(lagrangian, Float64, 2)
# Next build an FESpace on the mesh using the reference finite element
V₁ = FESpace(model, reffe)
Ω₁ = Triangulation(model)
```

to construct the finite element space. This `FESpace`

contains second
order Lagrangian finite elements. I will be using this space as the
*source*, where I will be interpolating `FEFunction`

defined on this
space to other Lagrangian `FESpaces`

. Simple enough!

The next step is to interpolate an exact solution onto the new
`FESpace`

. This is done by using the `interpolate_everywhere`

method
available in Gridap.

```
# An exact solution
f(x) = x[1] + x[2]
# A finite element (discrete) function defined on V₁.
fₕ = interpolate_everywhere(f, V₁)
```

Now, `fₕ`

is a finite element function defined cell-wise inside
`V₁`

. So now, running for example

```
julia> get_cell_coordinates(Ω₁)
20×20 LazyArray{FillArrays.Fill{Broadcasting{Reindex{Gridap.Geometry.CartesianCoordinates{2, Float64, typeof(identity)}}}, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, Vector{VectorValue{2, Float64}}, 2, Tuple{Gridap.Geometry.CartesianCellNodes{2}}}:
[(0.0, 0.0), (0.05, 0.0), (0.0, 0.05), (0.05, 0.05)] ⋯
⋮
julia> get_cell_dof_values(fₕ)
400-element LazyArray{FillArrays.Fill{Broadcasting{PosNegReindex{Vector{Float64}, Vector{Float64}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{Float64}, 1, Tuple{Table{Int32, Vector{Int32}, Vector{Int32}}}}:
[0.0, 0.05, 0.05, 0.1]
⋮
```

The second command generates a `LazyArray`

containing the cell-wise
degree of freedom values of `fₕ`

, which in case of Lagrangian Elements
is simply the point-value of functions. This can differ in case of
other elements, for example, the Raviart-Thomas elements which will be
considered later.

So far so good … Now we have built the associated

Sourceproblem. The next step is to define the target mesh and finite element space.

The target mesh in reality could be anything. But for our example
here, Gridap offers a way to create different meshes by providing an
option `map`

to specify a function that maps the Cartesian points to
the new mesh points.

```
partition = (40,40)
model₂ = CartesianDiscreteModel(domain,partition; map=sinusoidal_map)
```

The function `sinusoidal_map`

defines the map that defines the position of the
new mesh points. I considered two different maps,

- A sinusoidal map.
- A random perturbation.

The function `sinusoidal_map`

for the sinusoidal mapping defined below

```
function sinusoidal_map(p::Point)
r, s = p
x = r + 0.08*sin(2π*r)*sin(2π*s)
y = s + 0.08*sin(2π*r)*sin(2π*s)
Point(x,y)
end
```

The function takes the point `p`

from the old mesh, applies the
mapping and returns the new point. The map is taken from Figure 6.1
in this paper
here.

The random map is much cooler and is considered in the script here. Do take a look!

The three meshes are shown in the Figure below. The **leftmost** image is
the source, the **middle** is the sinusoidal mesh and the **rightmost**
image is the randomly perturbed mesh.

Now we are set to perform the interpolation. See my next blog post for Tasks #4 to #6.

tags: