You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Copy file name to clipboardExpand all lines: chapter1/membrane_code.py
+79-33Lines changed: 79 additions & 33 deletions
Original file line number
Diff line number
Diff line change
@@ -19,58 +19,70 @@
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
-
# - 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.fem.Function` at any point $x$
22
+
# - Create constant boundary conditions using a {py:func}`geometrical identifier<dolfinx.fem.locate_dofs_geometrical>`
23
+
# - Use {py:class}`ufl.SpatialCoordinate` to create a spatially varying function
24
+
# - Interpolate a {py:class}`ufl-Expression<ufl.core.expr.Expr>` into an appropriate function space
25
+
# - Evaluate a {py:class}`dolfinx.fem.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](https://gmsh.info/). We start by importing the gmsh-module and initializing it.
30
+
# To create the computational geometry, we use the Python-API of [GMSH](https://gmsh.info/).
31
+
# We start by importing the gmsh-module and initializing it.
31
32
32
33
# +
33
34
importgmsh
34
35
35
36
gmsh.initialize()
36
37
# -
37
38
38
-
# 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.
39
+
# The next step is to create the membrane and start the computations by the GMSH CAD kernel,
40
+
# to generate the relevant underlying data structures.
41
+
# The first arguments of `addDisk` are the x, y and z coordinate of the center of the circle,
42
+
# while the two last arguments are the x-radius and y-radius.
39
43
40
44
membrane=gmsh.model.occ.addDisk(0, 0, 0, 1, 1)
41
45
gmsh.model.occ.synchronize()
42
46
43
-
# 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.
47
+
# After that, we make the membrane a physical surface, such that it is recognized by `gmsh` when generating the mesh.
48
+
# As a surface is a two-dimensional entity, we add `2` as the first argument,
49
+
# the entity tag of the membrane as the second argument, and the physical tag as the last argument.
50
+
# In a later demo, we will get into when this tag matters.
44
51
45
52
gdim=2
46
53
gmsh.model.addPhysicalGroup(gdim, [membrane], 1)
47
54
48
-
# Finally, we generate the two-dimensional mesh. We set a uniform mesh size by modifying the GMSH options.
55
+
# Finally, we generate the two-dimensional mesh.
56
+
# We set a uniform mesh size by modifying the GMSH options.
# ## Create a Dirichlet boundary condition using geometrical conditions
109
-
# 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$.
124
+
# The next step is to create the homogeneous boundary condition.
125
+
# As opposed to the [first tutorial](./fundamentals_code.ipynb) we will use
126
+
# {py:func}`locate_dofs_geometrical<dolfinx.fem.locate_dofs_geometrical>` to locate the degrees of freedom on the boundary.
127
+
# As we know that our domain is a circle with radius 1, we know that any degree of freedom should be
128
+
# located at a coordinate $(x,y)$ such that $\sqrt{x^2+y^2}=1$.
# 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.
141
+
# As our Dirichlet condition is homogeneous (`u=0` on the whole boundary), we can initialize the
142
+
# {py:class}`dolfinx.fem.DirichletBC` with a constant value, the degrees of freedom and the function
143
+
# space to apply the boundary condition on. We use the constructor {py:func}`dolfinx.fem.dirichletbc`.
# 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.
163
+
# ## Interpolation of a UFL-expression
164
+
# As we previously defined the load `p` as a spatially varying function,
165
+
# we would like to interpolate this function into an appropriate function space for visualization.
166
+
# To do this we use the class {py:class}`Expression<dolfinx.fem.Expression>`.
167
+
# The expression takes in any UFL-expression, and a set of points on the reference element.
168
+
# We will use the {py:attr}`interpolation points<dolfinx.fem.FiniteElement.interpolation_points>`
169
+
# of the space we want to interpolate in to.
144
170
# We choose a high order function space to represent the function `p`, as it is rapidly varying in space.
145
171
146
172
Q=fem.functionspace(domain, ("Lagrange", 5))
@@ -154,8 +180,7 @@ def on_boundary(x):
154
180
fromdolfinx.plotimportvtk_mesh
155
181
importpyvista
156
182
157
-
158
-
# Extract topology from mesh and create pyvista mesh
183
+
# Extract topology from mesh and create {py:class}`pyvista.UnstructuredGrid`
# Another way to compare the deflection and the load is to make a plot along the line $x=0$.
193
-
# 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.
218
+
# This is just a matter of defining a set of points along the $y$-axis and evaluating the
219
+
# finite element functions $u$ and $p$ at these points.
194
220
195
221
tol=0.001# Avoid hitting the outside of the domain
196
222
y=np.linspace(-1+tol, 1-tol, 101)
@@ -199,22 +225,40 @@ def on_boundary(x):
199
225
u_values= []
200
226
p_values= []
201
227
202
-
# 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$.
203
-
# 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.
204
-
# 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.
228
+
# As a finite element function is the linear combination of all degrees of freedom,
229
+
# $u_h(x)=\sum_{i=1}^N c_i \phi_i(x)$ where $c_i$ are the coefficients of $u_h$ and $\phi_i$
230
+
# is the $i$-th basis function, we can compute the exact solution at any point in $\Omega$.
231
+
# However, as a mesh consists of a large set of degrees of freedom (i.e. $N$ is large),
232
+
# we want to reduce the number of evaluations of the basis function $\phi_i(x)$.
233
+
# We do this by identifying which cell of the mesh $x$ is in.
234
+
# This is efficiently done by creating a {py:class}`bounding box tree<dolfinx.geometry.BoundingBoxTree`
235
+
# of the cells of the mesh,
236
+
# allowing a quick recursive search through the mesh entities.
# 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)`.
213
-
# 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
214
-
# using `dolfinx.geometry.select_colliding_cells`, which measures the exact distance between the point and the cell (approximated as a convex hull for higher order geometries).
215
-
# 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.
244
+
# Now we can compute which cells the bounding box tree collides with using
# which measures the exact distance between the point and the cell
254
+
# (approximated as a convex hull for higher order geometries).
255
+
# This function also returns an adjacency-list, as the point might align with a facet,
256
+
# edge or vertex that is shared between multiple cells in the mesh.
216
257
#
217
-
# 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.
258
+
# Finally, we would like the code below to run in parallel,
259
+
# when the mesh is distributed over multiple processors.
260
+
# In that case, it is not guaranteed that every point in `points` is on each processor.
261
+
# Therefore we create a subset `points_on_proc` only containing the points found on the current processor.
218
262
219
263
cells= []
220
264
points_on_proc= []
@@ -227,14 +271,16 @@ def on_boundary(x):
227
271
points_on_proc.append(point)
228
272
cells.append(colliding_cells.links(i)[0])
229
273
230
-
# 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.
231
-
#
274
+
# We now have a list of points on the processor, on in which cell each point belongs.
275
+
# We can then call {py:meth}`uh.eval<dolfinx.fem.Function.eval>` and
276
+
# {py:meth}`pressure.eval<dolfinx.fem.Function.eval>` to obtain the set of values for all the points.
0 commit comments