diff --git a/.github/workflows/test_nightly.yml b/.github/workflows/test_nightly.yml index e4f4736e..933c28e4 100644 --- a/.github/workflows/test_nightly.yml +++ b/.github/workflows/test_nightly.yml @@ -67,7 +67,6 @@ jobs: mpirun -n 2 python3 nonlinpoisson_code.py mpirun -n 2 python3 ns_code1.py mpirun -n 2 python3 ns_code2.py - - name: Test chapter 3 working-directory: chapter3 run: | @@ -77,7 +76,6 @@ jobs: mpirun -n 2 python3 robin_neumann_dirichlet.py mpirun -n 2 python3 component_bc.py mpirun -n 2 python3 em.py - - name: Test chapter 4 working-directory: chapter4 run: | diff --git a/Changelog.md b/Changelog.md index 2817a0bb..6dd8d612 100644 --- a/Changelog.md +++ b/Changelog.md @@ -1,5 +1,15 @@ # Changelog +## v0.10.0 + +- Change how one reads in GMSH data with `gmshio`. See [the membrane code](./chapter1/membrane_code.ipynb) for more details. +- `dolfinx.fem.FiniteElement.interpolation_points()` -> `dolfinx.fem.FiniteElement.interpolation_points`. + +## v0.9.0 + +- `scale` in `apply_lifting` has been renamed to `alpha` +- Use `dolfinx.fem.Function.x.petsc_vec` as opposed to `dolfinx.fem.Function.vector` + ## v0.8.0 - Replace all `ufl.FiniteElement` and `ufl.VectorElement` with the appropriate `basix.ufl.element` diff --git a/chapter1/complex_mode.py b/chapter1/complex_mode.py index 034b019a..65981976 100644 --- a/chapter1/complex_mode.py +++ b/chapter1/complex_mode.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.4 +# jupytext_version: 1.16.5 # kernelspec: # display_name: Python 3 (DOLFINx complex) # language: python diff --git a/chapter1/fundamentals_code.py b/chapter1/fundamentals_code.py index 7dcd7629..7670ef10 100644 --- a/chapter1/fundamentals_code.py +++ b/chapter1/fundamentals_code.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.4 +# jupytext_version: 1.16.5 # kernelspec: # display_name: Python 3 (ipykernel) # language: python diff --git a/chapter1/membrane_code.ipynb b/chapter1/membrane_code.ipynb index ccaa69e8..ecff9393 100644 --- a/chapter1/membrane_code.ipynb +++ b/chapter1/membrane_code.ipynb @@ -2,6 +2,7 @@ "cells": [ { "cell_type": "markdown", + "id": "0", "metadata": {}, "source": [ "# Implementation\n", @@ -24,15 +25,26 @@ { "cell_type": "code", "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "import gmsh" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", "metadata": {}, "outputs": [], "source": [ - "import gmsh\n", "gmsh.initialize()" ] }, { "cell_type": "markdown", + "id": "3", "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." @@ -41,6 +53,7 @@ { "cell_type": "code", "execution_count": null, + "id": "4", "metadata": {}, "outputs": [], "source": [ @@ -50,6 +63,7 @@ }, { "cell_type": "markdown", + "id": "5", "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." @@ -58,6 +72,7 @@ { "cell_type": "code", "execution_count": null, + "id": "6", "metadata": {}, "outputs": [], "source": [ @@ -67,6 +82,7 @@ }, { "cell_type": "markdown", + "id": "7", "metadata": {}, "source": [ "Finally, we generate the two-dimensional mesh. We set a uniform mesh size by modifying the GMSH options." @@ -75,6 +91,7 @@ { "cell_type": "code", "execution_count": null, + "id": "8", "metadata": {}, "outputs": [], "source": [ @@ -85,16 +102,31 @@ }, { "cell_type": "markdown", + "id": "9", "metadata": {}, "source": [ "# Interfacing with GMSH in DOLFINx\n", "We will import the GMSH-mesh directly from GMSH into DOLFINx via the `dolfinx.io.gmshio` interface.\n", - "As in this example, we have not specified which process we have created the `gmsh` model on, a model has been created on each mpi process. However, we would like to be able to use a mesh distributed over all processes. Therefore, we take the model generated on rank 0 of `MPI.COMM_WORLD`, and distribute it over all available ranks. We will also get two mesh tags, one for cells marked with physical groups in the mesh and one for facets marked with physical groups. As we did not add any physical groups of dimension `gdim-1`, there will be no entities in the `facet_markers`." + "The `gmshio` module contains two functions\n", + "1. `gmshio.model_to_mesh` which takes in a `gmsh.model` and returns a `dolfinx.io.gmshio.MeshData` object.\n", + "2. `gmshio.read_from_msh` which takes in a path to a `.msh`-file and returns a `dolfinx.io.gmshio.MeshData` object.\n", + "\n", + "The `MeshData` object will contain a `dolfinx.mesh.Mesh`, under the attribute `mesh`.\n", + "This mesh will contain all GMSH Physical Groups of the highest topolgoical 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, that being `vertex_tags`, `edge_tags`, `facet_tags` and `cell_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, and which communicator rank should distribute the mesh.\n", + "The `model_to_mesh` will then load the mesh on the specified rank, and distribute it to the communicator using a mesh partitioner." ] }, { "cell_type": "code", "execution_count": null, + "id": "10", "metadata": {}, "outputs": [], "source": [ @@ -104,11 +136,15 @@ "\n", "gmsh_model_rank = 0\n", "mesh_comm = MPI.COMM_WORLD\n", - "domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)" + "mesh_data = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)\n", + "assert mesh_data.cell_tags is not None\n", + "cell_markers = mesh_data.cell_tags\n", + "domain = mesh_data.mesh" ] }, { "cell_type": "markdown", + "id": "11", "metadata": {}, "source": [ "We define the function space as in the previous tutorial" @@ -117,15 +153,26 @@ { "cell_type": "code", "execution_count": null, + "id": "12", + "metadata": {}, + "outputs": [], + "source": [ + "from dolfinx import fem" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", "metadata": {}, "outputs": [], "source": [ - "from dolfinx import fem\n", "V = fem.functionspace(domain, (\"Lagrange\", 1))" ] }, { "cell_type": "markdown", + "id": "14", "metadata": {}, "source": [ "## Defining a spatially varying load\n", @@ -135,19 +182,30 @@ { "cell_type": "code", "execution_count": null, + "id": "15", "metadata": {}, "outputs": [], "source": [ "import ufl\n", - "from dolfinx import default_scalar_type\n", + "from dolfinx import default_scalar_type" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16", + "metadata": {}, + "outputs": [], + "source": [ "x = ufl.SpatialCoordinate(domain)\n", "beta = fem.Constant(domain, default_scalar_type(12))\n", "R0 = fem.Constant(domain, default_scalar_type(0.3))\n", - "p = 4 * ufl.exp(-beta**2 * (x[0]**2 + (x[1] - R0)**2))" + "p = 4 * ufl.exp(-(beta**2) * (x[0] ** 2 + (x[1] - R0) ** 2))" ] }, { "cell_type": "markdown", + "id": "17", "metadata": {}, "source": [ "## Create a Dirichlet boundary condition using geometrical conditions\n", @@ -157,6 +215,7 @@ { "cell_type": "code", "execution_count": null, + "id": "18", "metadata": {}, "outputs": [], "source": [ @@ -164,7 +223,7 @@ "\n", "\n", "def on_boundary(x):\n", - " return np.isclose(np.sqrt(x[0]**2 + x[1]**2), 1)\n", + " return np.isclose(np.sqrt(x[0] ** 2 + x[1] ** 2), 1)\n", "\n", "\n", "boundary_dofs = fem.locate_dofs_geometrical(V, on_boundary)" @@ -172,6 +231,7 @@ }, { "cell_type": "markdown", + "id": "19", "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." @@ -180,6 +240,7 @@ { "cell_type": "code", "execution_count": null, + "id": "20", "metadata": {}, "outputs": [], "source": [ @@ -188,6 +249,7 @@ }, { "cell_type": "markdown", + "id": "21", "metadata": {}, "source": [ "## Defining the variational problem\n", @@ -197,6 +259,7 @@ { "cell_type": "code", "execution_count": null, + "id": "22", "metadata": {}, "outputs": [], "source": [ @@ -204,12 +267,15 @@ "v = ufl.TestFunction(V)\n", "a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx\n", "L = p * v * ufl.dx\n", - "problem = LinearProblem(a, L, bcs=[bc], petsc_options={\"ksp_type\": \"preonly\", \"pc_type\": \"lu\"})\n", + "problem = LinearProblem(\n", + " a, L, bcs=[bc], petsc_options={\"ksp_type\": \"preonly\", \"pc_type\": \"lu\"}\n", + ")\n", "uh = problem.solve()" ] }, { "cell_type": "markdown", + "id": "23", "metadata": {}, "source": [ "## Interpolation of a `ufl`-expression\n", @@ -220,17 +286,19 @@ { "cell_type": "code", "execution_count": null, + "id": "24", "metadata": {}, "outputs": [], "source": [ "Q = fem.functionspace(domain, (\"Lagrange\", 5))\n", - "expr = fem.Expression(p, Q.element.interpolation_points())\n", + "expr = fem.Expression(p, Q.element.interpolation_points)\n", "pressure = fem.Function(Q)\n", "pressure.interpolate(expr)" ] }, { "cell_type": "markdown", + "id": "25", "metadata": {}, "source": [ "## Plotting the solution over a line\n", @@ -240,11 +308,13 @@ { "cell_type": "code", "execution_count": null, + "id": "26", "metadata": {}, "outputs": [], "source": [ "from dolfinx.plot import vtk_mesh\n", "import pyvista\n", + "\n", "pyvista.start_xvfb()\n", "\n", "# Extract topology from mesh and create pyvista mesh\n", @@ -265,6 +335,7 @@ }, { "cell_type": "markdown", + "id": "27", "metadata": {}, "source": [ "We next plot the load on the domain" @@ -273,6 +344,7 @@ { "cell_type": "code", "execution_count": null, + "id": "28", "metadata": {}, "outputs": [], "source": [ @@ -291,16 +363,18 @@ }, { "cell_type": "markdown", + "id": "29", "metadata": {}, "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. " + "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." ] }, { "cell_type": "code", "execution_count": null, + "id": "30", "metadata": {}, "outputs": [], "source": [ @@ -314,25 +388,37 @@ }, { "cell_type": "markdown", + "id": "31", "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", + "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." ] }, { "cell_type": "code", "execution_count": null, + "id": "32", + "metadata": {}, + "outputs": [], + "source": [ + "from dolfinx import geometry" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "33", "metadata": {}, "outputs": [], "source": [ - "from dolfinx import geometry\n", "bb_tree = geometry.bb_tree(domain, domain.topology.dim)" ] }, { "cell_type": "markdown", + "id": "34", "metadata": {}, "source": [ "Now we can compute which cells the bounding box tree collides with using `dolfinx.geometry.compute_collisions_points`. This function returns a list of cells whose bounding box collide for each input point. As different points might have different number of cells, the data is stored in `dolfinx.cpp.graph.AdjacencyList_int32`, where one can access the cells for the `i`th point by calling `links(i)`.\n", @@ -346,6 +432,7 @@ { "cell_type": "code", "execution_count": null, + "id": "35", "metadata": {}, "outputs": [], "source": [ @@ -363,6 +450,7 @@ }, { "cell_type": "markdown", + "id": "36", "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" @@ -371,6 +459,7 @@ { "cell_type": "code", "execution_count": null, + "id": "37", "metadata": {}, "outputs": [], "source": [ @@ -381,6 +470,7 @@ }, { "cell_type": "markdown", + "id": "38", "metadata": {}, "source": [ "As we now have an array of coordinates and two arrays of function values, we can use `matplotlib` to plot them" @@ -389,12 +479,28 @@ { "cell_type": "code", "execution_count": null, + "id": "39", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "40", "metadata": {}, "outputs": [], "source": [ - "import matplotlib.pyplot as plt\n", "fig = plt.figure()\n", - "plt.plot(points_on_proc[:, 1], 50 * u_values, \"k\", linewidth=2, label=\"Deflection ($\\\\times 50$)\")\n", + "plt.plot(\n", + " points_on_proc[:, 1],\n", + " 50 * u_values,\n", + " \"k\",\n", + " linewidth=2,\n", + " label=\"Deflection ($\\\\times 50$)\",\n", + ")\n", "plt.plot(points_on_proc[:, 1], p_values, \"b--\", linewidth=2, label=\"Load\")\n", "plt.grid(True)\n", "plt.xlabel(\"y\")\n", @@ -405,6 +511,7 @@ }, { "cell_type": "markdown", + "id": "41", "metadata": {}, "source": [ "## Saving functions to file\n", @@ -414,27 +521,34 @@ { "cell_type": "code", "execution_count": null, + "id": "42", "metadata": {}, "outputs": [], "source": [ "import dolfinx.io\n", - "from pathlib import Path\n", - "pressure.name = \"Load\"\n", - "uh.name = \"Deflection\"\n", - "results_folder = Path(\"results\")\n", - "results_folder.mkdir(exist_ok=True, parents=True)\n", - "with dolfinx.io.VTXWriter(MPI.COMM_WORLD, results_folder / \"membrane_pressure.bp\", [pressure], engine=\"BP4\") as vtx:\n", - " vtx.write(0.0)\n", - "with dolfinx.io.VTXWriter(MPI.COMM_WORLD, results_folder / \"membrane_deflection.bp\", [uh], engine=\"BP4\") as vtx:\n", - " vtx.write(0.0)" + "from pathlib import Path" ] }, { "cell_type": "code", "execution_count": null, + "id": "43", "metadata": {}, "outputs": [], - "source": [] + "source": [ + "pressure.name = \"Load\"\n", + "uh.name = \"Deflection\"\n", + "results_folder = Path(\"results\")\n", + "results_folder.mkdir(exist_ok=True, parents=True)\n", + "with dolfinx.io.VTXWriter(\n", + " MPI.COMM_WORLD, results_folder / \"membrane_pressure.bp\", [pressure], engine=\"BP4\"\n", + ") as vtx:\n", + " vtx.write(0.0)\n", + "with dolfinx.io.VTXWriter(\n", + " MPI.COMM_WORLD, results_folder / \"membrane_deflection.bp\", [uh], engine=\"BP4\"\n", + ") as vtx:\n", + " vtx.write(0.0)" + ] } ], "metadata": { @@ -445,20 +559,8 @@ "display_name": "Python 3 (ipykernel)", "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.10.12" } }, "nbformat": 4, - "nbformat_minor": 4 + "nbformat_minor": 5 } diff --git a/chapter1/membrane_code.py b/chapter1/membrane_code.py index b87f3cec..57b03437 100644 --- a/chapter1/membrane_code.py +++ b/chapter1/membrane_code.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.4 +# jupytext_version: 1.16.5 # kernelspec: # display_name: Python 3 (ipykernel) # language: python @@ -30,6 +30,7 @@ # 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 + 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. @@ -50,7 +51,20 @@ # # Interfacing with GMSH in DOLFINx # We will import the GMSH-mesh directly from GMSH into DOLFINx via the `dolfinx.io.gmshio` interface. -# As in this example, we have not specified which process we have created the `gmsh` model on, a model has been created on each mpi process. However, we would like to be able to use a mesh distributed over all processes. Therefore, we take the model generated on rank 0 of `MPI.COMM_WORLD`, and distribute it over all available ranks. We will also get two mesh tags, one for cells marked with physical groups in the mesh and one for facets marked with physical groups. As we did not add any physical groups of dimension `gdim-1`, there will be no entities in the `facet_markers`. +# The `gmshio` module contains two functions +# 1. `gmshio.model_to_mesh` which takes in a `gmsh.model` and returns a `dolfinx.io.gmshio.MeshData` object. +# 2. `gmshio.read_from_msh` which takes in a path to a `.msh`-file and returns a `dolfinx.io.gmshio.MeshData` object. +# +# The `MeshData` object will contain a `dolfinx.mesh.Mesh`, under the attribute `mesh`. +# This mesh will contain all GMSH Physical Groups of the highest topolgoical 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`. +# 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, and distribute it to the communicator using a mesh partitioner. # + from dolfinx.io import gmshio @@ -59,12 +73,16 @@ gmsh_model_rank = 0 mesh_comm = MPI.COMM_WORLD -domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim) +mesh_data = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim) +assert mesh_data.cell_tags is not None +cell_markers = mesh_data.cell_tags +domain = mesh_data.mesh # - # We define the function space as in the previous tutorial from dolfinx import fem + V = fem.functionspace(domain, ("Lagrange", 1)) # ## Defining a spatially varying load @@ -72,10 +90,11 @@ import ufl from dolfinx import default_scalar_type + x = ufl.SpatialCoordinate(domain) 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)) +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$. @@ -85,7 +104,7 @@ def on_boundary(x): - return np.isclose(np.sqrt(x[0]**2 + x[1]**2), 1) + return np.isclose(np.sqrt(x[0] ** 2 + x[1] ** 2), 1) boundary_dofs = fem.locate_dofs_geometrical(V, on_boundary) @@ -102,7 +121,9 @@ def on_boundary(x): v = ufl.TestFunction(V) a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx L = p * v * ufl.dx -problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"}) +problem = LinearProblem( + a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"} +) uh = problem.solve() # ## Interpolation of a `ufl`-expression @@ -110,7 +131,7 @@ def on_boundary(x): # 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)) -expr = fem.Expression(p, Q.element.interpolation_points()) +expr = fem.Expression(p, Q.element.interpolation_points) pressure = fem.Function(Q) pressure.interpolate(expr) @@ -120,6 +141,7 @@ def on_boundary(x): # + from dolfinx.plot import vtk_mesh import pyvista + pyvista.start_xvfb() # Extract topology from mesh and create pyvista mesh @@ -153,8 +175,8 @@ def on_boundary(x): load_plotter.screenshot("load.png") # ## 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. +# 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. tol = 0.001 # Avoid hitting the outside of the domain y = np.linspace(-1 + tol, 1 - tol, 101) @@ -164,10 +186,11 @@ def on_boundary(x): 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. +# 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. from dolfinx import geometry + bb_tree = geometry.bb_tree(domain, domain.topology.dim) # Now we can compute which cells the bounding box tree collides with using `dolfinx.geometry.compute_collisions_points`. This function returns a list of cells whose bounding box collide for each input point. As different points might have different number of cells, the data is stored in `dolfinx.cpp.graph.AdjacencyList_int32`, where one can access the cells for the `i`th point by calling `links(i)`. @@ -198,8 +221,15 @@ def on_boundary(x): # As we now have an array of coordinates and two arrays of function values, we can use `matplotlib` to plot them import matplotlib.pyplot as plt + fig = plt.figure() -plt.plot(points_on_proc[:, 1], 50 * u_values, "k", linewidth=2, label="Deflection ($\\times 50$)") +plt.plot( + points_on_proc[:, 1], + 50 * u_values, + "k", + linewidth=2, + label="Deflection ($\\times 50$)", +) plt.plot(points_on_proc[:, 1], p_values, "b--", linewidth=2, label="Load") plt.grid(True) plt.xlabel("y") @@ -212,13 +242,16 @@ def on_boundary(x): import dolfinx.io from pathlib import Path + pressure.name = "Load" uh.name = "Deflection" results_folder = Path("results") results_folder.mkdir(exist_ok=True, parents=True) -with dolfinx.io.VTXWriter(MPI.COMM_WORLD, results_folder / "membrane_pressure.bp", [pressure], engine="BP4") as vtx: +with dolfinx.io.VTXWriter( + MPI.COMM_WORLD, results_folder / "membrane_pressure.bp", [pressure], engine="BP4" +) as vtx: vtx.write(0.0) -with dolfinx.io.VTXWriter(MPI.COMM_WORLD, results_folder / "membrane_deflection.bp", [uh], engine="BP4") as vtx: +with dolfinx.io.VTXWriter( + MPI.COMM_WORLD, results_folder / "membrane_deflection.bp", [uh], engine="BP4" +) as vtx: vtx.write(0.0) - - diff --git a/chapter1/nitsche.ipynb b/chapter1/nitsche.ipynb index e093878c..bc0fa6ec 100644 --- a/chapter1/nitsche.ipynb +++ b/chapter1/nitsche.ipynb @@ -51,7 +51,7 @@ "uD = fem.Function(V)\n", "x = SpatialCoordinate(domain)\n", "u_ex = 1 + x[0]**2 + 2 * x[1]**2\n", - "uD.interpolate(fem.Expression(u_ex, V.element.interpolation_points()))\n", + "uD.interpolate(fem.Expression(u_ex, V.element.interpolation_points))\n", "f = -div(grad(u_ex))" ] }, diff --git a/chapter1/nitsche.py b/chapter1/nitsche.py index a71f688e..03818d3d 100644 --- a/chapter1/nitsche.py +++ b/chapter1/nitsche.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.4 +# jupytext_version: 1.16.5 # kernelspec: # display_name: Python 3 (ipykernel) # language: python @@ -39,7 +39,7 @@ uD = fem.Function(V) x = SpatialCoordinate(domain) u_ex = 1 + x[0]**2 + 2 * x[1]**2 -uD.interpolate(fem.Expression(u_ex, V.element.interpolation_points())) +uD.interpolate(fem.Expression(u_ex, V.element.interpolation_points)) f = -div(grad(u_ex)) # As opposed to the first tutorial, we now have to have another look at the variational form. diff --git a/chapter2/diffusion_code.py b/chapter2/diffusion_code.py index ea42b3c4..ff693ddf 100644 --- a/chapter2/diffusion_code.py +++ b/chapter2/diffusion_code.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.4 +# jupytext_version: 1.16.5 # kernelspec: # display_name: Python 3 (ipykernel) # language: python diff --git a/chapter2/heat_code.py b/chapter2/heat_code.py index 556e4d43..9db0d9bc 100644 --- a/chapter2/heat_code.py +++ b/chapter2/heat_code.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.4 +# jupytext_version: 1.16.5 # kernelspec: # display_name: Python 3 (ipykernel) # language: python diff --git a/chapter2/hyperelasticity.ipynb b/chapter2/hyperelasticity.ipynb index 1da7becd..3906a2d4 100644 --- a/chapter2/hyperelasticity.ipynb +++ b/chapter2/hyperelasticity.ipynb @@ -343,7 +343,7 @@ "# Compute magnitude of displacement to visualize in GIF\n", "Vs = fem.functionspace(domain, (\"Lagrange\", 2))\n", "magnitude = fem.Function(Vs)\n", - "us = fem.Expression(ufl.sqrt(sum([u[i]**2 for i in range(len(u))])), Vs.element.interpolation_points())\n", + "us = fem.Expression(ufl.sqrt(sum([u[i]**2 for i in range(len(u))])), Vs.element.interpolation_points)\n", "magnitude.interpolate(us)\n", "warped[\"mag\"] = magnitude.x.array" ] diff --git a/chapter2/hyperelasticity.py b/chapter2/hyperelasticity.py index 0916f30a..0a9d4fd9 100644 --- a/chapter2/hyperelasticity.py +++ b/chapter2/hyperelasticity.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.4 +# jupytext_version: 1.16.5 # kernelspec: # display_name: Python 3 (ipykernel) # language: python @@ -173,7 +173,7 @@ def right(x): # Compute magnitude of displacement to visualize in GIF Vs = fem.functionspace(domain, ("Lagrange", 2)) magnitude = fem.Function(Vs) -us = fem.Expression(ufl.sqrt(sum([u[i]**2 for i in range(len(u))])), Vs.element.interpolation_points()) +us = fem.Expression(ufl.sqrt(sum([u[i]**2 for i in range(len(u))])), Vs.element.interpolation_points) magnitude.interpolate(us) warped["mag"] = magnitude.x.array # - diff --git a/chapter2/linearelasticity_code.ipynb b/chapter2/linearelasticity_code.ipynb index 99e27b3f..2f3e8c5e 100644 --- a/chapter2/linearelasticity_code.ipynb +++ b/chapter2/linearelasticity_code.ipynb @@ -270,7 +270,7 @@ "outputs": [], "source": [ "V_von_mises = fem.functionspace(domain, (\"DG\", 0))\n", - "stress_expr = fem.Expression(von_Mises, V_von_mises.element.interpolation_points())\n", + "stress_expr = fem.Expression(von_Mises, V_von_mises.element.interpolation_points)\n", "stresses = fem.Function(V_von_mises)\n", "stresses.interpolate(stress_expr)" ] diff --git a/chapter2/linearelasticity_code.py b/chapter2/linearelasticity_code.py index 038e0cdc..816d171d 100644 --- a/chapter2/linearelasticity_code.py +++ b/chapter2/linearelasticity_code.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.4 +# jupytext_version: 1.16.5 # kernelspec: # display_name: Python 3 (ipykernel) # language: python @@ -151,7 +151,7 @@ def sigma(u): # The `von_Mises` variable is now an expression that must be projected into an appropriate function space so that we can visualize it. As `uh` is a linear combination of first order piecewise continuous functions, the von Mises stresses will be a cell-wise constant function. V_von_mises = fem.functionspace(domain, ("DG", 0)) -stress_expr = fem.Expression(von_Mises, V_von_mises.element.interpolation_points()) +stress_expr = fem.Expression(von_Mises, V_von_mises.element.interpolation_points) stresses = fem.Function(V_von_mises) stresses.interpolate(stress_expr) diff --git a/chapter2/nonlinpoisson_code.py b/chapter2/nonlinpoisson_code.py index 09dcb091..f713c623 100644 --- a/chapter2/nonlinpoisson_code.py +++ b/chapter2/nonlinpoisson_code.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.4 +# jupytext_version: 1.16.5 # kernelspec: # display_name: Python 3 (ipykernel) # language: python diff --git a/chapter2/ns_code1.py b/chapter2/ns_code1.py index 7cc0dcee..bb235bd7 100644 --- a/chapter2/ns_code1.py +++ b/chapter2/ns_code1.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.4 +# jupytext_version: 1.16.5 # kernelspec: # display_name: Python 3 (ipykernel) # language: python diff --git a/chapter2/ns_code2.ipynb b/chapter2/ns_code2.ipynb index 237dca40..3ea55b76 100644 --- a/chapter2/ns_code2.ipynb +++ b/chapter2/ns_code2.ipynb @@ -2,6 +2,7 @@ "cells": [ { "cell_type": "markdown", + "id": "0", "metadata": {}, "source": [ "# Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark)\n", @@ -39,9 +40,8 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "tags": [] - }, + "id": "1", + "metadata": {}, "outputs": [], "source": [ "import gmsh\n", @@ -55,17 +55,41 @@ "\n", "from basix.ufl import element\n", "\n", - "from dolfinx.cpp.mesh import to_type, cell_entity_type\n", - "from dolfinx.fem import (Constant, Function, functionspace,\n", - " assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc)\n", - "from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,\n", - " create_vector, create_matrix, set_bc)\n", - "from dolfinx.graph import adjacencylist\n", + "from dolfinx.fem import (\n", + " Constant,\n", + " Function,\n", + " functionspace,\n", + " assemble_scalar,\n", + " dirichletbc,\n", + " form,\n", + " locate_dofs_topological,\n", + " set_bc,\n", + ")\n", + "from dolfinx.fem.petsc import (\n", + " apply_lifting,\n", + " assemble_matrix,\n", + " assemble_vector,\n", + " create_vector,\n", + " create_matrix,\n", + " set_bc,\n", + ")\n", "from dolfinx.geometry import bb_tree, compute_collisions_points, compute_colliding_cells\n", - "from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio)\n", - "from dolfinx.mesh import create_mesh, meshtags_from_entities\n", - "from ufl import (FacetNormal, Identity, Measure, TestFunction, TrialFunction,\n", - " as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym, system)\n", + "from dolfinx.io import VTXWriter, gmshio\n", + "from ufl import (\n", + " FacetNormal,\n", + " Measure,\n", + " TestFunction,\n", + " TrialFunction,\n", + " as_vector,\n", + " div,\n", + " dot,\n", + " dx,\n", + " inner,\n", + " lhs,\n", + " grad,\n", + " nabla_grad,\n", + " rhs,\n", + ")\n", "\n", "gmsh.initialize()\n", "\n", @@ -83,6 +107,7 @@ }, { "cell_type": "markdown", + "id": "2", "metadata": {}, "source": [ "The next step is to subtract the obstacle from the channel, such that we do not mesh the interior of the circle.\n" @@ -91,6 +116,7 @@ { "cell_type": "code", "execution_count": null, + "id": "3", "metadata": {}, "outputs": [], "source": [ @@ -101,6 +127,7 @@ }, { "cell_type": "markdown", + "id": "4", "metadata": {}, "source": [ "To get GMSH to mesh the fluid, we add a physical volume marker\n" @@ -109,19 +136,21 @@ { "cell_type": "code", "execution_count": null, + "id": "5", "metadata": {}, "outputs": [], "source": [ "fluid_marker = 1\n", "if mesh_comm.rank == model_rank:\n", " volumes = gmsh.model.getEntities(dim=gdim)\n", - " assert (len(volumes) == 1)\n", + " assert len(volumes) == 1\n", " gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)\n", " gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, \"Fluid\")" ] }, { "cell_type": "markdown", + "id": "6", "metadata": {}, "source": [ "To tag the different surfaces of the mesh, we tag the inflow (left hand side) with marker 2, the outflow (right hand side) with marker 3 and the fluid walls with 4 and obstacle with 5. We will do this by computing the center of mass for each geometrical entity.\n" @@ -130,6 +159,7 @@ { "cell_type": "code", "execution_count": null, + "id": "7", "metadata": {}, "outputs": [], "source": [ @@ -143,7 +173,9 @@ " inflow.append(boundary[1])\n", " elif np.allclose(center_of_mass, [L, H / 2, 0]):\n", " outflow.append(boundary[1])\n", - " elif np.allclose(center_of_mass, [L / 2, H, 0]) or np.allclose(center_of_mass, [L / 2, 0, 0]):\n", + " elif np.allclose(center_of_mass, [L / 2, H, 0]) or np.allclose(\n", + " center_of_mass, [L / 2, 0, 0]\n", + " ):\n", " walls.append(boundary[1])\n", " else:\n", " obstacle.append(boundary[1])\n", @@ -159,6 +191,7 @@ }, { "cell_type": "markdown", + "id": "8", "metadata": {}, "source": [ "In our previous meshes, we have used uniform mesh sizes. In this example, we will have variable mesh sizes to resolve the flow solution in the area of interest; close to the circular obstacle. To do this, we use GMSH Fields.\n" @@ -167,6 +200,7 @@ { "cell_type": "code", "execution_count": null, + "id": "9", "metadata": {}, "outputs": [], "source": [ @@ -194,6 +228,7 @@ }, { "cell_type": "markdown", + "id": "10", "metadata": {}, "source": [ "## Generating the mesh\n", @@ -204,6 +239,7 @@ { "cell_type": "code", "execution_count": null, + "id": "11", "metadata": {}, "outputs": [], "source": [ @@ -219,26 +255,33 @@ }, { "cell_type": "markdown", + "id": "12", "metadata": {}, "source": [ "## Loading mesh and boundary markers\n", "\n", "As we have generated the mesh, we now need to load the mesh and corresponding facet markers into DOLFINx.\n", - "To load the mesh, we follow the same structure as in [Deflection of a membrane](./../chapter1/membrane_code.ipynb), with the difference being that we will load in facet markers as well. To learn more about the specifics of the function below, see [A GMSH tutorial for DOLFINx](https://jsdokken.com/src/tutorial_gmsh.html).\n" + "To load the mesh, we follow the same structure as in [Deflection of a membrane](./../chapter1/membrane_code.ipynb), with the difference being that we will load in facet markers as well.\n", + "To learn more about the specifics of the function below, see [A GMSH tutorial for DOLFINx](https://jsdokken.com/src/tutorial_gmsh.html).\n" ] }, { "cell_type": "code", "execution_count": null, + "id": "13", "metadata": {}, "outputs": [], "source": [ - "mesh, _, ft = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)\n", + "mesh_data = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)\n", + "mesh = mesh_data.mesh\n", + "assert mesh_data.facet_tags is not None\n", + "ft = mesh_data.facet_tags\n", "ft.name = \"Facet markers\"" ] }, { "cell_type": "markdown", + "id": "14", "metadata": {}, "source": [ "## Physical and discretization parameters\n", @@ -249,20 +292,22 @@ { "cell_type": "code", "execution_count": null, + "id": "15", "metadata": {}, "outputs": [], "source": [ "t = 0\n", - "T = 8 # Final time\n", - "dt = 1 / 1600 # Time step size\n", + "T = 8 # Final time\n", + "dt = 1 / 1600 # Time step size\n", "num_steps = int(T / dt)\n", "k = Constant(mesh, PETSc.ScalarType(dt))\n", "mu = Constant(mesh, PETSc.ScalarType(0.001)) # Dynamic viscosity\n", - "rho = Constant(mesh, PETSc.ScalarType(1)) # Density" + "rho = Constant(mesh, PETSc.ScalarType(1)) # Density" ] }, { "cell_type": "markdown", + "id": "16", "metadata": {}, "source": [ "```{admonition} Reduced end-time of problem\n", @@ -277,12 +322,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "tags": [] - }, + "id": "17", + "metadata": {}, "outputs": [], "source": [ - "v_cg2 = element(\"Lagrange\", mesh.topology.cell_name(), 2, shape=(mesh.geometry.dim, ))\n", + "v_cg2 = element(\"Lagrange\", mesh.topology.cell_name(), 2, shape=(mesh.geometry.dim,))\n", "s_cg1 = element(\"Lagrange\", mesh.topology.cell_name(), 1)\n", "V = functionspace(mesh, v_cg2)\n", "Q = functionspace(mesh, s_cg1)\n", @@ -292,13 +336,15 @@ "# Define boundary conditions\n", "\n", "\n", - "class InletVelocity():\n", + "class InletVelocity:\n", " def __init__(self, t):\n", " self.t = t\n", "\n", " def __call__(self, x):\n", " values = np.zeros((gdim, x.shape[1]), dtype=PETSc.ScalarType)\n", - " values[0] = 4 * 1.5 * np.sin(self.t * np.pi / 8) * x[1] * (0.41 - x[1]) / (0.41**2)\n", + " values[0] = (\n", + " 4 * 1.5 * np.sin(self.t * np.pi / 8) * x[1] * (0.41 - x[1]) / (0.41**2)\n", + " )\n", " return values\n", "\n", "\n", @@ -306,20 +352,29 @@ "u_inlet = Function(V)\n", "inlet_velocity = InletVelocity(t)\n", "u_inlet.interpolate(inlet_velocity)\n", - "bcu_inflow = dirichletbc(u_inlet, locate_dofs_topological(V, fdim, ft.find(inlet_marker)))\n", + "bcu_inflow = dirichletbc(\n", + " u_inlet, locate_dofs_topological(V, fdim, ft.find(inlet_marker))\n", + ")\n", "# Walls\n", "u_nonslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)\n", - "bcu_walls = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(wall_marker)), V)\n", + "bcu_walls = dirichletbc(\n", + " u_nonslip, locate_dofs_topological(V, fdim, ft.find(wall_marker)), V\n", + ")\n", "# Obstacle\n", - "bcu_obstacle = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(obstacle_marker)), V)\n", + "bcu_obstacle = dirichletbc(\n", + " u_nonslip, locate_dofs_topological(V, fdim, ft.find(obstacle_marker)), V\n", + ")\n", "bcu = [bcu_inflow, bcu_obstacle, bcu_walls]\n", "# Outlet\n", - "bcp_outlet = dirichletbc(PETSc.ScalarType(0), locate_dofs_topological(Q, fdim, ft.find(outlet_marker)), Q)\n", + "bcp_outlet = dirichletbc(\n", + " PETSc.ScalarType(0), locate_dofs_topological(Q, fdim, ft.find(outlet_marker)), Q\n", + ")\n", "bcp = [bcp_outlet]" ] }, { "cell_type": "markdown", + "id": "18", "metadata": {}, "source": [ "## Variational form\n", @@ -366,6 +421,7 @@ { "cell_type": "code", "execution_count": null, + "id": "19", "metadata": {}, "outputs": [], "source": [ @@ -385,6 +441,7 @@ }, { "cell_type": "markdown", + "id": "20", "metadata": {}, "source": [ "Next, we define the variational formulation for the first step, where we have integrated the diffusion term, as well as the pressure term by parts.\n" @@ -393,6 +450,7 @@ { "cell_type": "code", "execution_count": null, + "id": "21", "metadata": {}, "outputs": [], "source": [ @@ -409,6 +467,7 @@ }, { "cell_type": "markdown", + "id": "22", "metadata": {}, "source": [ "Next we define the second step\n" @@ -417,6 +476,7 @@ { "cell_type": "code", "execution_count": null, + "id": "23", "metadata": {}, "outputs": [], "source": [ @@ -429,6 +489,7 @@ }, { "cell_type": "markdown", + "id": "24", "metadata": {}, "source": [ "We finally create the last step\n" @@ -437,6 +498,7 @@ { "cell_type": "code", "execution_count": null, + "id": "25", "metadata": {}, "outputs": [], "source": [ @@ -449,6 +511,7 @@ }, { "cell_type": "markdown", + "id": "26", "metadata": {}, "source": [ "As in the previous tutorials, we use PETSc as a linear algebra backend.\n" @@ -457,6 +520,7 @@ { "cell_type": "code", "execution_count": null, + "id": "27", "metadata": {}, "outputs": [], "source": [ @@ -485,6 +549,7 @@ }, { "cell_type": "markdown", + "id": "28", "metadata": {}, "source": [ "## Verification of the implementation compute known physical quantities\n", @@ -505,6 +570,7 @@ { "cell_type": "code", "execution_count": null, + "id": "29", "metadata": {}, "outputs": [], "source": [ @@ -522,6 +588,7 @@ }, { "cell_type": "markdown", + "id": "30", "metadata": {}, "source": [ "We will also evaluate the pressure at two points, one in front of the obstacle, $(0.15, 0.2)$, and one behind the obstacle, $(0.25, 0.2)$. To do this, we have to find which cell contains each of the points, so that we can create a linear combination of the local basis functions and coefficients.\n" @@ -530,6 +597,7 @@ { "cell_type": "code", "execution_count": null, + "id": "31", "metadata": {}, "outputs": [], "source": [ @@ -545,6 +613,7 @@ }, { "cell_type": "markdown", + "id": "32", "metadata": {}, "source": [ "## Solving the time-dependent problem\n", @@ -560,16 +629,24 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "tags": [] - }, + "id": "33", + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "34", + "metadata": {}, "outputs": [], "source": [ - "from pathlib import Path\n", "folder = Path(\"results\")\n", "folder.mkdir(exist_ok=True, parents=True)\n", - "vtx_u = VTXWriter(mesh.comm, \"dfg2D-3-u.bp\", [u_], engine=\"BP4\")\n", - "vtx_p = VTXWriter(mesh.comm, \"dfg2D-3-p.bp\", [p_], engine=\"BP4\")\n", + "vtx_u = VTXWriter(mesh.comm, folder / \"dfg2D-3-u.bp\", [u_], engine=\"BP4\")\n", + "vtx_p = VTXWriter(mesh.comm, folder / \"dfg2D-3-p.bp\", [p_], engine=\"BP4\")\n", "vtx_u.write(t)\n", "vtx_p.write(t)\n", "progress = tqdm.autonotebook.tqdm(desc=\"Solving PDE\", total=num_steps)\n", @@ -620,7 +697,11 @@ " vtx_p.write(t)\n", "\n", " # Update variable with solution form this time step\n", - " with u_.x.petsc_vec.localForm() as loc_, u_n.x.petsc_vec.localForm() as loc_n, u_n1.x.petsc_vec.localForm() as loc_n1:\n", + " with (\n", + " u_.x.petsc_vec.localForm() as loc_,\n", + " u_n.x.petsc_vec.localForm() as loc_n,\n", + " u_n1.x.petsc_vec.localForm() as loc_n1,\n", + " ):\n", " loc_n.copy(loc_n1)\n", " loc_.copy(loc_n)\n", "\n", @@ -658,6 +739,7 @@ }, { "cell_type": "markdown", + "id": "35", "metadata": {}, "source": [ "## Verification using data from FEATFLOW\n", @@ -668,9 +750,8 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "tags": [] - }, + "id": "36", + "metadata": {}, "outputs": [], "source": [ "if mesh.comm.rank == 0:\n", @@ -682,40 +763,68 @@ " turek = np.loadtxt(\"bdforces_lv4\")\n", " turek_p = np.loadtxt(\"pointvalues_lv4\")\n", " fig = plt.figure(figsize=(25, 8))\n", - " l1 = plt.plot(t_u, C_D, label=r\"FEniCSx ({0:d} dofs)\".format(num_velocity_dofs + num_pressure_dofs), linewidth=2)\n", - " l2 = plt.plot(turek[1:, 1], turek[1:, 3], marker=\"x\", markevery=50,\n", - " linestyle=\"\", markersize=4, label=\"FEATFLOW (42016 dofs)\")\n", + " l1 = plt.plot(\n", + " t_u,\n", + " C_D,\n", + " label=r\"FEniCSx ({0:d} dofs)\".format(num_velocity_dofs + num_pressure_dofs),\n", + " linewidth=2,\n", + " )\n", + " l2 = plt.plot(\n", + " turek[1:, 1],\n", + " turek[1:, 3],\n", + " marker=\"x\",\n", + " markevery=50,\n", + " linestyle=\"\",\n", + " markersize=4,\n", + " label=\"FEATFLOW (42016 dofs)\",\n", + " )\n", " plt.title(\"Drag coefficient\")\n", " plt.grid()\n", " plt.legend()\n", " plt.savefig(\"figures/drag_comparison.png\")\n", "\n", " fig = plt.figure(figsize=(25, 8))\n", - " l1 = plt.plot(t_u, C_L, label=r\"FEniCSx ({0:d} dofs)\".format(\n", - " num_velocity_dofs + num_pressure_dofs), linewidth=2)\n", - " l2 = plt.plot(turek[1:, 1], turek[1:, 4], marker=\"x\", markevery=50,\n", - " linestyle=\"\", markersize=4, label=\"FEATFLOW (42016 dofs)\")\n", + " l1 = plt.plot(\n", + " t_u,\n", + " C_L,\n", + " label=r\"FEniCSx ({0:d} dofs)\".format(num_velocity_dofs + num_pressure_dofs),\n", + " linewidth=2,\n", + " )\n", + " l2 = plt.plot(\n", + " turek[1:, 1],\n", + " turek[1:, 4],\n", + " marker=\"x\",\n", + " markevery=50,\n", + " linestyle=\"\",\n", + " markersize=4,\n", + " label=\"FEATFLOW (42016 dofs)\",\n", + " )\n", " plt.title(\"Lift coefficient\")\n", " plt.grid()\n", " plt.legend()\n", " plt.savefig(\"figures/lift_comparison.png\")\n", "\n", " fig = plt.figure(figsize=(25, 8))\n", - " l1 = plt.plot(t_p, p_diff, label=r\"FEniCSx ({0:d} dofs)\".format(num_velocity_dofs + num_pressure_dofs), linewidth=2)\n", - " l2 = plt.plot(turek[1:, 1], turek_p[1:, 6] - turek_p[1:, -1], marker=\"x\", markevery=50,\n", - " linestyle=\"\", markersize=4, label=\"FEATFLOW (42016 dofs)\")\n", + " l1 = plt.plot(\n", + " t_p,\n", + " p_diff,\n", + " label=r\"FEniCSx ({0:d} dofs)\".format(num_velocity_dofs + num_pressure_dofs),\n", + " linewidth=2,\n", + " )\n", + " l2 = plt.plot(\n", + " turek[1:, 1],\n", + " turek_p[1:, 6] - turek_p[1:, -1],\n", + " marker=\"x\",\n", + " markevery=50,\n", + " linestyle=\"\",\n", + " markersize=4,\n", + " label=\"FEATFLOW (42016 dofs)\",\n", + " )\n", " plt.title(\"Pressure difference\")\n", " plt.grid()\n", " plt.legend()\n", " plt.savefig(\"figures/pressure_comparison.png\")" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { @@ -726,20 +835,8 @@ "display_name": "Python 3 (ipykernel)", "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.10.12" } }, "nbformat": 4, - "nbformat_minor": 4 + "nbformat_minor": 5 } diff --git a/chapter2/ns_code2.py b/chapter2/ns_code2.py index 1c0be052..08cddb4a 100644 --- a/chapter2/ns_code2.py +++ b/chapter2/ns_code2.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.4 +# jupytext_version: 1.16.5 # kernelspec: # display_name: Python 3 (ipykernel) # language: python @@ -57,17 +57,41 @@ from basix.ufl import element -from dolfinx.cpp.mesh import to_type, cell_entity_type -from dolfinx.fem import (Constant, Function, functionspace, - assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc) -from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector, - create_vector, create_matrix, set_bc) -from dolfinx.graph import adjacencylist +from dolfinx.fem import ( + Constant, + Function, + functionspace, + assemble_scalar, + dirichletbc, + form, + locate_dofs_topological, + set_bc, +) +from dolfinx.fem.petsc import ( + apply_lifting, + assemble_matrix, + assemble_vector, + create_vector, + create_matrix, + set_bc, +) from dolfinx.geometry import bb_tree, compute_collisions_points, compute_colliding_cells -from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio) -from dolfinx.mesh import create_mesh, meshtags_from_entities -from ufl import (FacetNormal, Identity, Measure, TestFunction, TrialFunction, - as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym, system) +from dolfinx.io import VTXWriter, gmshio +from ufl import ( + FacetNormal, + Measure, + TestFunction, + TrialFunction, + as_vector, + div, + dot, + dx, + inner, + lhs, + grad, + nabla_grad, + rhs, +) gmsh.initialize() @@ -96,7 +120,7 @@ fluid_marker = 1 if mesh_comm.rank == model_rank: volumes = gmsh.model.getEntities(dim=gdim) - assert (len(volumes) == 1) + assert len(volumes) == 1 gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker) gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid") @@ -113,7 +137,9 @@ inflow.append(boundary[1]) elif np.allclose(center_of_mass, [L, H / 2, 0]): outflow.append(boundary[1]) - elif np.allclose(center_of_mass, [L / 2, H, 0]) or np.allclose(center_of_mass, [L / 2, 0, 0]): + elif np.allclose(center_of_mass, [L / 2, H, 0]) or np.allclose( + center_of_mass, [L / 2, 0, 0] + ): walls.append(boundary[1]) else: obstacle.append(boundary[1]) @@ -167,10 +193,14 @@ # ## Loading mesh and boundary markers # # As we have generated the mesh, we now need to load the mesh and corresponding facet markers into DOLFINx. -# To load the mesh, we follow the same structure as in [Deflection of a membrane](./../chapter1/membrane_code.ipynb), with the difference being that we will load in facet markers as well. To learn more about the specifics of the function below, see [A GMSH tutorial for DOLFINx](https://jsdokken.com/src/tutorial_gmsh.html). +# To load the mesh, we follow the same structure as in [Deflection of a membrane](./../chapter1/membrane_code.ipynb), with the difference being that we will load in facet markers as well. +# To learn more about the specifics of the function below, see [A GMSH tutorial for DOLFINx](https://jsdokken.com/src/tutorial_gmsh.html). # -mesh, _, ft = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim) +mesh_data = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim) +mesh = mesh_data.mesh +assert mesh_data.facet_tags is not None +ft = mesh_data.facet_tags ft.name = "Facet markers" # ## Physical and discretization parameters @@ -179,12 +209,12 @@ # t = 0 -T = 8 # Final time -dt = 1 / 1600 # Time step size +T = 8 # Final time +dt = 1 / 1600 # Time step size num_steps = int(T / dt) k = Constant(mesh, PETSc.ScalarType(dt)) mu = Constant(mesh, PETSc.ScalarType(0.001)) # Dynamic viscosity -rho = Constant(mesh, PETSc.ScalarType(1)) # Density +rho = Constant(mesh, PETSc.ScalarType(1)) # Density # ```{admonition} Reduced end-time of problem # In the current demo, we have reduced the run time to one second to make it easier to illustrate the concepts of the benchmark. By increasing the end-time `T` to 8, the runtime in a notebook is approximately 25 minutes. If you convert the notebook to a python file and use `mpirun`, you can reduce the runtime of the problem. @@ -196,7 +226,7 @@ # # + -v_cg2 = element("Lagrange", mesh.topology.cell_name(), 2, shape=(mesh.geometry.dim, )) +v_cg2 = element("Lagrange", mesh.topology.cell_name(), 2, shape=(mesh.geometry.dim,)) s_cg1 = element("Lagrange", mesh.topology.cell_name(), 1) V = functionspace(mesh, v_cg2) Q = functionspace(mesh, s_cg1) @@ -206,13 +236,15 @@ # Define boundary conditions -class InletVelocity(): +class InletVelocity: def __init__(self, t): self.t = t def __call__(self, x): values = np.zeros((gdim, x.shape[1]), dtype=PETSc.ScalarType) - values[0] = 4 * 1.5 * np.sin(self.t * np.pi / 8) * x[1] * (0.41 - x[1]) / (0.41**2) + values[0] = ( + 4 * 1.5 * np.sin(self.t * np.pi / 8) * x[1] * (0.41 - x[1]) / (0.41**2) + ) return values @@ -220,15 +252,23 @@ def __call__(self, x): u_inlet = Function(V) inlet_velocity = InletVelocity(t) u_inlet.interpolate(inlet_velocity) -bcu_inflow = dirichletbc(u_inlet, locate_dofs_topological(V, fdim, ft.find(inlet_marker))) +bcu_inflow = dirichletbc( + u_inlet, locate_dofs_topological(V, fdim, ft.find(inlet_marker)) +) # Walls u_nonslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType) -bcu_walls = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(wall_marker)), V) +bcu_walls = dirichletbc( + u_nonslip, locate_dofs_topological(V, fdim, ft.find(wall_marker)), V +) # Obstacle -bcu_obstacle = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(obstacle_marker)), V) +bcu_obstacle = dirichletbc( + u_nonslip, locate_dofs_topological(V, fdim, ft.find(obstacle_marker)), V +) bcu = [bcu_inflow, bcu_obstacle, bcu_walls] # Outlet -bcp_outlet = dirichletbc(PETSc.ScalarType(0), locate_dofs_topological(Q, fdim, ft.find(outlet_marker)), Q) +bcp_outlet = dirichletbc( + PETSc.ScalarType(0), locate_dofs_topological(Q, fdim, ft.find(outlet_marker)), Q +) bcp = [bcp_outlet] # - @@ -393,10 +433,11 @@ def __call__(self, x): # from pathlib import Path + folder = Path("results") folder.mkdir(exist_ok=True, parents=True) -vtx_u = VTXWriter(mesh.comm, "dfg2D-3-u.bp", [u_], engine="BP4") -vtx_p = VTXWriter(mesh.comm, "dfg2D-3-p.bp", [p_], engine="BP4") +vtx_u = VTXWriter(mesh.comm, folder / "dfg2D-3-u.bp", [u_], engine="BP4") +vtx_p = VTXWriter(mesh.comm, folder / "dfg2D-3-p.bp", [p_], engine="BP4") vtx_u.write(t) vtx_p.write(t) progress = tqdm.autonotebook.tqdm(desc="Solving PDE", total=num_steps) @@ -447,7 +488,11 @@ def __call__(self, x): vtx_p.write(t) # Update variable with solution form this time step - with u_.x.petsc_vec.localForm() as loc_, u_n.x.petsc_vec.localForm() as loc_n, u_n1.x.petsc_vec.localForm() as loc_n1: + with ( + u_.x.petsc_vec.localForm() as loc_, + u_n.x.petsc_vec.localForm() as loc_n, + u_n1.x.petsc_vec.localForm() as loc_n1, + ): loc_n.copy(loc_n1) loc_.copy(loc_n) @@ -496,31 +541,64 @@ def __call__(self, x): turek = np.loadtxt("bdforces_lv4") turek_p = np.loadtxt("pointvalues_lv4") fig = plt.figure(figsize=(25, 8)) - l1 = plt.plot(t_u, C_D, label=r"FEniCSx ({0:d} dofs)".format(num_velocity_dofs + num_pressure_dofs), linewidth=2) - l2 = plt.plot(turek[1:, 1], turek[1:, 3], marker="x", markevery=50, - linestyle="", markersize=4, label="FEATFLOW (42016 dofs)") + l1 = plt.plot( + t_u, + C_D, + label=r"FEniCSx ({0:d} dofs)".format(num_velocity_dofs + num_pressure_dofs), + linewidth=2, + ) + l2 = plt.plot( + turek[1:, 1], + turek[1:, 3], + marker="x", + markevery=50, + linestyle="", + markersize=4, + label="FEATFLOW (42016 dofs)", + ) plt.title("Drag coefficient") plt.grid() plt.legend() plt.savefig("figures/drag_comparison.png") fig = plt.figure(figsize=(25, 8)) - l1 = plt.plot(t_u, C_L, label=r"FEniCSx ({0:d} dofs)".format( - num_velocity_dofs + num_pressure_dofs), linewidth=2) - l2 = plt.plot(turek[1:, 1], turek[1:, 4], marker="x", markevery=50, - linestyle="", markersize=4, label="FEATFLOW (42016 dofs)") + l1 = plt.plot( + t_u, + C_L, + label=r"FEniCSx ({0:d} dofs)".format(num_velocity_dofs + num_pressure_dofs), + linewidth=2, + ) + l2 = plt.plot( + turek[1:, 1], + turek[1:, 4], + marker="x", + markevery=50, + linestyle="", + markersize=4, + label="FEATFLOW (42016 dofs)", + ) plt.title("Lift coefficient") plt.grid() plt.legend() plt.savefig("figures/lift_comparison.png") fig = plt.figure(figsize=(25, 8)) - l1 = plt.plot(t_p, p_diff, label=r"FEniCSx ({0:d} dofs)".format(num_velocity_dofs + num_pressure_dofs), linewidth=2) - l2 = plt.plot(turek[1:, 1], turek_p[1:, 6] - turek_p[1:, -1], marker="x", markevery=50, - linestyle="", markersize=4, label="FEATFLOW (42016 dofs)") + l1 = plt.plot( + t_p, + p_diff, + label=r"FEniCSx ({0:d} dofs)".format(num_velocity_dofs + num_pressure_dofs), + linewidth=2, + ) + l2 = plt.plot( + turek[1:, 1], + turek_p[1:, 6] - turek_p[1:, -1], + marker="x", + markevery=50, + linestyle="", + markersize=4, + label="FEATFLOW (42016 dofs)", + ) plt.title("Pressure difference") plt.grid() plt.legend() plt.savefig("figures/pressure_comparison.png") - - diff --git a/chapter3/component_bc.py b/chapter3/component_bc.py index a7b20855..0bd1afff 100644 --- a/chapter3/component_bc.py +++ b/chapter3/component_bc.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.4 +# jupytext_version: 1.16.5 # kernelspec: # display_name: Python 3 (ipykernel) # language: python diff --git a/chapter3/em.ipynb b/chapter3/em.ipynb index 12502202..64536f37 100644 --- a/chapter3/em.ipynb +++ b/chapter3/em.ipynb @@ -100,7 +100,7 @@ "from dolfinx.mesh import compute_midpoints, locate_entities_boundary\n", "from dolfinx.plot import vtk_mesh\n", "\n", - "from ufl import TestFunction, TrialFunction, as_vector, dot, dx, grad, inner\n", + "from ufl import TestFunction, TrialFunction, as_vector, dot, dx, grad\n", "from mpi4py import MPI\n", "\n", "import gmsh\n", @@ -216,7 +216,10 @@ }, "outputs": [], "source": [ - "mesh, ct, _ = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)\n", + "mesh_data = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)\n", + "mesh = mesh_data.mesh\n", + "assert mesh_data.cell_tags is not None\n", + "ct = mesh_data.cell_tags\n", "gmsh.finalize()" ] }, @@ -376,7 +379,7 @@ "source": [ "W = functionspace(mesh, (\"DG\", 0, (mesh.geometry.dim, )))\n", "B = Function(W)\n", - "B_expr = Expression(as_vector((A_z.dx(1), -A_z.dx(0))), W.element.interpolation_points())\n", + "B_expr = Expression(as_vector((A_z.dx(1), -A_z.dx(0))), W.element.interpolation_points)\n", "B.interpolate(B_expr)" ] }, diff --git a/chapter3/em.py b/chapter3/em.py index ec54d672..de37af0a 100644 --- a/chapter3/em.py +++ b/chapter3/em.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.4 +# jupytext_version: 1.16.5 # kernelspec: # display_name: Python 3 (ipykernel) # language: python @@ -96,7 +96,7 @@ from dolfinx.mesh import compute_midpoints, locate_entities_boundary from dolfinx.plot import vtk_mesh -from ufl import TestFunction, TrialFunction, as_vector, dot, dx, grad, inner +from ufl import TestFunction, TrialFunction, as_vector, dot, dx, grad from mpi4py import MPI import gmsh @@ -199,7 +199,10 @@ # As in [the Navier-Stokes tutorial](../chapter2/ns_code2) we load the mesh directly into DOLFINx, without writing it to file. -mesh, ct, _ = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2) +mesh_data = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2) +mesh = mesh_data.mesh +assert mesh_data.cell_tags is not None +ct = mesh_data.cell_tags gmsh.finalize() # To inspect the mesh, we use Paraview, and obtain the following mesh @@ -280,7 +283,7 @@ W = functionspace(mesh, ("DG", 0, (mesh.geometry.dim, ))) B = Function(W) -B_expr = Expression(as_vector((A_z.dx(1), -A_z.dx(0))), W.element.interpolation_points()) +B_expr = Expression(as_vector((A_z.dx(1), -A_z.dx(0))), W.element.interpolation_points) B.interpolate(B_expr) # Note that we used `ufl.as_vector` to interpret the `Python`-tuple `(A_z.dx(1), -A_z.dx(0))` as a vector in the unified form language (UFL). diff --git a/chapter3/multiple_dirichlet.py b/chapter3/multiple_dirichlet.py index 6dac99b0..628b41f3 100644 --- a/chapter3/multiple_dirichlet.py +++ b/chapter3/multiple_dirichlet.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.4 +# jupytext_version: 1.16.5 # kernelspec: # display_name: Python 3 (ipykernel) # language: python diff --git a/chapter3/neumann_dirichlet_code.py b/chapter3/neumann_dirichlet_code.py index dfb15636..e5ed35e1 100644 --- a/chapter3/neumann_dirichlet_code.py +++ b/chapter3/neumann_dirichlet_code.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.4 +# jupytext_version: 1.16.5 # kernelspec: # display_name: Python 3 (ipykernel) # language: python diff --git a/chapter3/robin_neumann_dirichlet.py b/chapter3/robin_neumann_dirichlet.py index 7cac576e..c2530a4c 100644 --- a/chapter3/robin_neumann_dirichlet.py +++ b/chapter3/robin_neumann_dirichlet.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.4 +# jupytext_version: 1.16.5 # kernelspec: # display_name: Python 3 (ipykernel) # language: python diff --git a/chapter3/subdomains.ipynb b/chapter3/subdomains.ipynb index d73d725f..7aaed31a 100644 --- a/chapter3/subdomains.ipynb +++ b/chapter3/subdomains.ipynb @@ -15,19 +15,28 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "lines_to_end_of_cell_marker": 2 + }, "outputs": [], "source": [ "from dolfinx import default_scalar_type\n", - "from dolfinx.fem import (Constant, dirichletbc, Function, functionspace, assemble_scalar,\n", - " form, locate_dofs_geometrical, locate_dofs_topological)\n", + "from dolfinx.fem import (\n", + " Constant,\n", + " dirichletbc,\n", + " Function,\n", + " functionspace,\n", + " assemble_scalar,\n", + " form,\n", + " locate_dofs_geometrical,\n", + " locate_dofs_topological,\n", + ")\n", "from dolfinx.fem.petsc import LinearProblem\n", "from dolfinx.io import XDMFFile, gmshio\n", "from dolfinx.mesh import create_unit_square, locate_entities\n", "from dolfinx.plot import vtk_mesh\n", "\n", - "from ufl import (SpatialCoordinate, TestFunction, TrialFunction,\n", - " dx, grad, inner)\n", + "from ufl import SpatialCoordinate, TestFunction, TrialFunction, dx, grad, inner\n", "\n", "from mpi4py import MPI\n", "\n", @@ -44,7 +53,9 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "lines_to_next_cell": 2 + }, "source": [ "We will use a simple example with two materials in two dimensions to demonstrate the idea. The whole domain will be $\\Omega=[0,1]\\times[0,1]$, which consists of two subdomains\n", "$\\Omega_0=[0,1]\\times [0,1/2]$ and $\\Omega_1=[0,1]\\times[1/2, 1]$. We start by creating two python functions, where each returns `True` if the input coordinate is inside its domain." @@ -153,7 +164,9 @@ }, "outputs": [], "source": [ - "problem = LinearProblem(a, L, bcs=bcs, petsc_options={\"ksp_type\": \"preonly\", \"pc_type\": \"lu\"})\n", + "problem = LinearProblem(\n", + " a, L, bcs=bcs, petsc_options={\"ksp_type\": \"preonly\", \"pc_type\": \"lu\"}\n", + ")\n", "uh = problem.solve()\n", "\n", "# Filter out ghosted cells\n", @@ -165,7 +178,9 @@ "marker[cells_0] = 1\n", "marker[cells_1] = 2\n", "mesh.topology.create_connectivity(tdim, tdim)\n", - "topology, cell_types, x = vtk_mesh(mesh, tdim, np.arange(num_cells_local, dtype=np.int32))\n", + "topology, cell_types, x = vtk_mesh(\n", + " mesh, tdim, np.arange(num_cells_local, dtype=np.int32)\n", + ")\n", "\n", "p = pyvista.Plotter(window_size=[800, 800])\n", "grid = pyvista.UnstructuredGrid(topology, cell_types, x)\n", @@ -180,7 +195,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "lines_to_next_cell": 2 + }, "outputs": [], "source": [ "p2 = pyvista.Plotter(window_size=[800, 800])\n", @@ -203,7 +220,9 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "lines_to_next_cell": 2 + }, "source": [ "## Interpolation with Python-function\n", "As we saw in the first approach, in many cases, we can use the geometrical coordinates to determine which coefficient we should use. Using the unstructured mesh from the previous example, we illustrate an alternative approach using interpolation:" @@ -250,7 +269,7 @@ "outputs": [], "source": [ "# Difference in kappa's\n", - "error = mesh.comm.allreduce(assemble_scalar(form((kappa - kappa2)**2 * dx)))\n", + "error = mesh.comm.allreduce(assemble_scalar(form((kappa - kappa2) ** 2 * dx)))\n", "print(error)" ] }, @@ -317,7 +336,12 @@ "metadata": {}, "outputs": [], "source": [ - "mesh, cell_markers, facet_markers = gmshio.read_from_msh(\"mesh.msh\", MPI.COMM_WORLD, gdim=2)" + "mesh_data = gmshio.read_from_msh(\"mesh.msh\", MPI.COMM_WORLD, gdim=2)\n", + "mesh = mesh_data.mesh\n", + "assert mesh_data.cell_tags is not None\n", + "cell_markers = mesh_data.cell_tags\n", + "assert mesh_data.facet_tags is not None\n", + "facet_markers = mesh_data.facet_tags" ] }, { @@ -347,7 +371,11 @@ " cells = mesh.get_cells_type(cell_type)\n", " cell_data = mesh.get_cell_data(\"gmsh:physical\", cell_type)\n", " points = mesh.points[:, :2] if prune_z else mesh.points\n", - " out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={\"name_to_read\": [cell_data.astype(np.int32)]})\n", + " out_mesh = meshio.Mesh(\n", + " points=points,\n", + " cells={cell_type: cells},\n", + " cell_data={\"name_to_read\": [cell_data.astype(np.int32)]},\n", + " )\n", " return out_mesh" ] }, @@ -461,7 +489,9 @@ "x = SpatialCoordinate(mesh)\n", "L = Constant(mesh, default_scalar_type(1)) * v * dx\n", "\n", - "problem = LinearProblem(a, L, bcs=bcs, petsc_options={\"ksp_type\": \"preonly\", \"pc_type\": \"lu\"})\n", + "problem = LinearProblem(\n", + " a, L, bcs=bcs, petsc_options={\"ksp_type\": \"preonly\", \"pc_type\": \"lu\"}\n", + ")\n", "uh = problem.solve()\n", "\n", "# As the dolfinx.MeshTag contains a value for every cell in the\n", @@ -499,13 +529,6 @@ "else:\n", " p2.screenshot(\"unstructured_u.png\")" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/chapter3/subdomains.py b/chapter3/subdomains.py index d6a09d96..2bf8b046 100644 --- a/chapter3/subdomains.py +++ b/chapter3/subdomains.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.4 +# jupytext_version: 1.16.5 # kernelspec: # display_name: Python 3 (ipykernel) # language: python @@ -22,15 +22,22 @@ # + from dolfinx import default_scalar_type -from dolfinx.fem import (Constant, dirichletbc, Function, functionspace, assemble_scalar, - form, locate_dofs_geometrical, locate_dofs_topological) +from dolfinx.fem import ( + Constant, + dirichletbc, + Function, + functionspace, + assemble_scalar, + form, + locate_dofs_geometrical, + locate_dofs_topological, +) from dolfinx.fem.petsc import LinearProblem from dolfinx.io import XDMFFile, gmshio from dolfinx.mesh import create_unit_square, locate_entities from dolfinx.plot import vtk_mesh -from ufl import (SpatialCoordinate, TestFunction, TrialFunction, - dx, grad, inner) +from ufl import SpatialCoordinate, TestFunction, TrialFunction, dx, grad, inner from mpi4py import MPI @@ -50,6 +57,7 @@ # We will use a simple example with two materials in two dimensions to demonstrate the idea. The whole domain will be $\Omega=[0,1]\times[0,1]$, which consists of two subdomains # $\Omega_0=[0,1]\times [0,1/2]$ and $\Omega_1=[0,1]\times[1/2, 1]$. We start by creating two python functions, where each returns `True` if the input coordinate is inside its domain. + # + def Omega_0(x): return x[1] <= 0.5 @@ -102,7 +110,9 @@ def Omega_1(x): # We can now solve and visualize the solution of the problem # + -problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"}) +problem = LinearProblem( + a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"} +) uh = problem.solve() # Filter out ghosted cells @@ -114,7 +124,9 @@ def Omega_1(x): marker[cells_0] = 1 marker[cells_1] = 2 mesh.topology.create_connectivity(tdim, tdim) -topology, cell_types, x = vtk_mesh(mesh, tdim, np.arange(num_cells_local, dtype=np.int32)) +topology, cell_types, x = vtk_mesh( + mesh, tdim, np.arange(num_cells_local, dtype=np.int32) +) p = pyvista.Plotter(window_size=[800, 800]) grid = pyvista.UnstructuredGrid(topology, cell_types, x) @@ -142,6 +154,7 @@ def Omega_1(x): # ## Interpolation with Python-function # As we saw in the first approach, in many cases, we can use the geometrical coordinates to determine which coefficient we should use. Using the unstructured mesh from the previous example, we illustrate an alternative approach using interpolation: + def eval_kappa(x): values = np.zeros(x.shape[1], dtype=default_scalar_type) # Create a boolean array indicating which dofs (corresponding to cell centers) @@ -159,7 +172,7 @@ def eval_kappa(x): # We verify this by assembling the error between this new function and the old one # Difference in kappa's -error = mesh.comm.allreduce(assemble_scalar(form((kappa - kappa2)**2 * dx))) +error = mesh.comm.allreduce(assemble_scalar(form((kappa - kappa2) ** 2 * dx))) print(error) # ## Subdomains defined from external mesh data @@ -202,7 +215,12 @@ def eval_kappa(x): # ## Read in MSH files with DOLFINx # You can read in MSH files with DOLFINx, which will read them in on a single process, and then distribute them over the available ranks in the MPI communicator. -mesh, cell_markers, facet_markers = gmshio.read_from_msh("mesh.msh", MPI.COMM_WORLD, gdim=2) +mesh_data = gmshio.read_from_msh("mesh.msh", MPI.COMM_WORLD, gdim=2) +mesh = mesh_data.mesh +assert mesh_data.cell_tags is not None +cell_markers = mesh_data.cell_tags +assert mesh_data.facet_tags is not None +facet_markers = mesh_data.facet_tags # ## Convert msh-files to XDMF using meshio # We will use `meshio` to read in the `msh` file, and convert it to a more suitable IO format. Meshio requires `h5py`, and can be installed on linux with the following commands: @@ -219,7 +237,11 @@ def create_mesh(mesh, cell_type, prune_z=False): cells = mesh.get_cells_type(cell_type) cell_data = mesh.get_cell_data("gmsh:physical", cell_type) points = mesh.points[:, :2] if prune_z else mesh.points - out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data.astype(np.int32)]}) + out_mesh = meshio.Mesh( + points=points, + cells={cell_type: cells}, + cell_data={"name_to_read": [cell_data.astype(np.int32)]}, + ) return out_mesh @@ -273,7 +295,9 @@ def create_mesh(mesh, cell_type, prune_z=False): x = SpatialCoordinate(mesh) L = Constant(mesh, default_scalar_type(1)) * v * dx -problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"}) +problem = LinearProblem( + a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"} +) uh = problem.solve() # As the dolfinx.MeshTag contains a value for every cell in the @@ -303,5 +327,3 @@ def create_mesh(mesh, cell_type, prune_z=False): p2.show() else: p2.screenshot("unstructured_u.png") - - diff --git a/chapter4/compiler_parameters.ipynb b/chapter4/compiler_parameters.ipynb index 77e42ab6..c067f1f7 100644 --- a/chapter4/compiler_parameters.ipynb +++ b/chapter4/compiler_parameters.ipynb @@ -50,7 +50,7 @@ "id": "2", "metadata": {}, "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. For codes using `dolfinx.LinearProblem`, you can supply `jit_options` as a keyword argument." + "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." ] }, { diff --git a/chapter4/compiler_parameters.py b/chapter4/compiler_parameters.py index e5811c4d..ff20725f 100644 --- a/chapter4/compiler_parameters.py +++ b/chapter4/compiler_parameters.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.4 +# jupytext_version: 1.16.5 # kernelspec: # display_name: Python 3 (ipykernel) # language: python @@ -46,7 +46,7 @@ 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.LinearProblem`, you can supply `jit_options` as a keyword argument. +# 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. # + diff --git a/chapter4/convergence.ipynb b/chapter4/convergence.ipynb index b2b20018..cf0babe4 100644 --- a/chapter4/convergence.ipynb +++ b/chapter4/convergence.ipynb @@ -137,7 +137,7 @@ " # is a ufl expression or a python lambda function\n", " u_ex_W = Function(W)\n", " if isinstance(u_ex, ufl.core.expr.Expr):\n", - " u_expr = Expression(u_ex, W.element.interpolation_points())\n", + " u_expr = Expression(u_ex, W.element.interpolation_points)\n", " u_ex_W.interpolate(u_expr)\n", " else:\n", " u_ex_W.interpolate(u_ex)\n", @@ -253,7 +253,7 @@ " comm = u_h.function_space.mesh.comm\n", " u_ex_V = Function(u_h.function_space)\n", " if isinstance(u_ex, ufl.core.expr.Expr):\n", - " u_expr = Expression(u_ex, u_h.function_space.element.interpolation_points())\n", + " u_expr = Expression(u_ex, u_h.function_space.element.interpolation_points)\n", " u_ex_V.interpolate(u_expr)\n", " else:\n", " u_ex_V.interpolate(u_ex)\n", diff --git a/chapter4/convergence.py b/chapter4/convergence.py index 24374792..dfa22927 100644 --- a/chapter4/convergence.py +++ b/chapter4/convergence.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.4 +# jupytext_version: 1.16.5 # kernelspec: # display_name: Python 3 (ipykernel) # language: python @@ -104,7 +104,7 @@ def error_L2(uh, u_ex, degree_raise=3): # is a ufl expression or a python lambda function u_ex_W = Function(W) if isinstance(u_ex, ufl.core.expr.Expr): - u_expr = Expression(u_ex, W.element.interpolation_points()) + u_expr = Expression(u_ex, W.element.interpolation_points) u_ex_W.interpolate(u_expr) else: u_ex_W.interpolate(u_ex) @@ -174,7 +174,7 @@ def error_infinity(u_h, u_ex): comm = u_h.function_space.mesh.comm u_ex_V = Function(u_h.function_space) if isinstance(u_ex, ufl.core.expr.Expr): - u_expr = Expression(u_ex, u_h.function_space.element.interpolation_points()) + u_expr = Expression(u_ex, u_h.function_space.element.interpolation_points) u_ex_V.interpolate(u_expr) else: u_ex_V.interpolate(u_ex) diff --git a/chapter4/newton-solver.ipynb b/chapter4/newton-solver.ipynb index 1d37295c..01c1ee34 100644 --- a/chapter4/newton-solver.ipynb +++ b/chapter4/newton-solver.ipynb @@ -564,18 +564,6 @@ "display_name": "Python 3 (ipykernel)", "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.10.12" } }, "nbformat": 4, diff --git a/chapter4/newton-solver.py b/chapter4/newton-solver.py index b76f771a..8c2f5e2f 100644 --- a/chapter4/newton-solver.py +++ b/chapter4/newton-solver.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.4 +# jupytext_version: 1.16.5 # kernelspec: # display_name: Python 3 (ipykernel) # language: python diff --git a/chapter4/solvers.py b/chapter4/solvers.py index ed460c7b..ceee9393 100644 --- a/chapter4/solvers.py +++ b/chapter4/solvers.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.4 +# jupytext_version: 1.16.5 # kernelspec: # display_name: Python 3 (ipykernel) # language: python diff --git a/fem.md b/fem.md index f2f86710..14ad6afd 100644 --- a/fem.md +++ b/fem.md @@ -40,13 +40,13 @@ The tutorial uses several dependencies for meshing, plotting and timings. A comp To use the notebooks in this tutorial with DOLFINx on your own computer, you should use the docker image obtained using the following command: ```bash - docker run --init -p 8888:8888 -v "$(pwd)":/root/shared ghcr.io/jorgensd/dolfinx-tutorial:v0.7.2 + docker run --init -p 8888:8888 -v "$(pwd)":/root/shared ghcr.io/jorgensd/dolfinx-tutorial:release ``` This image can also be used as a normal docker container by adding: ```bash - docker run --ti -v "$(pwd)":/root/shared --entrypoint="/bin/bash" ghcr.io/jorgensd/dolfinx-tutorial:v0.7.2 + docker run --ti -v "$(pwd)":/root/shared --entrypoint="/bin/bash" ghcr.io/jorgensd/dolfinx-tutorial:release ``` The tutorials can also be exported as an IPython notebook or PDF by clicking the ![Download](save.png)-symbol in the top right corner of the relevant tutorial. The notebook can in turn be used with a Python kernel which has DOLFINx. @@ -58,7 +58,7 @@ The [Dockerfile](https://github.com/FEniCS/dolfinx/blob/main/docker/Dockerfile) provides a definitive build recipe. As the DOLFINx docker images are hosted at Docker-hub, one can directly access the image using: ```bash -docker run dolfinx/dolfinx:v0.8.0 +docker run dolfinx/dolfinx:stable ``` There are several ways of customizing a docker container, such as mounting volumes/sharing folder, setting a working directory, sharing graphical interfaces etc. See `docker run --help` for an extensive list.