|
6 | 6 | # extension: .py
|
7 | 7 | # format_name: light
|
8 | 8 | # format_version: '1.5'
|
9 |
| -# jupytext_version: 1.14.1 |
| 9 | +# jupytext_version: 1.14.4 |
10 | 10 | # kernelspec:
|
11 | 11 | # display_name: Python 3 (ipykernel)
|
12 | 12 | # language: python
|
|
19 | 19 | # In this section, we will solve the deflection of the membrane problem.
|
20 | 20 | # After finishing this section, you should be able to:
|
21 | 21 | # - Create a simple mesh using the GMSH Python API and load it into DOLFINx
|
22 |
| -# - How to create a constant boundary conditions using a geometrical identifier |
23 |
| -# - Using s`ufl.SpatialCoordinate` to create a spatially varying function |
24 |
| -# - How to interpolate a `ufl.Expression` into an appropriate function space |
25 |
| -# - How to evaluate a `dolfinx.Function` at any point $x$ |
| 22 | +# - Create constant boundary conditions using a geometrical identifier |
| 23 | +# - Use `ufl.SpatialCoordinate` to create a spatially varying function |
| 24 | +# - Interpolate a `ufl.Expression` into an appropriate function space |
| 25 | +# - Evaluate a `dolfinx.Function` at any point $x$ |
26 | 26 | # - Use Paraview to visualize the solution of a PDE
|
27 | 27 | #
|
28 | 28 | # ## Creating the mesh
|
29 | 29 | #
|
30 |
| -# To create the computational geometry, we use the python-API of [GMSH](http://gmsh.info/). We start by import the gmsh-module, and initalizing it. |
| 30 | +# To create the computational geometry, we use the python-API of [GMSH](http://gmsh.info/). We start by importing the gmsh-module and initializing it. |
31 | 31 |
|
32 | 32 | import gmsh
|
33 | 33 | gmsh.initialize()
|
34 | 34 |
|
35 |
| -# The next step is to create the membrane and starting the computations by the GMSH CAD kernel, to generate the relevant underlying data structures. The arguments into `addDisk` is the x, y and z coordinate of the center of the circle, while the to last argument is the x-radius and y-radius. |
| 35 | +# 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. |
36 | 36 |
|
37 | 37 | membrane = gmsh.model.occ.addDisk(0, 0, 0, 1, 1)
|
38 | 38 | gmsh.model.occ.synchronize()
|
39 | 39 |
|
40 |
| -# The next step is to 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 two as the first argument, the entity tag of the membrane as the second argument, and the last argument is the physical tag. In a later demo, we will get into when this tag matters. |
| 40 | +# 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. |
41 | 41 |
|
42 | 42 | gdim = 2
|
43 | 43 | gmsh.model.addPhysicalGroup(gdim, [membrane], 1)
|
44 | 44 |
|
45 |
| -# Finally, we generate the two-dimensional mesh. We set a uniform mesh size by modifying the GMSH options |
| 45 | +# Finally, we generate the two-dimensional mesh. We set a uniform mesh size by modifying the GMSH options. |
46 | 46 |
|
47 | 47 | gmsh.option.setNumber("Mesh.CharacteristicLengthMin",0.05)
|
48 | 48 | gmsh.option.setNumber("Mesh.CharacteristicLengthMax",0.05)
|
49 | 49 | gmsh.model.mesh.generate(gdim)
|
50 | 50 |
|
51 | 51 | # # Interfacing with GMSH in DOLFINx
|
52 |
| -# We will import the GMSH-mesh directly from GMSH, using the `dolfinx.io.gmshio` interface in DOLFINx. |
53 |
| -# As we in the example 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. We therefore 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 not add any physical groups of dimension `gdim-1`, there will be no entities in the `facet_markers`. |
| 52 | +# We will import the GMSH-mesh directly from GMSH into DOLFINx via the `dolfinx.io.gmshio` interface. |
| 53 | +# 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`. |
54 | 54 |
|
55 | 55 | # +
|
56 | 56 | from dolfinx.io import gmshio
|
|
67 | 67 | V = fem.FunctionSpace(domain, ("CG", 1))
|
68 | 68 |
|
69 | 69 | # ## Defining a spatially varying load
|
70 |
| -# The right hand side pressure function is represented using `ufl.SpatialCoordinate` and a two constants, one for $\beta$ and one for $R_0$. |
| 70 | +# The right hand side pressure function is represented using `ufl.SpatialCoordinate` and two constants, one for $\beta$ and one for $R_0$. |
71 | 71 |
|
72 | 72 | import ufl
|
73 | 73 | from petsc4py.PETSc import ScalarType
|
|
77 | 77 | p = 4 * ufl.exp(-beta**2 * (x[0]**2 + (x[1] - R0)**2))
|
78 | 78 |
|
79 | 79 | # ## Create a Dirichlet boundary condition using geometrical conditions
|
80 |
| -# The next step is to create the homogenous 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$. |
| 80 | +# 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$. |
81 | 81 |
|
82 | 82 | import numpy as np
|
83 | 83 | def on_boundary(x):
|
84 | 84 | return np.isclose(np.sqrt(x[0]**2 + x[1]**2), 1)
|
85 | 85 | boundary_dofs = fem.locate_dofs_geometrical(V, on_boundary)
|
86 | 86 |
|
87 |
| -# As our Dirichlet condition is homogenous (`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. |
| 87 | +# 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. |
88 | 88 |
|
89 | 89 | bc = fem.dirichletbc(ScalarType(0), boundary_dofs, V)
|
90 | 90 |
|
@@ -159,7 +159,7 @@ def on_boundary(x):
|
159 | 159 | u_values = []
|
160 | 160 | p_values = []
|
161 | 161 |
|
162 |
| -# 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$, $\phi_i$ the $i$th basis function, we can compute the exact solution at any point in $\Omega$. |
| 162 | +# 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$. |
163 | 163 | # 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.
|
164 | 164 | # 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.
|
165 | 165 |
|
@@ -190,7 +190,7 @@ def on_boundary(x):
|
190 | 190 | u_values = uh.eval(points_on_proc, cells)
|
191 | 191 | p_values = pressure.eval(points_on_proc, cells)
|
192 | 192 |
|
193 |
| -# As we now have an array of coordinates and two arrays of function values, we use matplotlib to plot them |
| 193 | +# As we now have an array of coordinates and two arrays of function values, we can use `matplotlib` to plot them |
194 | 194 |
|
195 | 195 | import matplotlib.pyplot as plt
|
196 | 196 | fig = plt.figure()
|
|
0 commit comments