diff --git a/_config.yml b/_config.yml index 2f4c9428..32047621 100644 --- a/_config.yml +++ b/_config.yml @@ -42,19 +42,20 @@ sphinx: # navigation_with_keys: false codeautolink_concat_default: True intersphinx_mapping: - basix : ["https://docs.fenicsproject.org/basix/main/python/", null] - ffcx : ["https://docs.fenicsproject.org/ffcx/main/", null] - ufl : ["https://docs.fenicsproject.org/ufl/main/", null] - dolfinx : ["https://docs.fenicsproject.org/dolfinx/main/python", null] - petsc4py: ["https://petsc.org/release/petsc4py", null] - mpi4py: ["https://mpi4py.readthedocs.io/en/stable", null] - numpy: ["https://numpy.org/doc/stable/", null] - pyvista: ["https://docs.pyvista.org/", null] - packaging: ["https://packaging.pypa.io/en/stable/", null] + basix: ["https://docs.fenicsproject.org/basix/main/python/", null] + ffcx: ["https://docs.fenicsproject.org/ffcx/main/", null] + ufl: ["https://docs.fenicsproject.org/ufl/main/", null] + dolfinx: ["https://docs.fenicsproject.org/dolfinx/main/python", null] + petsc4py: ["https://petsc.org/release/petsc4py", null] + mpi4py: ["https://mpi4py.readthedocs.io/en/stable", null] + numpy: ["https://numpy.org/doc/stable/", null] + pyvista: ["https://docs.pyvista.org/", null] + packaging: ["https://packaging.pypa.io/en/stable/", null] + matplotlib: ["https://matplotlib.org/stable/", null] extra_extensions: - - 'sphinx.ext.autodoc' - - 'sphinx.ext.intersphinx' + - "sphinx.ext.autodoc" + - "sphinx.ext.intersphinx" - "sphinx_codeautolink" parse: diff --git a/chapter1/membrane_code.ipynb b/chapter1/membrane_code.ipynb index f5f7d4cf..fcc6945a 100644 --- a/chapter1/membrane_code.ipynb +++ b/chapter1/membrane_code.ipynb @@ -11,15 +11,16 @@ "In this section, we will solve the deflection of the membrane problem.\n", "After finishing this section, you should be able to:\n", "- Create a simple mesh using the GMSH Python API and load it into DOLFINx\n", - "- Create constant boundary conditions using a geometrical identifier\n", - "- Use `ufl.SpatialCoordinate` to create a spatially varying function\n", - "- Interpolate a `ufl.Expression` into an appropriate function space\n", - "- Evaluate a `dolfinx.fem.Function` at any point $x$\n", + "- Create constant boundary conditions using a {py:func}`geometrical identifier`\n", + "- Use {py:class}`ufl.SpatialCoordinate` to create a spatially varying function\n", + "- Interpolate a {py:class}`ufl-Expression` into an appropriate function space\n", + "- Evaluate a {py:class}`dolfinx.fem.Function` at any point $x$\n", "- Use Paraview to visualize the solution of a PDE\n", "\n", "## Creating the mesh\n", "\n", - "To create the computational geometry, we use the Python-API of [GMSH](https://gmsh.info/). We start by importing the gmsh-module and initializing it." + "To create the computational geometry, we use the Python-API of [GMSH](https://gmsh.info/).\n", + "We start by importing the gmsh-module and initializing it." ] }, { @@ -39,7 +40,10 @@ "id": "2", "metadata": {}, "source": [ - "The next step is to create the membrane and start the computations by the GMSH CAD kernel, to generate the relevant underlying data structures. The first arguments of `addDisk` are the x, y and z coordinate of the center of the circle, while the two last arguments are the x-radius and y-radius." + "The next step is to create the membrane and start the computations by the GMSH CAD kernel,\n", + "to generate the relevant underlying data structures.\n", + "The first arguments of `addDisk` are the x, y and z coordinate of the center of the circle,\n", + "while the two last arguments are the x-radius and y-radius." ] }, { @@ -58,7 +62,10 @@ "id": "4", "metadata": {}, "source": [ - "After that, we make the membrane a physical surface, such that it is recognized by `gmsh` when generating the mesh. As a surface is a two-dimensional entity, we add `2` as the first argument, the entity tag of the membrane as the second argument, and the physical tag as the last argument. In a later demo, we will get into when this tag matters." + "After that, we make the membrane a physical surface, such that it is recognized by `gmsh` when generating the mesh.\n", + "As a surface is a two-dimensional entity, we add `2` as the first argument,\n", + "the entity tag of the membrane as the second argument, and the physical tag as the last argument.\n", + "In a later demo, we will get into when this tag matters." ] }, { @@ -77,7 +84,8 @@ "id": "6", "metadata": {}, "source": [ - "Finally, we generate the two-dimensional mesh. We set a uniform mesh size by modifying the GMSH options." + "Finally, we generate the two-dimensional mesh.\n", + "We set a uniform mesh size by modifying the GMSH options." ] }, { @@ -98,25 +106,29 @@ "metadata": {}, "source": [ "# Interfacing with GMSH in DOLFINx\n", - "We will import the GMSH-mesh directly from GMSH into DOLFINx via the `dolfinx.io.gmsh` interface.\n", + "We will import the GMSH-mesh directly from GMSH into DOLFINx via the {py:mod}`dolfinx.io.gmsh` interface.\n", "The {py:mod}`dolfinx.io.gmsh` module contains two functions\n", "1. {py:func}`model_to_mesh` which takes in a `gmsh.model`\n", " and returns a {py:class}`dolfinx.io.gmsh.MeshData` object.\n", "2. {py:func}`read_from_msh` which takes in a path to a `.msh`-file\n", " and returns a {py:class}`dolfinx.io.gmsh.MeshData` object.\n", "\n", - "The `MeshData` object will contain a `dolfinx.mesh.Mesh`, under the attribute `mesh`.\n", + "The {py:class}`MeshData` object will contain a {py:class}`dolfinx.mesh.Mesh`,\n", + "under the attribute {py:attr}`mesh`.\n", "This mesh will contain all GMSH Physical Groups of the highest topological dimension.\n", "```{note}\n", "If you do not use `gmsh.model.addPhysicalGroup` when creating the mesh with GMSH, it can not be read into DOLFINx.\n", "```\n", - "The `MeshData` object can also contain tags for all other `PhysicalGroups` that has been added to the mesh,\n", - "that being `vertex_tags`, `edge_tags`, `facet_tags` and `cell_tags`.\n", + "The {py:class}`MeshData` object can also contain tags for\n", + "all other `PhysicalGroups` that has been added to the mesh, that being\n", + "{py:attr}`cell_tags`, {py:attr}`facet_tags`,\n", + "{py:attr}`ridge_tags` and\n", + "{py:attr}`peak_tags`.\n", "To read either `gmsh.model` or a `.msh`-file, one has to distribute the mesh to all processes used by DOLFINx.\n", "As GMSH does not support mesh creation with MPI, we currently have a `gmsh.model.mesh` on each process.\n", "To distribute the mesh, we have to specify which process the mesh was created on,\n", "and which communicator rank should distribute the mesh.\n", - "The `model_to_mesh` will then load the mesh on the specified rank,\n", + "The {py:func}`model_to_mesh` will then load the mesh on the specified rank,\n", "and distribute it to the communicator using a mesh partitioner." ] }, @@ -165,7 +177,8 @@ "metadata": {}, "source": [ "## Defining a spatially varying load\n", - "The right hand side pressure function is represented using `ufl.SpatialCoordinate` and two constants, one for $\\beta$ and one for $R_0$." + "The right hand side pressure function is represented using {py:class}`ufl.SpatialCoordinate` and two constants,\n", + "one for $\\beta$ and one for $R_0$." ] }, { @@ -176,16 +189,8 @@ "outputs": [], "source": [ "import ufl\n", - "from dolfinx import default_scalar_type" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "14", - "metadata": {}, - "outputs": [], - "source": [ + "from dolfinx import default_scalar_type\n", + "\n", "x = ufl.SpatialCoordinate(domain)\n", "beta = fem.Constant(domain, default_scalar_type(12))\n", "R0 = fem.Constant(domain, default_scalar_type(0.3))\n", @@ -194,11 +199,15 @@ }, { "cell_type": "markdown", - "id": "15", + "id": "9f059bd4", "metadata": {}, "source": [ "## Create a Dirichlet boundary condition using geometrical conditions\n", - "The next step is to create the homogeneous boundary condition. As opposed to the [first tutorial](./fundamentals_code.ipynb) we will use `dolfinx.fem.locate_dofs_geometrical` to locate the degrees of freedom on the boundary. As we know that our domain is a circle with radius 1, we know that any degree of freedom should be located at a coordinate $(x,y)$ such that $\\sqrt{x^2+y^2}=1$." + "The next step is to create the homogeneous boundary condition.\n", + "As opposed to the [first tutorial](./fundamentals_code.ipynb) we will use\n", + "{py:func}`locate_dofs_geometrical` to locate the degrees of freedom on the boundary.\n", + "As we know that our domain is a circle with radius 1, we know that any degree of freedom should be\n", + "located at a coordinate $(x,y)$ such that $\\sqrt{x^2+y^2}=1$." ] }, { @@ -223,7 +232,9 @@ "id": "17", "metadata": {}, "source": [ - "As our Dirichlet condition is homogeneous (`u=0` on the whole boundary), we can initialize the `dolfinx.fem.dirichletbc` with a constant value, the degrees of freedom and the function space to apply the boundary condition on." + "As our Dirichlet condition is homogeneous (`u=0` on the whole boundary), we can initialize the\n", + "{py:class}`dolfinx.fem.DirichletBC` with a constant value, the degrees of freedom and the function\n", + "space to apply the boundary condition on. We use the constructor {py:func}`dolfinx.fem.dirichletbc`." ] }, { @@ -271,8 +282,13 @@ "id": "21", "metadata": {}, "source": [ - "## Interpolation of a `ufl`-expression\n", - "As we previously defined the load `p` as a spatially varying function, we would like to interpolate this function into an appropriate function space for visualization. To do this we use the `dolfinx.Expression`. The expression takes in any `ufl`-expression, and a set of points on the reference element. We will use the interpolation points of the space we want to interpolate in to.\n", + "## Interpolation of a UFL-expression\n", + "As we previously defined the load `p` as a spatially varying function,\n", + "we would like to interpolate this function into an appropriate function space for visualization.\n", + "To do this we use the class {py:class}`Expression`.\n", + "The expression takes in any UFL-expression, and a set of points on the reference element.\n", + "We will use the {py:attr}`interpolation points`\n", + "of the space we want to interpolate in to.\n", "We choose a high order function space to represent the function `p`, as it is rapidly varying in space." ] }, @@ -306,15 +322,15 @@ "outputs": [], "source": [ "from dolfinx.plot import vtk_mesh\n", - "import pyvista\n" + "import pyvista" ] }, { "cell_type": "markdown", - "id": "25", + "id": "55b28bcd", "metadata": {}, "source": [ - "Extract topology from mesh and create pyvista mesh" + "Extract topology from mesh and create {py:class}`pyvista.UnstructuredGrid`" ] }, { @@ -330,7 +346,7 @@ }, { "cell_type": "markdown", - "id": "27", + "id": "b2b183e2", "metadata": {}, "source": [ "Set deflection values and add it to plotter" @@ -389,7 +405,8 @@ "source": [ "## Making curve plots throughout the domain\n", "Another way to compare the deflection and the load is to make a plot along the line $x=0$.\n", - "This is just a matter of defining a set of points along the $y$-axis and evaluating the finite element functions $u$ and $p$ at these points." + "This is just a matter of defining a set of points along the $y$-axis and evaluating the\n", + "finite element functions $u$ and $p$ at these points." ] }, { @@ -412,9 +429,15 @@ "id": "33", "metadata": {}, "source": [ - "As a finite element function is the linear combination of all degrees of freedom, $u_h(x)=\\sum_{i=1}^N c_i \\phi_i(x)$ where $c_i$ are the coefficients of $u_h$ and $\\phi_i$ is the $i$-th basis function, we can compute the exact solution at any point in $\\Omega$.\n", - "However, as a mesh consists of a large set of degrees of freedom (i.e. $N$ is large), we want to reduce the number of evaluations of the basis function $\\phi_i(x)$. We do this by identifying which cell of the mesh $x$ is in.\n", - "This is efficiently done by creating a bounding box tree of the cells of the mesh, allowing a quick recursive search through the mesh entities." + "As a finite element function is the linear combination of all degrees of freedom,\n", + "$u_h(x)=\\sum_{i=1}^N c_i \\phi_i(x)$ where $c_i$ are the coefficients of $u_h$ and $\\phi_i$\n", + "is the $i$-th basis function, we can compute the exact solution at any point in $\\Omega$.\n", + "However, as a mesh consists of a large set of degrees of freedom (i.e. $N$ is large),\n", + "we want to reduce the number of evaluations of the basis function $\\phi_i(x)$.\n", + "We do this by identifying which cell of the mesh $x$ is in.\n", + "This is efficiently done by creating a {py:class}`bounding box tree`.\n", + "However, as the bounding box of a cell spans more of $\\mathbb{R}^n$ than the actual cell,\n", + "we check that the actual cell collides with the input point using\n", + "{py:func}`dolfinx.geometry.compute_colliding_cells`,\n", + "which measures the exact distance between the point and the cell\n", + "(approximated as a convex hull for higher order geometries).\n", + "This function also returns an adjacency-list, as the point might align with a facet,\n", + "edge or vertex that is shared between multiple cells in the mesh.\n", "\n", - "Finally, we would like the code below to run in parallel, when the mesh is distributed over multiple processors. In that case, it is not guaranteed that every point in `points` is on each processor. Therefore we create a subset `points_on_proc` only containing the points found on the current processor." + "Finally, we would like the code below to run in parallel,\n", + "when the mesh is distributed over multiple processors.\n", + "In that case, it is not guaranteed that every point in `points` is on each processor.\n", + "Therefore we create a subset `points_on_proc` only containing the points found on the current processor." ] }, { @@ -466,7 +501,9 @@ "id": "37", "metadata": {}, "source": [ - "We now have a list of points on the processor, on in which cell each point belongs. We can then call `uh.eval` and `pressure.eval` to obtain the set of values for all the points.\n" + "We now have a list of points on the processor, on in which cell each point belongs.\n", + "We can then call {py:meth}`uh.eval` and\n", + "{py:meth}`pressure.eval` to obtain the set of values for all the points." ] }, { @@ -486,7 +523,8 @@ "id": "39", "metadata": {}, "source": [ - "As we now have an array of coordinates and two arrays of function values, we can use `matplotlib` to plot them" + "As we now have an array of coordinates and two arrays of function values,\n", + "we can use {py:mod}`matplotlib` to plot them" ] }, { @@ -517,7 +555,7 @@ "id": "41", "metadata": {}, "source": [ - "If executed in parallel as a python file, we save a plot per processor" + "If executed in parallel as a Python file, we save a plot per processor" ] }, { diff --git a/chapter1/membrane_code.py b/chapter1/membrane_code.py index 0c4a5b3d..09a58d8f 100644 --- a/chapter1/membrane_code.py +++ b/chapter1/membrane_code.py @@ -19,15 +19,16 @@ # In this section, we will solve the deflection of the membrane problem. # After finishing this section, you should be able to: # - Create a simple mesh using the GMSH Python API and load it into DOLFINx -# - Create constant boundary conditions using a geometrical identifier -# - Use `ufl.SpatialCoordinate` to create a spatially varying function -# - Interpolate a `ufl.Expression` into an appropriate function space -# - Evaluate a `dolfinx.fem.Function` at any point $x$ +# - Create constant boundary conditions using a {py:func}`geometrical identifier` +# - Use {py:class}`ufl.SpatialCoordinate` to create a spatially varying function +# - Interpolate a {py:class}`ufl-Expression` into an appropriate function space +# - Evaluate a {py:class}`dolfinx.fem.Function` at any point $x$ # - Use Paraview to visualize the solution of a PDE # # ## Creating the mesh # -# To create the computational geometry, we use the Python-API of [GMSH](https://gmsh.info/). We start by importing the gmsh-module and initializing it. +# To create the computational geometry, we use the Python-API of [GMSH](https://gmsh.info/). +# We start by importing the gmsh-module and initializing it. # + import gmsh @@ -35,42 +36,53 @@ gmsh.initialize() # - -# The next step is to create the membrane and start the computations by the GMSH CAD kernel, to generate the relevant underlying data structures. The first arguments of `addDisk` are the x, y and z coordinate of the center of the circle, while the two last arguments are the x-radius and y-radius. +# The next step is to create the membrane and start the computations by the GMSH CAD kernel, +# to generate the relevant underlying data structures. +# The first arguments of `addDisk` are the x, y and z coordinate of the center of the circle, +# while the two last arguments are the x-radius and y-radius. membrane = gmsh.model.occ.addDisk(0, 0, 0, 1, 1) gmsh.model.occ.synchronize() -# After that, we make the membrane a physical surface, such that it is recognized by `gmsh` when generating the mesh. As a surface is a two-dimensional entity, we add `2` as the first argument, the entity tag of the membrane as the second argument, and the physical tag as the last argument. In a later demo, we will get into when this tag matters. +# After that, we make the membrane a physical surface, such that it is recognized by `gmsh` when generating the mesh. +# As a surface is a two-dimensional entity, we add `2` as the first argument, +# the entity tag of the membrane as the second argument, and the physical tag as the last argument. +# In a later demo, we will get into when this tag matters. gdim = 2 gmsh.model.addPhysicalGroup(gdim, [membrane], 1) -# Finally, we generate the two-dimensional mesh. We set a uniform mesh size by modifying the GMSH options. +# Finally, we generate the two-dimensional mesh. +# We set a uniform mesh size by modifying the GMSH options. gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.05) gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.05) gmsh.model.mesh.generate(gdim) # # Interfacing with GMSH in DOLFINx -# We will import the GMSH-mesh directly from GMSH into DOLFINx via the `dolfinx.io.gmsh` interface. +# We will import the GMSH-mesh directly from GMSH into DOLFINx via the {py:mod}`dolfinx.io.gmsh` interface. # The {py:mod}`dolfinx.io.gmsh` module contains two functions # 1. {py:func}`model_to_mesh` which takes in a `gmsh.model` # and returns a {py:class}`dolfinx.io.gmsh.MeshData` object. # 2. {py:func}`read_from_msh` which takes in a path to a `.msh`-file # and returns a {py:class}`dolfinx.io.gmsh.MeshData` object. # -# The `MeshData` object will contain a `dolfinx.mesh.Mesh`, under the attribute `mesh`. +# The {py:class}`MeshData` object will contain a {py:class}`dolfinx.mesh.Mesh`, +# under the attribute {py:attr}`mesh`. # This mesh will contain all GMSH Physical Groups of the highest topological dimension. # ```{note} # If you do not use `gmsh.model.addPhysicalGroup` when creating the mesh with GMSH, it can not be read into DOLFINx. # ``` -# The `MeshData` object can also contain tags for all other `PhysicalGroups` that has been added to the mesh, -# that being `vertex_tags`, `edge_tags`, `facet_tags` and `cell_tags`. +# The {py:class}`MeshData` object can also contain tags for +# all other `PhysicalGroups` that has been added to the mesh, that being +# {py:attr}`cell_tags`, {py:attr}`facet_tags`, +# {py:attr}`ridge_tags` and +# {py:attr}`peak_tags`. # To read either `gmsh.model` or a `.msh`-file, one has to distribute the mesh to all processes used by DOLFINx. # As GMSH does not support mesh creation with MPI, we currently have a `gmsh.model.mesh` on each process. # To distribute the mesh, we have to specify which process the mesh was created on, # and which communicator rank should distribute the mesh. -# The `model_to_mesh` will then load the mesh on the specified rank, +# The {py:func}`model_to_mesh` will then load the mesh on the specified rank, # and distribute it to the communicator using a mesh partitioner. # + @@ -95,8 +107,10 @@ # - # ## Defining a spatially varying load -# The right hand side pressure function is represented using `ufl.SpatialCoordinate` and two constants, one for $\beta$ and one for $R_0$. +# The right hand side pressure function is represented using {py:class}`ufl.SpatialCoordinate` and two constants, +# one for $\beta$ and one for $R_0$. +# + import ufl from dolfinx import default_scalar_type @@ -104,9 +118,14 @@ beta = fem.Constant(domain, default_scalar_type(12)) R0 = fem.Constant(domain, default_scalar_type(0.3)) p = 4 * ufl.exp(-(beta**2) * (x[0] ** 2 + (x[1] - R0) ** 2)) +# - # ## Create a Dirichlet boundary condition using geometrical conditions -# The next step is to create the homogeneous boundary condition. As opposed to the [first tutorial](./fundamentals_code.ipynb) we will use `dolfinx.fem.locate_dofs_geometrical` to locate the degrees of freedom on the boundary. As we know that our domain is a circle with radius 1, we know that any degree of freedom should be located at a coordinate $(x,y)$ such that $\sqrt{x^2+y^2}=1$. +# The next step is to create the homogeneous boundary condition. +# As opposed to the [first tutorial](./fundamentals_code.ipynb) we will use +# {py:func}`locate_dofs_geometrical` to locate the degrees of freedom on the boundary. +# As we know that our domain is a circle with radius 1, we know that any degree of freedom should be +# located at a coordinate $(x,y)$ such that $\sqrt{x^2+y^2}=1$. # + import numpy as np @@ -119,7 +138,9 @@ def on_boundary(x): boundary_dofs = fem.locate_dofs_geometrical(V, on_boundary) # - -# As our Dirichlet condition is homogeneous (`u=0` on the whole boundary), we can initialize the `dolfinx.fem.dirichletbc` with a constant value, the degrees of freedom and the function space to apply the boundary condition on. +# As our Dirichlet condition is homogeneous (`u=0` on the whole boundary), we can initialize the +# {py:class}`dolfinx.fem.DirichletBC` with a constant value, the degrees of freedom and the function +# space to apply the boundary condition on. We use the constructor {py:func}`dolfinx.fem.dirichletbc`. bc = fem.dirichletbc(default_scalar_type(0), boundary_dofs, V) @@ -139,8 +160,13 @@ def on_boundary(x): ) uh = problem.solve() -# ## Interpolation of a `ufl`-expression -# As we previously defined the load `p` as a spatially varying function, we would like to interpolate this function into an appropriate function space for visualization. To do this we use the `dolfinx.Expression`. The expression takes in any `ufl`-expression, and a set of points on the reference element. We will use the interpolation points of the space we want to interpolate in to. +# ## Interpolation of a UFL-expression +# As we previously defined the load `p` as a spatially varying function, +# we would like to interpolate this function into an appropriate function space for visualization. +# To do this we use the class {py:class}`Expression`. +# The expression takes in any UFL-expression, and a set of points on the reference element. +# We will use the {py:attr}`interpolation points` +# of the space we want to interpolate in to. # We choose a high order function space to represent the function `p`, as it is rapidly varying in space. Q = fem.functionspace(domain, ("Lagrange", 5)) @@ -154,8 +180,7 @@ def on_boundary(x): from dolfinx.plot import vtk_mesh import pyvista - -# Extract topology from mesh and create pyvista mesh +# Extract topology from mesh and create {py:class}`pyvista.UnstructuredGrid` topology, cell_types, x = vtk_mesh(V) grid = pyvista.UnstructuredGrid(topology, cell_types, x) @@ -190,7 +215,8 @@ def on_boundary(x): # ## Making curve plots throughout the domain # Another way to compare the deflection and the load is to make a plot along the line $x=0$. -# This is just a matter of defining a set of points along the $y$-axis and evaluating the finite element functions $u$ and $p$ at these points. +# This is just a matter of defining a set of points along the $y$-axis and evaluating the +# finite element functions $u$ and $p$ at these points. tol = 0.001 # Avoid hitting the outside of the domain y = np.linspace(-1 + tol, 1 - tol, 101) @@ -199,9 +225,15 @@ def on_boundary(x): u_values = [] p_values = [] -# As a finite element function is the linear combination of all degrees of freedom, $u_h(x)=\sum_{i=1}^N c_i \phi_i(x)$ where $c_i$ are the coefficients of $u_h$ and $\phi_i$ is the $i$-th basis function, we can compute the exact solution at any point in $\Omega$. -# However, as a mesh consists of a large set of degrees of freedom (i.e. $N$ is large), we want to reduce the number of evaluations of the basis function $\phi_i(x)$. We do this by identifying which cell of the mesh $x$ is in. -# This is efficiently done by creating a bounding box tree of the cells of the mesh, allowing a quick recursive search through the mesh entities. +# As a finite element function is the linear combination of all degrees of freedom, +# $u_h(x)=\sum_{i=1}^N c_i \phi_i(x)$ where $c_i$ are the coefficients of $u_h$ and $\phi_i$ +# is the $i$-th basis function, we can compute the exact solution at any point in $\Omega$. +# However, as a mesh consists of a large set of degrees of freedom (i.e. $N$ is large), +# we want to reduce the number of evaluations of the basis function $\phi_i(x)$. +# We do this by identifying which cell of the mesh $x$ is in. +# This is efficiently done by creating a {py:class}`bounding box tree`. +# However, as the bounding box of a cell spans more of $\mathbb{R}^n$ than the actual cell, +# we check that the actual cell collides with the input point using +# {py:func}`dolfinx.geometry.compute_colliding_cells`, +# which measures the exact distance between the point and the cell +# (approximated as a convex hull for higher order geometries). +# This function also returns an adjacency-list, as the point might align with a facet, +# edge or vertex that is shared between multiple cells in the mesh. # -# Finally, we would like the code below to run in parallel, when the mesh is distributed over multiple processors. In that case, it is not guaranteed that every point in `points` is on each processor. Therefore we create a subset `points_on_proc` only containing the points found on the current processor. +# Finally, we would like the code below to run in parallel, +# when the mesh is distributed over multiple processors. +# In that case, it is not guaranteed that every point in `points` is on each processor. +# Therefore we create a subset `points_on_proc` only containing the points found on the current processor. cells = [] points_on_proc = [] @@ -227,14 +271,16 @@ def on_boundary(x): points_on_proc.append(point) cells.append(colliding_cells.links(i)[0]) -# We now have a list of points on the processor, on in which cell each point belongs. We can then call `uh.eval` and `pressure.eval` to obtain the set of values for all the points. -# +# We now have a list of points on the processor, on in which cell each point belongs. +# We can then call {py:meth}`uh.eval` and +# {py:meth}`pressure.eval` to obtain the set of values for all the points. points_on_proc = np.array(points_on_proc, dtype=np.float64) u_values = uh.eval(points_on_proc, cells) p_values = pressure.eval(points_on_proc, cells) -# As we now have an array of coordinates and two arrays of function values, we can use `matplotlib` to plot them +# As we now have an array of coordinates and two arrays of function values, +# we can use {py:mod}`matplotlib` to plot them # + import matplotlib.pyplot as plt @@ -253,7 +299,7 @@ def on_boundary(x): plt.legend() # - -# If executed in parallel as a python file, we save a plot per processor +# If executed in parallel as a Python file, we save a plot per processor plt.savefig(f"membrane_rank{MPI.COMM_WORLD.rank:d}.png") diff --git a/chapter2/amr.ipynb b/chapter2/amr.ipynb index cb0c05af..2ec33082 100644 --- a/chapter2/amr.ipynb +++ b/chapter2/amr.ipynb @@ -138,7 +138,6 @@ }, "outputs": [], "source": [ - "\n", "grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh))\n", "grid.cell_data[\"ct\"] = ct.values\n", "\n", diff --git a/chapter2/amr.py b/chapter2/amr.py index 005b66e0..0ed1b555 100644 --- a/chapter2/amr.py +++ b/chapter2/amr.py @@ -72,7 +72,6 @@ # We use pyvista to visualize the mesh. # + tags=["hide-input"] - grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh)) grid.cell_data["ct"] = ct.values diff --git a/chapter2/helmholtz_code.ipynb b/chapter2/helmholtz_code.ipynb index 1bab721e..be8527f7 100644 --- a/chapter2/helmholtz_code.ipynb +++ b/chapter2/helmholtz_code.ipynb @@ -95,7 +95,7 @@ " assert mesh_data.facet_tags is not None\n", " facet_tags = mesh_data.facet_tags\n", "else:\n", - " domain, _, facet_tags = mesh_data\n" + " domain, _, facet_tags = mesh_data" ] }, { @@ -113,7 +113,6 @@ "metadata": {}, "outputs": [], "source": [ - "\n", "V = fem.functionspace(domain, (\"Lagrange\", 1))\n", "\n", "# Discrete frequency range\n", @@ -121,7 +120,7 @@ "\n", "# Air parameters\n", "rho0 = 1.225 # kg/m^3\n", - "c = 340 # m/s\n" + "c = 340 # m/s" ] }, { diff --git a/chapter2/helmholtz_code.py b/chapter2/helmholtz_code.py index 0fa0b03f..7122998e 100644 --- a/chapter2/helmholtz_code.py +++ b/chapter2/helmholtz_code.py @@ -89,13 +89,11 @@ facet_tags = mesh_data.facet_tags else: domain, _, facet_tags = mesh_data - # - # We define the function space for our unknown $p$ and define the range of frequencies we want to solve the Helmholtz equation for. # + - V = fem.functionspace(domain, ("Lagrange", 1)) # Discrete frequency range @@ -104,7 +102,6 @@ # Air parameters rho0 = 1.225 # kg/m^3 c = 340 # m/s - # - # ## Boundary conditions diff --git a/chapter2/nonlinpoisson_code.ipynb b/chapter2/nonlinpoisson_code.ipynb index 2a66c7fd..b29362aa 100644 --- a/chapter2/nonlinpoisson_code.ipynb +++ b/chapter2/nonlinpoisson_code.ipynb @@ -28,9 +28,8 @@ "import numpy\n", "\n", "from mpi4py import MPI\n", - "from petsc4py import PETSc\n", "\n", - "from dolfinx import mesh, fem, log\n", + "from dolfinx import mesh, fem\n", "from dolfinx.fem.petsc import NonlinearProblem\n", "\n", "\n", diff --git a/chapter2/nonlinpoisson_code.py b/chapter2/nonlinpoisson_code.py index db11896c..0b241c25 100644 --- a/chapter2/nonlinpoisson_code.py +++ b/chapter2/nonlinpoisson_code.py @@ -31,9 +31,8 @@ import numpy from mpi4py import MPI -from petsc4py import PETSc -from dolfinx import mesh, fem, log +from dolfinx import mesh, fem from dolfinx.fem.petsc import NonlinearProblem diff --git a/chapter2/ns_code1.ipynb b/chapter2/ns_code1.ipynb index 3bff8e41..2ca7bdc1 100644 --- a/chapter2/ns_code1.ipynb +++ b/chapter2/ns_code1.ipynb @@ -501,7 +501,9 @@ "cell_type": "code", "execution_count": null, "id": "26", - "metadata": {}, + "metadata": { + "lines_to_end_of_cell_marker": 2 + }, "outputs": [], "source": [ "from pathlib import Path\n", @@ -516,7 +518,9 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "lines_to_next_cell": 2 + }, "source": [ "We also interpolate the analytical solution into our function-space and create a variational formulation for the $L^2$-error.\n" ] @@ -527,8 +531,6 @@ "metadata": {}, "outputs": [], "source": [ - "\n", - "\n", "def u_exact(x):\n", " values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)\n", " values[0] = 4 * x[1] * (1.0 - x[1])\n", diff --git a/chapter2/ns_code1.py b/chapter2/ns_code1.py index 0b71f369..627a4e30 100644 --- a/chapter2/ns_code1.py +++ b/chapter2/ns_code1.py @@ -364,14 +364,15 @@ def sigma(u, p): vtx_p = VTXWriter(mesh.comm, folder / "poiseuille_p.bp", p_n, engine="BP4") vtx_u.write(t) vtx_p.write(t) + + # - # We also interpolate the analytical solution into our function-space and create a variational formulation for the $L^2$-error. # -# + - +# + def u_exact(x): values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType) values[0] = 4 * x[1] * (1.0 - x[1]) diff --git a/chapter3/component_bc.ipynb b/chapter3/component_bc.ipynb index bc2c5de2..2036e62d 100644 --- a/chapter3/component_bc.ipynb +++ b/chapter3/component_bc.ipynb @@ -65,7 +65,6 @@ "from dolfinx.fem import (\n", " Constant,\n", " dirichletbc,\n", - " Function,\n", " functionspace,\n", " locate_dofs_geometrical,\n", " locate_dofs_topological,\n", diff --git a/chapter3/component_bc.py b/chapter3/component_bc.py index 8b6bd726..5787ad54 100644 --- a/chapter3/component_bc.py +++ b/chapter3/component_bc.py @@ -66,7 +66,6 @@ from dolfinx.fem import ( Constant, dirichletbc, - Function, functionspace, locate_dofs_geometrical, locate_dofs_topological, diff --git a/chapter3/em.ipynb b/chapter3/em.ipynb index d286266a..4a983877 100644 --- a/chapter3/em.ipynb +++ b/chapter3/em.ipynb @@ -302,7 +302,6 @@ }, "outputs": [], "source": [ - "\n", "Q = functionspace(mesh, (\"DG\", 0))\n", "material_tags = np.unique(ct.values)\n", "mu = Function(Q)\n", diff --git a/chapter3/em.py b/chapter3/em.py index e04c41a8..6e9f9e90 100644 --- a/chapter3/em.py +++ b/chapter3/em.py @@ -244,8 +244,6 @@ # Next, we define the discontinous functions for the permeability $\mu$ and current $J_z$ using the `MeshTags` as in [Defining material parameters through subdomains](./subdomains) # -# + - Q = functionspace(mesh, ("DG", 0)) material_tags = np.unique(ct.values) mu = Function(Q) @@ -266,7 +264,6 @@ J.x.array[cells] = np.full_like(cells, 1, dtype=default_scalar_type) elif tag in range(2 + N, 2 * N + 2): J.x.array[cells] = np.full_like(cells, -1, dtype=default_scalar_type) -# - # In the code above, we have used a somewhat less extreme value for the magnetic permability of iron. This is to make the solution a little more interesting. It would otherwise be completely dominated by the field in the iron cylinder. # diff --git a/chapter3/robin_neumann_dirichlet.ipynb b/chapter3/robin_neumann_dirichlet.ipynb index 4d2d8189..ba455964 100644 --- a/chapter3/robin_neumann_dirichlet.ipynb +++ b/chapter3/robin_neumann_dirichlet.ipynb @@ -312,7 +312,7 @@ " BoundaryCondition(\"Dirichlet\", 2, u_ex),\n", " BoundaryCondition(\"Robin\", 3, (r, s)),\n", " BoundaryCondition(\"Neumann\", 4, g),\n", - "]\n" + "]" ] }, { diff --git a/chapter3/robin_neumann_dirichlet.py b/chapter3/robin_neumann_dirichlet.py index 113fcc01..c2d4371a 100644 --- a/chapter3/robin_neumann_dirichlet.py +++ b/chapter3/robin_neumann_dirichlet.py @@ -240,7 +240,6 @@ def type(self): BoundaryCondition("Robin", 3, (r, s)), BoundaryCondition("Neumann", 4, g), ] - # - # We can now loop through the boundary condition and append them to `L(v)` or the list of Dirichlet boundary conditions diff --git a/chapter4/compiler_parameters.ipynb b/chapter4/compiler_parameters.ipynb index 899e2366..39feeb91 100644 --- a/chapter4/compiler_parameters.ipynb +++ b/chapter4/compiler_parameters.ipynb @@ -23,7 +23,9 @@ "cell_type": "code", "execution_count": null, "id": "1", - "metadata": {}, + "metadata": { + "lines_to_end_of_cell_marker": 2 + }, "outputs": [], "source": [ "import pandas as pd\n", @@ -46,7 +48,9 @@ { "cell_type": "markdown", "id": "2", - "metadata": {}, + "metadata": { + "lines_to_next_cell": 2 + }, "source": [ "Next we generate a general function to assemble the mass matrix for a unit cube. Note that we use `dolfinx.fem.form` to compile the variational form.\n", "For codes using `dolfinx.fem.petsc.LinearProblem`, you can supply `jit_options` as a keyword argument." @@ -59,8 +63,6 @@ "metadata": {}, "outputs": [], "source": [ - "\n", - "\n", "def compile_form(space: str, degree: int, jit_options: Dict):\n", " N = 10\n", " mesh = create_unit_cube(MPI.COMM_WORLD, N, N, N)\n", diff --git a/chapter4/compiler_parameters.py b/chapter4/compiler_parameters.py index 53ce95cc..6536b670 100644 --- a/chapter4/compiler_parameters.py +++ b/chapter4/compiler_parameters.py @@ -42,13 +42,13 @@ cache_dir = f"{str(Path.cwd())}/.cache" print(f"Directory to put C files in: {cache_dir}") + + # - # Next we generate a general function to assemble the mass matrix for a unit cube. Note that we use `dolfinx.fem.form` to compile the variational form. # For codes using `dolfinx.fem.petsc.LinearProblem`, you can supply `jit_options` as a keyword argument. -# + - def compile_form(space: str, degree: int, jit_options: Dict): N = 10 @@ -64,8 +64,6 @@ def compile_form(space: str, degree: int, jit_options: Dict): return end - start -# - - # We start by considering the different levels of optimization that the C compiler can use on the optimized code. A list of optimization options and explanations can be found [here](https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html) optimization_options = ["-O1", "-O2", "-O3", "-Ofast"] diff --git a/chapter4/mixed_poisson.ipynb b/chapter4/mixed_poisson.ipynb index 9e5e1c66..0d751e8f 100644 --- a/chapter4/mixed_poisson.ipynb +++ b/chapter4/mixed_poisson.ipynb @@ -84,7 +84,6 @@ "metadata": {}, "outputs": [], "source": [ - "\n", "import numpy as np\n", "\n", "\n", diff --git a/chapter4/mixed_poisson.py b/chapter4/mixed_poisson.py index 0598168a..551afe4b 100644 --- a/chapter4/mixed_poisson.py +++ b/chapter4/mixed_poisson.py @@ -50,7 +50,6 @@ def u_ex(mod, x): # and $\Gamma_N = \{(0, y) \vert y \in [0, 1]\}\cup\{(1, y) \vert y \in [0, 1]\}$. # + - import numpy as np diff --git a/chapter4/newton-solver.ipynb b/chapter4/newton-solver.ipynb index f180907a..f979e173 100644 --- a/chapter4/newton-solver.ipynb +++ b/chapter4/newton-solver.ipynb @@ -388,7 +388,7 @@ "F = q(uh) * ufl.dot(ufl.grad(uh), ufl.grad(v)) * ufl.dx - f * v * ufl.dx\n", "J = ufl.derivative(F, uh)\n", "residual = dolfinx.fem.form(F)\n", - "jacobian = dolfinx.fem.form(J)\n" + "jacobian = dolfinx.fem.form(J)" ] }, { diff --git a/chapter4/newton-solver.py b/chapter4/newton-solver.py index ecdc37f6..b452b755 100644 --- a/chapter4/newton-solver.py +++ b/chapter4/newton-solver.py @@ -219,7 +219,6 @@ def u_exact(x): J = ufl.derivative(F, uh) residual = dolfinx.fem.form(F) jacobian = dolfinx.fem.form(J) - # - # Next, we define the matrix `A`, right hand side vector `L` and the correction function `du` diff --git a/index.ipynb b/index.ipynb index 1cce4752..92713239 100644 --- a/index.ipynb +++ b/index.ipynb @@ -52,7 +52,9 @@ "source": [ "from dolfinx.fem import functionspace # Click `functionspace`\n", "from ufl import div, grad # Click `div` and `grad`\n", - "import numpy as np # Click `numpy`" + "import numpy as np # Click `numpy`\n", + "\n", + "print(functionspace, div, grad, np)" ] }, {