{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Discontinuous Galerkin heat equation\n",
"\n",
"\n",
"\n",
"*Figure 1*: Temperature field on the unit square with an internal uniform heat source\n",
"solved with inhomogeneous Dirichlet boundary conditions on the left and right boundaries\n",
"and flux on the top and bottom boundaries."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"This example was developed\n",
"as part of the Google summer of code funded project\n",
"[\"Discontinuous Galerkin Infrastructure For the finite element toolbox Ferrite.jl\"](https://summerofcode.withgoogle.com/programs/2023/projects/SLGbRNI5)\n",
"\n",
"## Introduction\n",
"\n",
"This tutorial extends [Tutorial 1: Heat equation](heat_equation.md) by using the\n",
"discontinuous Galerkin method. The reader is expected to have gone through [Tutorial 1:\n",
"Heat equation](heat_equation.md) before proceeding with this tutorial. The main\n",
"differences between the two tutorials are the interface integral terms in the weak form,\n",
"the boundary conditions, and some implementation differences explained in the commented\n",
"program below.\n",
"\n",
"The strong form considered in this tutorial is given as follows\n",
"$$\n",
" -\\boldsymbol{\\nabla} \\cdot [\\boldsymbol{\\nabla}(u)] = 1 \\quad \\textbf{x} \\in \\Omega,\n",
"$$\n",
"\n",
"with the inhomogeneous Dirichlet boundary conditions\n",
"$$\n",
"u(\\textbf{x}) = 1 \\quad \\textbf{x} \\in \\partial \\Omega_D^+ = \\lbrace\\textbf{x} : x_1 = 1.0\\rbrace, \\\\\n",
"u(\\textbf{x}) = -1 \\quad \\textbf{x} \\in \\partial \\Omega_D^- = \\lbrace\\textbf{x} : x_1 = -1.0\\rbrace,\n",
"$$\n",
"and Neumann boundary conditions\n",
"$$\n",
"[\\boldsymbol{\\nabla} (u(\\textbf{x}))] \\cdot \\boldsymbol{n} = 1 \\quad \\textbf{x} \\in \\partial \\Omega_N^+ = \\lbrace\\textbf{x} : x_2 = 1.0\\rbrace, \\\\\n",
"[\\boldsymbol{\\nabla} (u(\\textbf{x}))] \\cdot \\boldsymbol{n} = -1 \\quad \\textbf{x} \\in \\partial \\Omega_N^- = \\lbrace\\textbf{x} : x_2 = -1.0\\rbrace,\n",
"$$\n",
"\n",
"The following definitions of average and jump on interfaces between elements are adopted\n",
"in this tutorial:\n",
"$$\n",
" \\{u\\} = \\frac{1}{2}(u^+ + u^-),\\quad [\\![ u]\\!] = u^+ \\boldsymbol{n}^+ + u^- \\boldsymbol{n}^-\\\\\n",
"$$\n",
"where $u^+$ and $u^-$ are the temperature on the two sides of the interface.\n",
"\n",
"\n",
"> **Derivation of the weak form for homogeneous Dirichlet boundary condition**\n",
">\n",
"> Defining $\\boldsymbol{\\sigma}$ as the gradient of the temperature field the equation can be expressed as\n",
"> $$\n",
"> \\boldsymbol{\\sigma} = \\boldsymbol{\\nabla} (u),\\\\\n",
"> -\\boldsymbol{\\nabla} \\cdot \\boldsymbol{\\sigma} = 1,\n",
"> $$\n",
"> Multiplying by test functions $ \\boldsymbol{\\tau} $ and $ \\delta u $ respectively and integrating\n",
"> over the domain,\n",
"> $$\n",
"> \\int_\\Omega \\boldsymbol{\\sigma} \\cdot \\boldsymbol{\\tau} \\,\\mathrm{d}\\Omega = \\int_\\Omega [\\boldsymbol{\\nabla} (u)] \\cdot \\boldsymbol{\\tau} \\,\\mathrm{d}\\Omega,\\\\\n",
"> -\\int_\\Omega \\boldsymbol{\\nabla} \\cdot \\boldsymbol{\\sigma} \\delta u \\,\\mathrm{d}\\Omega = \\int_\\Omega \\delta u \\,\\mathrm{d}\\Omega,\n",
"> $$\n",
"> Integrating by parts and applying divergence theorem,\n",
"> $$\n",
"> \\int_\\Omega \\boldsymbol{\\sigma} \\cdot \\boldsymbol{\\tau} \\,\\mathrm{d}\\Omega = -\\int_\\Omega u (\\boldsymbol{\\nabla} \\cdot \\boldsymbol{\\tau}) \\,\\mathrm{d}\\Omega + \\int_\\Gamma \\hat{u} \\boldsymbol{\\tau} \\cdot \\boldsymbol{n} \\,\\mathrm{d}\\Gamma,\\\\\n",
"> \\int_\\Omega \\boldsymbol{\\sigma} \\cdot [\\boldsymbol{\\nabla} (\\delta u)] \\,\\mathrm{d}\\Omega = \\int_\\Omega \\delta u \\,\\mathrm{d}\\Omega + \\int_\\Gamma \\delta u \\boldsymbol{\\hat{\\sigma}} \\cdot \\boldsymbol{n} \\,\\mathrm{d}\\Gamma,\n",
"> $$\n",
"> Where $\\boldsymbol{n}$ is the outwards pointing normal, $\\Gamma$ is the union of the elements' boundaries, and $\\hat{u}, \\, \\hat{\\sigma}$ are the numerical fluxes.\n",
"> Substituting the integrals of form\n",
"> $$\n",
"> \\int_\\Gamma q \\boldsymbol{\\phi} \\cdot \\boldsymbol{n} \\,\\mathrm{d}\\Gamma = \\int_\\Gamma [\\![ q]\\!] \\cdot \\{\\boldsymbol{\\phi}\\} \\,\\mathrm{d}\\Gamma + \\int_{\\Gamma^0} \\{q\\} [\\![ \\boldsymbol{\\phi}]\\!] \\,\\mathrm{d}\\Gamma^0,\n",
"> $$\n",
"> where $\\Gamma^0 : \\Gamma \\setminus \\partial \\Omega$, and the jump of the vector-valued field $\\boldsymbol{\\phi}$ is defined as\n",
"> $$\n",
"> [\\![ \\boldsymbol{\\phi}]\\!] = \\boldsymbol{\\phi}^+ \\cdot \\boldsymbol{n}^+ + \\boldsymbol{\\phi}^- \\cdot \\boldsymbol{n}^-\\\\\n",
"> $$\n",
"> with the jumps and averages results in\n",
"> $$\n",
"> \\int_\\Omega \\boldsymbol{\\sigma} \\cdot \\boldsymbol{\\tau} \\,\\mathrm{d}\\Omega = -\\int_\\Omega u (\\boldsymbol{\\nabla} \\cdot \\boldsymbol{\\tau}) \\,\\mathrm{d}\\Omega + \\int_\\Gamma [\\![ \\hat{u}]\\!] \\cdot \\{\\boldsymbol{\\tau}\\} \\,\\mathrm{d}\\Gamma + \\int_{\\Gamma^0} \\{\\hat{u}\\} [\\![ \\boldsymbol{\\tau}]\\!] \\,\\mathrm{d}\\Gamma^0,\\\\\n",
"> \\int_\\Omega \\boldsymbol{\\sigma} \\cdot [\\boldsymbol{\\nabla} (\\delta u)] \\,\\mathrm{d}\\Omega = \\int_\\Omega \\delta u \\,\\mathrm{d}\\Omega + \\int_\\Gamma [\\![ \\delta u]\\!] \\cdot \\{\\hat{\\boldsymbol{\\sigma}}\\} \\,\\mathrm{d}\\Gamma + \\int_{\\Gamma^0} \\{\\delta u\\} [\\![ \\hat{\\boldsymbol{\\sigma}}]\\!] \\,\\mathrm{d}\\Gamma^0,\n",
"> $$\n",
"> Integrating $ \\int_\\Omega [\\boldsymbol{\\nabla} (u)] \\cdot \\boldsymbol{\\tau} \\,\\mathrm{d}\\Omega $ by parts and applying divergence theorem\n",
"> without using numerical flux, then substitute in the equation to obtain a weak form.\n",
"> $$\n",
"> \\int_\\Omega \\boldsymbol{\\sigma} \\cdot \\boldsymbol{\\tau} \\,\\mathrm{d}\\Omega = \\int_\\Omega [\\boldsymbol{\\nabla} (u)] \\cdot \\boldsymbol{\\tau} \\,\\mathrm{d}\\Omega + \\int_\\Gamma [\\![ \\hat{u} - u]\\!] \\cdot \\{\\boldsymbol{\\tau}\\} \\,\\mathrm{d}\\Gamma + \\int_{\\Gamma^0} \\{\\hat{u} - u\\} [\\![ \\boldsymbol{\\tau}]\\!] \\,\\mathrm{d}\\Gamma^0,\\\\\n",
"> \\int_\\Omega \\boldsymbol{\\sigma} \\cdot [\\boldsymbol{\\nabla} (\\delta u)] \\,\\mathrm{d}\\Omega = \\int_\\Omega \\delta u \\,\\mathrm{d}\\Omega + \\int_\\Gamma [\\![ \\delta u]\\!] \\cdot \\{\\hat{\\boldsymbol{\\sigma}}\\} \\,\\mathrm{d}\\Gamma + \\int_{\\Gamma^0} \\{\\delta u\\} [\\![ \\hat{\\boldsymbol{\\sigma}}]\\!] \\,\\mathrm{d}\\Gamma^0,\n",
"> $$\n",
"> Substituting\n",
"> $$\n",
"> \\boldsymbol{\\tau} = \\boldsymbol{\\nabla} (\\delta u),\\\\\n",
"> $$\n",
"> results in\n",
"> $$\n",
"> \\int_\\Omega \\boldsymbol{\\sigma} \\cdot [\\boldsymbol{\\nabla} (\\delta u)] \\,\\mathrm{d}\\Omega = \\int_\\Omega [\\boldsymbol{\\nabla} (u)] \\cdot [\\boldsymbol{\\nabla} (\\delta u)] \\,\\mathrm{d}\\Omega + \\int_\\Gamma [\\![ \\hat{u} - u]\\!] \\cdot \\{\\boldsymbol{\\nabla} (\\delta u)\\} \\,\\mathrm{d}\\Gamma + \\int_{\\Gamma^0} \\{\\hat{u} - u\\} [\\![ \\boldsymbol{\\nabla} (\\delta u)]\\!] \\,\\mathrm{d}\\Gamma^0,\\\\\n",
"> \\int_\\Omega \\boldsymbol{\\sigma} \\cdot [\\boldsymbol{\\nabla} (\\delta u)] \\,\\mathrm{d}\\Omega = \\int_\\Omega \\delta u \\,\\mathrm{d}\\Omega + \\int_\\Gamma [\\![ \\delta u]\\!] \\cdot \\{\\hat{\\boldsymbol{\\sigma}}\\} \\,\\mathrm{d}\\Gamma + \\int_{\\Gamma^0} \\{\\delta u\\} [\\![ \\hat{\\boldsymbol{\\sigma}}]\\!] \\,\\mathrm{d}\\Gamma^0,\n",
"> $$\n",
"> Combining the two equations,\n",
"> $$\n",
"> \\int_\\Omega [\\boldsymbol{\\nabla} (u)] \\cdot [\\boldsymbol{\\nabla} (\\delta u)] \\,\\mathrm{d}\\Omega + \\int_\\Gamma [\\![ \\hat{u} - u]\\!] \\cdot \\{\\boldsymbol{\\nabla} (\\delta u)\\} \\,\\mathrm{d}\\Gamma + \\int_{\\Gamma^0} \\{\\hat{u} - u\\} [\\![ \\boldsymbol{\\nabla} (\\delta u)]\\!] \\,\\mathrm{d}\\Gamma^0 - \\int_\\Gamma [\\![ \\delta u]\\!] \\cdot \\{\\hat{\\boldsymbol{\\sigma}}\\} \\,\\mathrm{d}\\Gamma - \\int_{\\Gamma^0} \\{\\delta u\\} [\\![ \\hat{\\boldsymbol{\\sigma}}]\\!] \\,\\mathrm{d}\\Gamma^0 = \\int_\\Omega \\delta u \\,\\mathrm{d}\\Omega,\\\\\n",
"> $$\n",
"> The numerical fluxes chosen for the interior penalty method are $\\boldsymbol{\\hat{\\sigma}} = \\{\\boldsymbol{\\nabla} (u)\\} - \\alpha([\\![ u]\\!])$ on $\\Gamma$, $\\hat{u} = \\{u\\}$ on the interfaces between elements $\\Gamma^0 : \\Gamma \\setminus \\partial \\Omega$,\n",
"> and $\\hat{u} = 0$ on $\\partial \\Omega$. Such choice results in $\\{\\hat{\\boldsymbol{\\sigma}}\\} = \\{\\boldsymbol{\\nabla} (u)\\} - \\alpha([\\![ u]\\!])$, $[\\![ \\hat{u}]\\!] = 0$, $\\{\\hat{u}\\} = \\{u\\}$, $[\\![ \\hat{\\boldsymbol{\\sigma}}]\\!] = 0$ and the equation becomes\n",
"> $$\n",
"> \\int_\\Omega [\\boldsymbol{\\nabla} (u)] \\cdot [\\boldsymbol{\\nabla} (\\delta u)] \\,\\mathrm{d}\\Omega - \\int_\\Gamma [\\![ u]\\!] \\cdot \\{\\boldsymbol{\\nabla} (\\delta u)\\} \\,\\mathrm{d}\\Gamma - \\int_\\Gamma [\\![ \\delta u]\\!] \\cdot \\{\\boldsymbol{\\nabla} (u)\\} - [\\![ \\delta u]\\!] \\cdot \\alpha([\\![ u]\\!]) \\,\\mathrm{d}\\Gamma = \\int_\\Omega \\delta u \\,\\mathrm{d}\\Omega,\\\\\n",
"> $$\n",
"> Where\n",
"> $$\n",
"> \\alpha([\\![ u]\\!]) = \\mu [\\![ u]\\!]\n",
"> $$\n",
"> Where $\\mu = \\eta h_e^{-1}$, the weak form becomes\n",
"> $$\n",
"> \\int_\\Omega [\\boldsymbol{\\nabla} (u)] \\cdot [\\boldsymbol{\\nabla}] (\\delta u) \\,\\mathrm{d}\\Omega - \\int_\\Gamma [\\![ u ]\\!] \\cdot \\{\\boldsymbol{\\nabla} (\\delta u)\\} + [\\![ \\delta u ]\\!] \\cdot \\{\\boldsymbol{\\nabla} (u)\\} \\,\\mathrm{d}\\Gamma + \\int_\\Gamma \\frac{\\eta}{h_e} [\\![ u]\\!] \\cdot [\\![ \\delta u]\\!] \\,\\mathrm{d}\\Gamma = \\int_\\Omega \\delta u \\,\\mathrm{d}\\Omega,\\\\\n",
"> $$\n",
"Since $\\partial \\Omega$ is constrained with both Dirichlet and Neumann boundary conditions the term $\\int_{\\partial \\Omega} [\\boldsymbol{\\nabla} (u)] \\cdot \\boldsymbol{n} \\delta u \\,\\mathrm{d} \\Omega$ can be expressed as an integral over $\\partial \\Omega_N$, where $\\partial \\Omega_N$ is the boundaries with only prescribed Neumann boundary condition,\n",
"The resulting weak form is given as follows: Find $u \\in \\mathbb{U}$ such that\n",
"$$\n",
" \\int_\\Omega [\\boldsymbol{\\nabla} (u)] \\cdot [\\boldsymbol{\\nabla} (\\delta u)] \\,\\mathrm{d}\\Omega - \\int_{\\Gamma^0} [\\![ u]\\!] \\cdot \\{\\boldsymbol{\\nabla} (\\delta u)\\} + [\\![ \\delta u]\\!] \\cdot \\{\\boldsymbol{\\nabla} (u)\\} \\,\\mathrm{d}\\Gamma^0 + \\int_{\\Gamma^0} \\frac{\\eta}{h_e} [\\![ u]\\!] \\cdot [\\![ \\delta u]\\!] \\,\\mathrm{d}\\Gamma^0 = \\int_\\Omega \\delta u \\,\\mathrm{d}\\Omega + \\int_{\\partial \\Omega_N} ([\\boldsymbol{\\nabla} (u)] \\cdot \\boldsymbol{n}) \\delta u \\,\\mathrm{d} \\partial \\Omega_N,\\\\\n",
"$$\n",
"where $h_e$ is the characteristic size (the diameter of the interface), and $\\eta$ is a large enough positive number independent of $h_e$ [Mu:2014:IP](@cite),\n",
"$\\delta u \\in \\mathbb{T}$ is a test function, and where $\\mathbb{U}$ and $\\mathbb{T}$ are suitable\n",
"trial and test function sets, respectively.\n",
"We use the value $\\eta = (1 + O)^{D}$, where $O$ is the polynomial order and $D$ the\n",
"dimension, in this tutorial.\n",
"\n",
"More details on DG formulations for elliptic problems can be found in [Cockburn:2002:unifiedanalysis](@cite)."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"## Commented Program\n",
"\n",
"Now we solve the problem in Ferrite. What follows is a program spliced with comments.\n",
"\n",
"First we load Ferrite and other packages, and generate grid just like the [heat equation tutorial](heat_equation.md)"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using Ferrite, SparseArrays\n",
"dim = 2;\n",
"grid = generate_grid(Quadrilateral, ntuple(_ -> 20, dim));"
],
"metadata": {},
"execution_count": 1
},
{
"cell_type": "markdown",
"source": [
"We construct the topology information which is used later for generating the sparsity pattern for stiffness matrix."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"topology = ExclusiveTopology(grid);"
],
"metadata": {},
"execution_count": 2
},
{
"cell_type": "markdown",
"source": [
"### Trial and test functions\n",
"`CellValues`, `FacetValues`, and `InterfaceValues` facilitate the process of evaluating values and gradients of\n",
"test and trial functions (among other things). To define\n",
"these we need to specify an interpolation space for the shape functions.\n",
"We use `DiscontinuousLagrange` functions\n",
"based on the two-dimensional reference quadrilateral. We also define a quadrature rule based on\n",
"the same reference element. We combine the interpolation and the quadrature rule\n",
"to `CellValues` and `InterfaceValues` object. Note that `InterfaceValues` object contains two `FacetValues` objects which can be used individually."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"order = 1;\n",
"ip = DiscontinuousLagrange{RefQuadrilateral, order}();\n",
"qr = QuadratureRule{RefQuadrilateral}(2);"
],
"metadata": {},
"execution_count": 3
},
{
"cell_type": "markdown",
"source": [
"For `FacetValues` and `InterfaceValues` we use `FacetQuadratureRule`"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"facet_qr = FacetQuadratureRule{RefQuadrilateral}(2);\n",
"cellvalues = CellValues(qr, ip);\n",
"facetvalues = FacetValues(facet_qr, ip);\n",
"interfacevalues = InterfaceValues(facet_qr, ip);"
],
"metadata": {},
"execution_count": 4
},
{
"cell_type": "markdown",
"source": [
"### Penalty term parameters\n",
"We define functions to calculate the diameter of a set of points, used to calculate the characteristic size $h_e$ in the assembly routine."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"getdistance(p1::Vec{N, T}, p2::Vec{N, T}) where {N, T} = norm(p1 - p2);\n",
"getdiameter(cell_coords::Vector{Vec{N, T}}) where {N, T} = maximum(getdistance.(cell_coords, reshape(cell_coords, (1, :))));"
],
"metadata": {},
"execution_count": 5
},
{
"cell_type": "markdown",
"source": [
"### Degrees of freedom\n",
"Degrees of freedom distribution is handled using `DofHandler` as usual"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"dh = DofHandler(grid)\n",
"add!(dh, :u, ip)\n",
"close!(dh);"
],
"metadata": {},
"execution_count": 6
},
{
"cell_type": "markdown",
"source": [
"However, when generating the sparsity pattern we need to pass the topology and the cross-element coupling matrix when we're using\n",
"discontinuous interpolations. The cross-element coupling matrix is of size [1,1] in this case as\n",
"we have only one field and one DofHandler."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"K = allocate_matrix(dh, topology = topology, interface_coupling = trues(1, 1));"
],
"metadata": {},
"execution_count": 7
},
{
"cell_type": "markdown",
"source": [
"### Boundary conditions\n",
"The Dirichlet boundary conditions are treated\n",
"as usual by a `ConstraintHandler`."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"ch = ConstraintHandler(dh)\n",
"add!(ch, Dirichlet(:u, getfacetset(grid, \"right\"), (x, t) -> 1.0))\n",
"add!(ch, Dirichlet(:u, getfacetset(grid, \"left\"), (x, t) -> -1.0))\n",
"close!(ch);"
],
"metadata": {},
"execution_count": 8
},
{
"cell_type": "markdown",
"source": [
"Furthermore, we define $\\partial \\Omega_N$ as the `union` of the facet sets with Neumann boundary conditions for later use"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"∂Ωₙ = union(\n",
" getfacetset(grid, \"top\"),\n",
" getfacetset(grid, \"bottom\"),\n",
");"
],
"metadata": {},
"execution_count": 9
},
{
"cell_type": "markdown",
"source": [
"### Assembling the linear system\n",
"\n",
"Now we have all the pieces needed to assemble the linear system, $K u = f$. Assembling of\n",
"the global system is done by looping over i) all the elements in order to compute the\n",
"element contributions $K_e$ and $f_e$, ii) all the interfaces to compute their\n",
"contributions $K_i$, and iii) all the Neumann boundary facets to compute their\n",
"contributions $f_e$. All these local contributions are then assembled into the\n",
"appropriate place in the global $K$ and $f$.\n",
"\n",
"#### Local assembly\n",
"We define the functions\n",
"* `assemble_element!` to compute the contributions $K_e$ and $f_e$ of volume integrals\n",
" over an element using `cellvalues`.\n",
"* `assemble_interface!` to compute the contribution $K_i$ of surface integrals over an\n",
" interface using `interfacevalues`.\n",
"* `assemble_boundary!` to compute the contribution $f_e$ of surface integrals over a\n",
" boundary facet using `FacetValues`."
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "assemble_boundary! (generic function with 1 method)"
},
"metadata": {},
"execution_count": 10
}
],
"cell_type": "code",
"source": [
"function assemble_element!(Ke::Matrix, fe::Vector, cellvalues::CellValues)\n",
" n_basefuncs = getnbasefunctions(cellvalues)\n",
" # Reset to 0\n",
" fill!(Ke, 0)\n",
" fill!(fe, 0)\n",
" # Loop over quadrature points\n",
" for q_point in 1:getnquadpoints(cellvalues)\n",
" # Quadrature weight\n",
" dΩ = getdetJdV(cellvalues, q_point)\n",
" # Loop over test shape functions\n",
" for i in 1:n_basefuncs\n",
" δu = shape_value(cellvalues, q_point, i)\n",
" ∇δu = shape_gradient(cellvalues, q_point, i)\n",
" # Add contribution to fe\n",
" fe[i] += δu * dΩ\n",
" # Loop over trial shape functions\n",
" for j in 1:n_basefuncs\n",
" ∇u = shape_gradient(cellvalues, q_point, j)\n",
" # Add contribution to Ke\n",
" Ke[i, j] += (∇δu ⋅ ∇u) * dΩ\n",
" end\n",
" end\n",
" end\n",
" return Ke, fe\n",
"end\n",
"\n",
"function assemble_interface!(Ki::Matrix, iv::InterfaceValues, μ::Float64)\n",
" # Reset to 0\n",
" fill!(Ki, 0)\n",
" # Loop over quadrature points\n",
" for q_point in 1:getnquadpoints(iv)\n",
" # Get the normal to facet A\n",
" normal = getnormal(iv, q_point)\n",
" # Get the quadrature weight\n",
" dΓ = getdetJdV(iv, q_point)\n",
" # Loop over test shape functions\n",
" for i in 1:getnbasefunctions(iv)\n",
" # Multiply the jump by the negative normal to get the definition from the theory section.\n",
" δu_jump = shape_value_jump(iv, q_point, i) * (-normal)\n",
" ∇δu_avg = shape_gradient_average(iv, q_point, i)\n",
" # Loop over trial shape functions\n",
" for j in 1:getnbasefunctions(iv)\n",
" # Multiply the jump by the negative normal to get the definition from the theory section.\n",
" u_jump = shape_value_jump(iv, q_point, j) * (-normal)\n",
" ∇u_avg = shape_gradient_average(iv, q_point, j)\n",
" # Add contribution to Ki\n",
" Ki[i, j] += -(δu_jump ⋅ ∇u_avg + ∇δu_avg ⋅ u_jump) * dΓ + μ * (δu_jump ⋅ u_jump) * dΓ\n",
" end\n",
" end\n",
" end\n",
" return Ki\n",
"end\n",
"\n",
"function assemble_boundary!(fe::Vector, fv::FacetValues)\n",
" # Reset to 0\n",
" fill!(fe, 0)\n",
" # Loop over quadrature points\n",
" for q_point in 1:getnquadpoints(fv)\n",
" # Get the normal to facet A\n",
" normal = getnormal(fv, q_point)\n",
" # Get the quadrature weight\n",
" ∂Ω = getdetJdV(fv, q_point)\n",
" # Loop over test shape functions\n",
" for i in 1:getnbasefunctions(fv)\n",
" δu = shape_value(fv, q_point, i)\n",
" boundary_flux = normal[2]\n",
" fe[i] = boundary_flux * δu * ∂Ω\n",
" end\n",
" end\n",
" return fe\n",
"end"
],
"metadata": {},
"execution_count": 10
},
{
"cell_type": "markdown",
"source": [
"#### Global assembly\n",
"\n",
"We define the function `assemble_global` to loop over all elements and internal facets\n",
"(interfaces), as well as the external facets involved in Neumann boundary conditions."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"function assemble_global(cellvalues::CellValues, facetvalues::FacetValues, interfacevalues::InterfaceValues, K::SparseMatrixCSC, dh::DofHandler, order::Int, dim::Int)\n",
" # Allocate the element stiffness matrix and element force vector\n",
" n_basefuncs = getnbasefunctions(cellvalues)\n",
" Ke = zeros(n_basefuncs, n_basefuncs)\n",
" fe = zeros(n_basefuncs)\n",
" Ki = zeros(n_basefuncs * 2, n_basefuncs * 2)\n",
" # Allocate global force vector f\n",
" f = zeros(ndofs(dh))\n",
" # Create an assembler\n",
" assembler = start_assemble(K, f)\n",
" # Loop over all cells\n",
" for cell in CellIterator(dh)\n",
" # Reinitialize cellvalues for this cell\n",
" reinit!(cellvalues, cell)\n",
" # Compute volume integral contribution\n",
" assemble_element!(Ke, fe, cellvalues)\n",
" # Assemble Ke and fe into K and f\n",
" assemble!(assembler, celldofs(cell), Ke, fe)\n",
" end\n",
" # Loop over all interfaces\n",
" for ic in InterfaceIterator(dh)\n",
" # Reinitialize interfacevalues for this interface\n",
" reinit!(interfacevalues, ic)\n",
" # Calculate the characteristic size hₑ as the face diameter\n",
" interfacecoords = ∩(getcoordinates(ic)...)\n",
" hₑ = getdiameter(interfacecoords)\n",
" # Calculate μ\n",
" μ = (1 + order)^dim / hₑ\n",
" # Compute interface surface integrals contribution\n",
" assemble_interface!(Ki, interfacevalues, μ)\n",
" # Assemble Ki into K\n",
" assemble!(assembler, interfacedofs(ic), Ki)\n",
" end\n",
" # Loop over domain boundaries with Neumann boundary conditions\n",
" for fc in FacetIterator(dh, ∂Ωₙ)\n",
" # Reinitialize facetvalues for this boundary facet\n",
" reinit!(facetvalues, fc)\n",
" # Compute boundary facet surface integrals contribution\n",
" assemble_boundary!(fe, facetvalues)\n",
" # Assemble fe into f\n",
" assemble!(f, celldofs(fc), fe)\n",
" end\n",
" return K, f\n",
"end\n",
"K, f = assemble_global(cellvalues, facetvalues, interfacevalues, K, dh, order, dim);"
],
"metadata": {},
"execution_count": 11
},
{
"cell_type": "markdown",
"source": [
"### Solution of the system\n",
"\n",
"The solution of the system is independent of the discontinuous discretization and the\n",
"application of constraints, linear solve, and exporting is done as usual."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"apply!(K, f, ch)\n",
"u = K \\ f;\n",
"VTKGridFile(\"dg_heat_equation\", dh) do vtk\n",
" write_solution(vtk, dh, u)\n",
"end;"
],
"metadata": {},
"execution_count": 12
},
{
"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.11.3"
},
"kernelspec": {
"name": "julia-1.11",
"display_name": "Julia 1.11.3",
"language": "julia"
}
},
"nbformat": 4
}