Dokken/update mesh interface

.github/workflows/build-publish.yml

@@ -20,7 +20,7 @@ jobs:
   build:
     # The type of runner that the job will run on
     runs-on: ubuntu-latest
-    container: dokken92/dolfinx_custom:28112021
+    container: dokken92/dolfinx_custom:18122021
 
     env:
       HDF5_MPI: "ON"

.github/workflows/docker-image.yml

@@ -37,4 +37,4 @@ jobs:
         run: echo ${{ secrets.DOCKERHUB_TOKEN }} | docker login -u ${{ secrets.DOCKERHUB_USERNAME }} --password-stdin
       - name: Push to the DockerHub registry
         run: |
-          docker push dokken92/dolfinx_custom:28112021
+          docker push dokken92/dolfinx_custom:18122021

Dockerfile

@@ -1,4 +1,4 @@
-FROM dokken92/dolfinx_custom:28112021
+FROM dokken92/dolfinx_custom:18122021
 
 # create user with a home directory
 ARG NB_USER

chapter1/fundamentals_code.ipynb

@@ -49,20 +49,7 @@
     "we usually know about the numerical errors _asymptotic properties_. For instance that it is proportional to $h^2$ if $h$ is the size of a cell in the mesh. We can then compare the error on meshes with different $h$-values to see if the asymptotic behavior is correct. This technique will be explained in detail in the chapter [Improving your fenics code](./../chapter4/convergence).\n",
     "\n",
     "However, in cases where we have a solution we know that should have no approximation error, we know that the solution should\n",
