{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Tutorial 15: Interpolation of CellFields"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "In this tutorial, we will look at how to\n",
    "- Evaluate `CellFields` at arbitrary points\n",
    "- Interpolate finite element functions defined on different\n",
    "triangulations. We will consider examples for\n",
    "   - Lagrangian finite element spaces\n",
    "   - Raviart Thomas finite element spaces\n",
    "   - Vector-Valued Spaces\n",
    "   - Multifield finite element spaces"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Problem Statement\n",
    "Let $\\mathcal{T}_1$ and $\\mathcal{T}_2$ be two triangulations of a\n",
    "domain $\\Omega$. Let $V_i$ be the finite element space defined on\n",
    "the triangulation $\\mathcal{T}_i$ for $i=1,2$. Let $f_h \\in V_1$. The\n",
    "interpolation problem is to find $g_h \\in V_2$ such that\n",
    "\n",
    "$$\n",
    "dof_k^{V_2}(g_h) = dof_k^{V_2}(f_h),\\quad \\forall k \\in\n",
    "\\{1,\\dots,N_{dof}^{V_2}\\}\n",
    "$$"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Setup\n",
    "For the purpose of this tutorial we require `Test`, `Gridap` along with the\n",
    "following submodules of `Gridap`"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Test\n",
    "using Gridap\n",
    "using Gridap.CellData\n",
    "using Gridap.Visualization"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We now create a computational domain on the unit square $[0,1]^2$ consisting\n",
    "of 5 cells per direction"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "domain = (0,1,0,1)\n",
    "partition = (5,5)\n",
    "𝒯₁ = CartesianDiscreteModel(domain, partition)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Background\n",
    "`Gridap` offers the feature to evaluate functions at arbitrary\n",
    "points in the domain. This will be shown in the next\n",
    "section. Interpolation then takes advantage of this feature to\n",
    "obtain the `FEFunction` in the new space from the old one by\n",
    "evaluating the appropriate degrees of freedom. Interpolation works\n",
    "using the composite type `Interpolable` to tell `Gridap` that the\n",
    "argument can be interpolated between triangulations."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Interpolating between Lagrangian FE Spaces"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Let us define the infinite dimensional function"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "f(x) = x[1] + x[2]"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "This function will be interpolated to the source `FESpace`\n",
    "$V_1$. The space can be built using"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "reffe₁ = ReferenceFE(lagrangian, Float64, 1)\n",
    "V₁ = FESpace(𝒯₁, reffe₁)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Finally to build the function $f_h$, we do"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "fₕ = interpolate_everywhere(f,V₁)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "To construct arbitrary points in the domain, we use `Random` package:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Random\n",
    "pt = Point(rand(2))\n",
    "pts = [Point(rand(2)) for i in 1:3]"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The finite element function $f_h$ can be evaluated at arbitrary points (or\n",
    "array of points) by"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "fₕ(pt), fₕ.(pts)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can also check our results using"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "@test fₕ(pt) ≈ f(pt)\n",
    "@test fₕ.(pts) ≈ f.(pts)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now let us define the new triangulation $\\mathcal{T}_2$ of\n",
    "$\\Omega$. We build the new triangulation using a partition of 20 cells per\n",
    "direction. The map can be passed as an argument to\n",
    "`CartesianDiscreteModel` to define the position of the vertices in\n",
    "the new mesh."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "partition = (20,20)\n",
    "𝒯₂ = CartesianDiscreteModel(domain,partition)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "As before, we define the new `FESpace` consisting of second order\n",
    "elements"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "reffe₂ = ReferenceFE(lagrangian, Float64, 2)\n",
    "V₂ = FESpace(𝒯₂, reffe₂)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now we interpolate $f_h$ onto $V_2$ to obtain the new function\n",
    "$g_h$. The first step is to create the `Interpolable` version of\n",
    "$f_h$."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "ifₕ = Interpolable(fₕ)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Then to obtain $g_h$, we dispatch `ifₕ` and the new `FESpace` $V_2$\n",
    "to the `interpolate_everywhere` method of `Gridap`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "gₕ = interpolate_everywhere(ifₕ, V₂)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can also use\n",
    "`interpolate` if interpolating only on the free dofs or\n",
    "`interpolate_dirichlet` if interpolating the Dirichlet dofs of the\n",
    "`FESpace`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "ḡₕ = interpolate(ifₕ, V₂)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The finite element function $\\bar{g}_h$ is the same as $g_h$ in this\n",
    "example since all the dofs are free."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "@test gₕ.cell_dof_values ==  ḡₕ.cell_dof_values"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now we obtain a finite element function using `interpolate_dirichlet`"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "g̃ₕ = interpolate_dirichlet(ifₕ, V₂)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now $\\tilde{g}_h$ will be equal to 0 since there are\n",
    "no Dirichlet nodes defined in the `FESpace`. We can check by running"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "g̃ₕ.cell_dof_values"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Like earlier we can check our results for `gₕ`:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "@test fₕ(pt) ≈ gₕ(pt) ≈ f(pt)\n",
    "@test fₕ.(pts) ≈ gₕ.(pts) ≈ f.(pts)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can visualize the results using Paraview"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "writevtk(get_triangulation(fₕ), \"source\", cellfields=[\"fₕ\"=>fₕ])\n",
    "writevtk(get_triangulation(gₕ), \"target\", cellfields=[\"gₕ\"=>gₕ])"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "which produces the following output"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "![Target](../assets/interpolation_fe/source_and_target.png)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Interpolating between Raviart-Thomas FESpaces"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The procedure is identical to Lagrangian finite element spaces, as\n",
    "discussed in the previous section. The extra thing here is that\n",
    "functions in Raviart-Thomas spaces are vector-valued. The degrees of\n",
    "freedom of the RT spaces are fluxes of the function across the edge\n",
    "of the element. Refer to the\n",
    "tutorial\n",
    "on Darcy equation with RT for more information on the RT\n",
    "elements."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Assuming a function"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "f(x) = VectorValue([x[1], x[2]])"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "on the domain, we build the associated finite dimensional version\n",
    "$f_h \\in V_1$."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "reffe₁ = ReferenceFE(raviart_thomas, Float64, 1) # RT space of order 1\n",
    "V₁ = FESpace(𝒯₁, reffe₁)\n",
    "fₕ = interpolate_everywhere(f, V₁)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "As before, we can evaluate the RT function on any arbitrary point in\n",
    "the domain."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "fₕ(pt), fₕ.(pts)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Constructing the target RT space and building the `Interpolable`\n",
    "object,"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "reffe₂ = ReferenceFE(raviart_thomas, Float64, 1) # RT space of order 1\n",
    "V₂ = FESpace(𝒯₂, reffe₂)\n",
    "ifₕ = Interpolable(fₕ)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "we can construct the new `FEFunction` $g_h \\in V_2$ from $f_h$"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "gₕ = interpolate_everywhere(ifₕ, V₂)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Like earlier we can check our results"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "@test gₕ(pt) ≈ f(pt) ≈ fₕ(pt)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Interpolating vector-valued functions"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can also interpolate vector-valued functions across\n",
    "triangulations. First, we define a vector-valued function on a\n",
    "two-dimensional mesh."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "f(x) = VectorValue([x[1], x[1]+x[2]])"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We then create a vector-valued reference element containing linear\n",
    "elements along with the source finite element space $V_1$."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "reffe₁ = ReferenceFE(lagrangian, VectorValue{2,Float64}, 1)\n",
    "V₁ = FESpace(𝒯₁, reffe₁)\n",
    "fₕ = interpolate_everywhere(f, V₁)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The target finite element space $V_2$ can be defined in a similar manner."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "reffe₂ = ReferenceFE(lagrangian, VectorValue{2,Float64}, 2)\n",
    "V₂ = FESpace(𝒯₂, reffe₂)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The rest of the process is similar to the previous sections, i.e.,\n",
    "define the `Interpolable` version of $f_h$ and use\n",
    "`interpolate_everywhere` to find $g_h \\in V₂$."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "ifₕ = Interpolable(fₕ)\n",
    "gₕ = interpolate_everywhere(ifₕ, V₂)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can then check the results"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "@test gₕ(pt) ≈ f(pt) ≈ fₕ(pt)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Interpolating Multi-field Functions"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Similarly, it is possible to interpolate between multi-field finite element\n",
    "functions. First, we define the components $h_1(x), h_2(x)$ of a\n",
    "multi-field function $h(x)$ as follows."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "h₁(x) = x[1]+x[2]\n",
    "h₂(x) = x[1]"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Next we create a Lagrangian finite element space containing linear\n",
    "elements."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "reffe₁ = ReferenceFE(lagrangian, Float64, 1)\n",
    "V₁ = FESpace(𝒯₁, reffe₁)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Next we create a `MultiFieldFESpace` $V_1 \\times V_1$ and\n",
    "interpolate the function $h(x)$ to the source space $V_1$."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "V₁xV₁ = MultiFieldFESpace([V₁,V₁])\n",
    "fₕ = interpolate_everywhere([h₁, h₂], V₁xV₁)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Similarly, the target multi-field finite element space is created\n",
    "using $\\Omega_2$."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "reffe₂ = ReferenceFE(lagrangian, Float64, 2)\n",
    "V₂ = FESpace(𝒯₂, reffe₂)\n",
    "V₂xV₂ = MultiFieldFESpace([V₂,V₂])"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now, to find $g_h \\in V_2 \\times V_2$, we first extract the components of\n",
    "$f_h$ and obtain the `Interpolable` version of the components."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "fₕ¹, fₕ² = fₕ\n",
    "ifₕ¹ = Interpolable(fₕ¹)\n",
    "ifₕ² = Interpolable(fₕ²)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can then use `interpolate_everywhere` on the `Interpolable`\n",
    "version of the components and obtain $g_h \\in V_2 \\times V_2$ as\n",
    "follows."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "gₕ = interpolate_everywhere([ifₕ¹,ifₕ²], V₂xV₂)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can then check the results of the interpolation, component-wise."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "gₕ¹, gₕ² = gₕ\n",
    "@test fₕ¹(pt) ≈ gₕ¹(pt)\n",
    "@test fₕ²(pt) ≈ gₕ²(pt)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Acknowledgements"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Gridap contributors acknowledge support received from Google,\n",
    "Inc. through the Google Summer of Code 2021 project [A fast finite\n",
    "element interpolator in\n",
    "Gridap.jl](https://summerofcode.withgoogle.com/projects/#6175012823760896)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "---\n",
    "\n",
    "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
   ],
   "metadata": {}
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.8"
  },
  "kernelspec": {
   "name": "julia-1.10",
   "display_name": "Julia 1.10.8",
   "language": "julia"
  }
 },
 "nbformat": 4
}