{
 "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](https://gridap.github.io/Tutorials/dev/pages/t007_darcy/)\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.6.7"
  },
  "kernelspec": {
   "name": "julia-1.6",
   "display_name": "Julia 1.6.7",
   "language": "julia"
  }
 },
 "nbformat": 4
}