{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Weakly imposing a Dirichlet boundary condition\n", "\n", "This tutorial shows how to implement the weak imposition of a Dirichlet boundary condition, as proposed in the paper Boundary Element Methods with Weakly Imposed Boundary Conditions (2019)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we import Bempp and NumPy." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import bempp.api\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we define the grid for our problem, and the function spaces that we will use. In this example, we use a sphere with P1 and DUAL0 function spaces." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "h = 0.3\n", "grid = bempp.api.shapes.sphere(h=h)\n", "p1 = bempp.api.function_space(grid, \"P\", 1)\n", "dual0 = bempp.api.function_space(grid, \"DUAL\", 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we define the blocked operators proposed in the paper:\n", "$$\\left(\\left(\\begin{array}{cc}-\\mathsf{K}&\\mathsf{V}\\\\\\mathsf{W}&\\mathsf{K}'\\end{array}\\right)+\\left(\\begin{array}{cc}\\tfrac12\\mathsf{Id}&0\\\\\\beta\\mathsf{Id}&-\\tfrac12\\mathsf{Id}\\end{array}\\right)\\right)\\left(\\begin{array}{c}u\\\\\\lambda\\end{array}\\right)=\\left(\\begin{array}{c}g_\\text{D}\\\\\\beta g_\\text{D}\\end{array}\\right),$$\n", "where $\\beta>0$ is a parameter of our choice. In this example, we use $\\beta=0.1$." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "beta = 0.1\n", "multi = bempp.api.BlockedOperator(2,2)\n", "multi[0,0] = -bempp.api.operators.boundary.laplace.double_layer(p1, p1, dual0, assembler=\"fmm\")\n", "multi[0,1] = bempp.api.operators.boundary.laplace.single_layer(dual0, p1, dual0, assembler=\"fmm\")\n", "multi[1,0] = bempp.api.operators.boundary.laplace.hypersingular(p1, dual0, p1, assembler=\"fmm\")\n", "multi[1,1] = bempp.api.operators.boundary.laplace.adjoint_double_layer(dual0, dual0, p1, assembler=\"fmm\")\n", "\n", "diri = bempp.api.BlockedOperator(2,2)\n", "diri[0,0] = 0.5 * bempp.api.operators.boundary.sparse.identity(p1, p1, dual0)\n", "diri[1,0] = beta * bempp.api.operators.boundary.sparse.identity(p1, dual0, p1)\n", "diri[1,1] = -0.5 * bempp.api.operators.boundary.sparse.identity(dual0, dual0, p1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we define the function $g_\\text{D}$, and define the right hand side.\n", "\n", "Here, we use $$g_\\text{D}=\\sin(\\pi x)\\sin(\\pi y)\\sinh(\\sqrt2\\pi z),$$ as in section 5 of the paper." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Numba kernel time: 3.8878002166748047\n", "Numba kernel time: 2.7627720832824707\n", "Numba kernel time: 0.00045013427734375\n" ] } ], "source": [ "@bempp.api.real_callable\n", "def f(x, n, d, res):\n", " res[0] = np.sin(np.pi*x[0]) * np.sin(np.pi*x[1]) * np.sinh(np.sqrt(2)*np.pi*x[2])\n", "\n", "f_fun = bempp.api.GridFunction(p1, fun=f)\n", "\n", "rhs = [2*diri[0,0]*f_fun, diri[1,0]*f_fun]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we solve the system. We set `use_strong_form=True` to activate mass matrix preconditioning." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Numba kernel time: 0.002165555953979492\n", "Numba kernel time: 1.853149175643921\n", "Numba kernel time: 0.0011739730834960938\n", "Numba kernel time: 0.0011909008026123047\n", "Numba kernel time: 0.001819610595703125\n", "Solution took 18 iterations\n" ] } ], "source": [ "sol, info, it_count = bempp.api.linalg.gmres(multi+diri, rhs, return_iteration_count=True, use_strong_form=True)\n", "print(f\"Solution took {it_count} iterations\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For this problem, we know the analytic solution. We compute the error in the $\\mathcal{B}_\\text{D}$ norm:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Numba kernel time: 0.0009734630584716797\n", "Numba kernel time: 0.0009589195251464844\n", "Error: 7.761257842012365\n" ] } ], "source": [ "@bempp.api.real_callable\n", "def g(x, n, d, res):\n", " grad = np.array([\n", " np.cos(np.pi*x[0]) * np.sin(np.pi*x[1]) * np.sinh(np.sqrt(2)*np.pi*x[2]) * np.pi,\n", " np.sin(np.pi*x[0]) * np.cos(np.pi*x[1]) * np.sinh(np.sqrt(2)*np.pi*x[2]) * np.pi,\n", " np.sin(np.pi*x[0]) * np.sin(np.pi*x[1]) * np.cosh(np.sqrt(2)*np.pi*x[2]) * np.pi * np.sqrt(2)\n", " ])\n", " res[0] = np.dot(grad, n)\n", "\n", "g_fun = bempp.api.GridFunction(dual0, fun=g)\n", "\n", "e_fun = [sol[0]-f_fun,sol[1]-g_fun]\n", "\n", "error = 0\n", "# V norm\n", "slp = bempp.api.operators.boundary.laplace.single_layer(dual0, p1, dual0, assembler=\"fmm\")\n", "hyp = bempp.api.operators.boundary.laplace.hypersingular(p1, dual0, p1, assembler=\"fmm\")\n", "error += np.sqrt(np.dot(e_fun[1].coefficients.conjugate(),(slp * e_fun[1]).projections(dual0)))\n", "error += np.sqrt(np.dot(e_fun[0].coefficients.conjugate(),(hyp * e_fun[0]).projections(p1)))\n", "# D part\n", "error += beta**.5 * e_fun[0].l2_norm()\n", "\n", "print(f\"Error: {error}\")" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" } }, "nbformat": 4, "nbformat_minor": 2 }