-    "be produced to machine precision by the program.\n",
-    "\n",
-    "We start by importing external modules used alongside `dolfinx`."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import numpy\n",
-    "import ufl\n",
-    "from mpi4py import MPI"
+    "be produced to machine precision by the program."
    ]
   },
   {
@@ -71,44 +58,43 @@
    "source": [
     "A major difference between a traditional FEniCS code and a FEniCSx code, is that one is not advised to use the wildcard import. We will see this throughout this first example.\n",
     "## Generating  simple meshes\n",
-    "The next step is to define the discrete domain, _the mesh_. We do this by importing the built-in `dolfinx.generation.UnitSquareMesh` function, that can build a $[0,1]\\times[0,1]$ mesh, consisting of either triangles or quadrilaterals."
+    "The next step is to define the discrete domain, _the mesh_. We do this by importing one of the built-in mesh generators. We will build a unit square mesh, i.e. a mesh spanning $[0,1]\\times[0,1]$. It can consist of either triangles or quadrilaterals."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [],
    "source": [
-    "from dolfinx.generation import UnitSquareMesh\n",
-    "from dolfinx.mesh import CellType\n",
-    "mesh = UnitSquareMesh(MPI.COMM_WORLD, 8, 8, CellType.quadrilateral)"
+    "from mpi4py import MPI\n",
+    "from dolfinx.mesh import CellType, create_unit_square\n",
+    "mesh = create_unit_square(MPI.COMM_WORLD, 8, 8, CellType.quadrilateral)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "which defines a uniform finite element mesh over the unit square $[0,1]\\times[0, 1]$. The mesh consists of _cells_, which in 2D can be either triangles of quadrilaterals. Note that in addition to give how many elements we would like to have in each direction. \n",
-    "Note that we also supply the _MPI-communicator_. \n",
-    "This is to specify how we would like the program to behave in parallel. If we supply `MPI.COMM_WORLD` we create a single mesh,\n",
-    "whose data is distributed over the number of processors we \n",
-    "would like to use. We can for instance run the program in \n",
-    "parallel on two processors by using `mpirun`, as: \n",
+    "Note that in addition to give how many elements we would like to have in each direction. \n",
+    "We also have to supply the _MPI-communicator_. \n",
+    "This is to specify how we would like the program to behave in parallel. \n",
+    "If we supply `MPI.COMM_WORLD` we create a single mesh, whose data is distributed over the number of processors we \n",
+    "would like to use. We can for instance run the program in  parallel on two processors by using `mpirun`, as: \n",
     "``` bash\n",
     " mpirun -n 2 python3 t1.py\n",
     "```\n",
-    "However, if we would like to create a separate mesh on each processor, we can use `MPI.COMM_SELF`. This is for instance \n",
-    "useful if we run a small problem, and would like to run it with\n",
-    "multiple parameters.\n",
+    "However, if we would like to create a separate mesh on each processor, we can use `MPI.COMM_SELF`.\n",
+    "This is for instance  useful if we run a small problem, and would like to run it with multiple parameters.\n",
     "\n",
     "## Defining the finite element function space\n",
-    " Once the mesh has been created, we can create the finite element function space $V$:\n"
+    " Once the mesh has been created, we can create the finite element function space $V$.\n",
+    " We import the function space initializer from the `dolfinx.fem` module."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -133,13 +119,13 @@
     "\n",
     "The next step is to specify the boundary condition $u=u_D$ on $\\partial\\Omega_D$, which is done by over several steps. \n",
     "The first step is to define the function $u_D$. Into this function, we would like to interpolate the boundary condition $1 + x^2+2y^2$.\n",
-    "We do this by first defining a `dolfinx.Function`, and then using a [lambda-function](https://docs.python.org/3/tutorial/controlflow.html#lambda-expressions) in Python to define the \n",
+    "We do this by first defining a `dolfinx.fem.Function`, and then using a [lambda-function](https://docs.python.org/3/tutorial/controlflow.html#lambda-expressions) in Python to define the \n",
     "spatially varying function.\n"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -154,21 +140,23 @@
    "source": [
     "We now have the boundary data (and in this case the solution of \n",
     "the finite element problem) represented in the discrete function space.\n",
-    "Next we would like to apply the boundary values to all degrees of freedom that are on the boundary of the discrete domain. We start by identifying the facets (line-segments) representing the outer boundary, using `dolfinx.cpp.mesh.compute_boundary_facets`.\n",
-    "This function returns an array of booleans of the same size as the number of facets on this processor, where `True` indicates that the local facet $i$ is on the boundary. To reduce this to only the indices that are `True`, we use `numpy.where`."
+    "Next we would like to apply the boundary values to all degrees of freedom that are on the boundary of the discrete domain. We start by identifying the facets (line-segments) representing the outer boundary, using `dolfinx.mesh.compute_boundary_facets`.\n",
+    "This function returns an array of booleans of the same size as the number of facets on this processor, where `True` indicates that the local facet $i$ is on the boundary. To reduce this to only the indices that are `True`, we use [`numpy.flatnonzero`](https://numpy.org/doc/stable/reference/generated/numpy.flatnonzero.html)."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [],
    "source": [
-    "fdim = mesh.topology.dim - 1\n",
+    "import numpy\n",
+    "from dolfinx.mesh import compute_boundary_facets\n",
     "# Create facet to cell connectivity required to determine boundary facets\n",
-    "from dolfinx.cpp.mesh import compute_boundary_facets\n",
-    "mesh.topology.create_connectivity(fdim, mesh.topology.dim)\n",
-    "boundary_facets = numpy.where(numpy.array(compute_boundary_facets(mesh.topology)) == 1)[0]"
+    "tdim = mesh.topology.dim\n",
+    "fdim = tdim - 1\n",
+    "mesh.topology.create_connectivity(fdim, tdim)\n",
+    "boundary_facets = numpy.flatnonzero(compute_boundary_facets(mesh.topology))"
    ]
   },
   {
@@ -187,7 +175,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 5,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -203,15 +191,18 @@
     "## Defining the trial and test function\n",
     "\n",
     "In mathematics, we distinguish between trial and test spaces $V$ and $\\hat{V}$. The only difference in the present problem is the boundary conditions.\n",
-    "In FEniCSx, we do not specify boundary conditions as part of the function space, so it is sufficient to use a common space for the trial and test function:"
+    "In FEniCSx, we do not specify boundary conditions as part of the function space, so it is sufficient to use a common space for the trial and test function.\n",
+    "\n",
+    "We use the [Unified Form Language](https://github.com/FEniCS/ufl/) (UFL) to specify the varitional formulations."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 6,
    "metadata": {},
    "outputs": [],
    "source": [
+    "import ufl\n",
     "u = ufl.TrialFunction(V)\n",
     "v = ufl.TestFunction(V)"
    ]
@@ -226,12 +217,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 11,
    "metadata": {},
    "outputs": [],
    "source": [
     "from dolfinx.fem import Constant\n",
-    "f = Constant(mesh, -6)"
+    "from petsc4py.PETSc import ScalarType\n",
+    "f = Constant(mesh, ScalarType(-6))"
    ]
   },
   {
@@ -248,7 +240,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": 12,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -267,7 +259,14 @@
     "This is the key strength of FEniCSx: the formulas in the variational formulation translate directly to very similar Python code, a feature that makes it easy to specify and solve complicated PDE problems. The language used to express weak forms is the Unified Form Language [UFL](https://doi.org/10.1145/2566630).\n",
     "\n",
     "## Expressing inner products\n",
-    "The inner product $\\int_\\Omega \\nabla u \\cdot \\nabla v ~\\mathrm{d} x$ can be expressed in various ways in UFL. We have used the notation `ufl.dot(ufl.grad(u), ufl.grad(v))*uf.dx`. The dot product in UFL computes the sum (contraction) over the last index of the first factor and first index of the second factor. In this case, both factors are tensors of rank one (vectors) and so the sum is just over the single index of both $\\nabla u$ and $\\nabla v$. To compute an inner product of matrices (with two indices), one must instead of `ufl.dot` use the function `ufl.inner`. For vectors, `ufl.dot` and `ufl.inner` are equivalent.\n",
+    "The inner product $\\int_\\Omega \\nabla u \\cdot \\nabla v ~\\mathrm{d} x$ can be expressed in various ways in UFL. We have used the notation `ufl.dot(ufl.grad(u), ufl.grad(v))*ufl.dx`. The dot product in UFL computes the sum (contraction) over the last index of the first factor and first index of the second factor. In this case, both factors are tensors of rank one (vectors) and so the sum is just over the single index of both $\\nabla u$ and $\\nabla v$. To compute an inner product of matrices (with two indices), one must instead of `ufl.dot` use the function `ufl.inner`. For vectors, `ufl.dot` and `ufl.inner` are equivalent.\n",
+    "\n",
+    "```{admonition} Complex numbers\n",
+    "In DOLFINx, one can solve complex number problems by using an installation of PETSc using complex numbers.\n",
+    "For variational formulations with complex numbers, one cannot use `ufl.dot` to compute inner products.\n",
+    "One has to use `ufl.inner`, with the test-function as the second input argument for `ufl.inner`.\n",
+    "```\n",
+    "\n",
     "\n",
     "## Forming and solving the linear system\n",
     "\n",
@@ -277,7 +276,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": 13,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -297,7 +296,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": 14,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -310,7 +309,8 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "We compute the error in two different ways. First, we compute the $L^2$-norm of the error, defined by $E=\\sqrt{\\int_\\Omega (u_D-u_h)^2\\mathrm{d} x}$. We use UFL to express the $L^2$-error, and use `dolfinx.fem.assemble_scalar` to compute the scalar value."
+    "We compute the error in two different ways. First, we compute the $L^2$-norm of the error, defined by $E=\\sqrt{\\int_\\Omega (u_D-u_h)^2\\mathrm{d} x}$. We use UFL to express the $L^2$-error, and use `dolfinx.fem.assemble_scalar` to compute the scalar value. In DOLFINx, `assemble_scalar` only assembles over the cells on the local process. This means that if we use 2 processes to solve our problem, we need to gather the solution to one (or all the processes.\n",
+    "We can do this with the `MPI.allreduce` function."
    ]
   },
   {
@@ -321,7 +321,8 @@
    "source": [
     "from dolfinx.fem import assemble_scalar\n",
     "L2_error = ufl.inner(uh - uex, uh - uex) * ufl.dx\n",
-    "error_L2 = numpy.sqrt(assemble_scalar(L2_error))"
+    "error_local = assemble_scalar(L2_error)\n",
+    "error_L2 = numpy.sqrt(mesh.comm.allreduce(error_local, op=MPI.SUM))"
    ]
   },
   {
@@ -338,24 +339,26 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 17,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
       "Error_L2 : 8.24e-03\n",
-      "Error_max : 3.11e-15\n"
+      "Error_max : 2.22e-15\n"
      ]
     }
    ],
    "source": [
     "u_vertex_values = uh.compute_point_values()\n",
     "u_ex_vertex_values = uex.compute_point_values()\n",
     "error_max = numpy.max(numpy.abs(u_vertex_values - u_ex_vertex_values))\n",
-    "print(f\"Error_L2 : {error_L2:.2e}\")\n",
-    "print(f\"Error_max : {error_max:.2e}\")"
+    "# Only print the error on one process\n",
+    "if mesh.comm.rank == 0:\n",
+    "    print(f\"Error_L2 : {error_L2:.2e}\")\n",
+    "    print(f\"Error_max : {error_max:.2e}\")"
    ]
   },
   {
@@ -370,17 +373,17 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": 22,
    "metadata": {},
    "outputs": [],
    "source": [
-    "import dolfinx.plot\n",
-    "topology, cell_types = dolfinx.plot.create_vtk_topology(mesh, mesh.topology.dim)"
+    "from dolfinx.plot import create_vtk_topology\n",
+    "topology, cell_types = create_vtk_topology(mesh, mesh.topology.dim)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 18,
+   "execution_count": 23,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -397,7 +400,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 19,
+   "execution_count": 24,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -418,13 +421,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 20,
+   "execution_count": 26,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "c8d829e8954c46699aca8c42cd8b5b12",
+       "model_id": "cee1088755b04b99bdedbc3bf73335a3",
        "version_major": 2,
        "version_minor": 0
       },
@@ -457,13 +460,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 21,
+   "execution_count": 27,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "5a5bdbffdc324a548691f67bb09a5a4c",
+       "model_id": "9d3ee0699fb744d8b619d2cf7f4af5d2",
        "version_major": 2,
        "version_minor": 0
       },
@@ -493,7 +496,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 22,
+   "execution_count": 28,
    "metadata": {},
    "outputs": [],
    "source": [