Skip to content
Merged
20 changes: 12 additions & 8 deletions python/demo/demo_axis.py
Original file line number Diff line number Diff line change
Expand Up @@ -154,8 +154,8 @@ def generate_mesh_sphere_axis(
# \cdot (\nabla \times \bar{\mathbf{v}})+\varepsilon_{r} k_{0}^{2}
# \mathbf{E}_s \cdot \bar{\mathbf{v}}+k_{0}^{2}\left(\varepsilon_{r}
# -\varepsilon_b\right)\mathbf{E}_b \cdot \bar{\mathbf{v}}~\mathrm{d} x\\
# +\int_{\Omega_{pml}}\left[\boldsymbol{\mu}^{-1}_{pml} \nabla \times \mathbf{E}_s
# \right]\cdot \nabla \times \bar{\mathbf{v}}-k_{0}^{2}
# +\int_{\Omega_{pml}}\left[\boldsymbol{\mu}^{-1}_{pml} \nabla \times
# \mathbf{E}_s \right]\cdot \nabla \times \bar{\mathbf{v}}-k_{0}^{2}
# \left[\boldsymbol{\varepsilon}_{pml} \mathbf{E}_s \right]\cdot
# \bar{\mathbf{v}}~ d x=0
# \end{split}
Expand All @@ -172,8 +172,8 @@ def generate_mesh_sphere_axis(
# -\varepsilon_b\right)\mathbf{E}_b \cdot
# \bar{\mathbf{v}}~ \rho d\rho dz d \phi\\
# +\int_{\Omega_{cs}}
# \int_{0}^{2\pi}\left[\boldsymbol{\mu}^{-1}_{pml} \nabla \times \mathbf{E}_s
# \right]\cdot \nabla \times \bar{\mathbf{v}}-k_{0}^{2}
# \int_{0}^{2\pi}\left[\boldsymbol{\mu}^{-1}_{pml} \nabla \times
# \mathbf{E}_s \right]\cdot \nabla \times \bar{\mathbf{v}}-k_{0}^{2}
# \left[\boldsymbol{\varepsilon}_{pml} \mathbf{E}_s \right]\cdot
# \bar{\mathbf{v}}~ \rho d\rho dz d \phi=0
# \end{split}
Expand All @@ -184,8 +184,10 @@ def generate_mesh_sphere_axis(
#
# $$
# \begin{align}
# \mathbf{E}_s(\rho, z, \phi) &= \sum_m\mathbf{E}^{(m)}_s(\rho, z)e^{-jm\phi} \\
# \mathbf{E}_b(\rho, z, \phi) &= \sum_m\mathbf{E}^{(m)}_b(\rho, z)e^{-jm\phi} \\
# \mathbf{E}_s(\rho, z, \phi) &= \sum_m\mathbf{E}^{(m)}_s(\rho, z)
# e^{-jm\phi} \\
# \mathbf{E}_b(\rho, z, \phi) &= \sum_m\mathbf{E}^{(m)}_b(\rho, z)
# e^{-jm\phi} \\
# \bar{\mathbf{v}}(\rho, z, \phi) &=
# \sum_m\bar{\mathbf{v}}^{(m)}(\rho, z)e^{+jm\phi}
# \end{align}
Expand Down Expand Up @@ -749,8 +751,10 @@ def create_eps_mu(pml, rho, eps_bkg, mu_bkg):
# m = np.sqrt(eps_au)/n_bkg
# x = 2*np.pi*radius_sph/wl0*n_bkg
#
# q_sca_analyt, q_abs_analyt = scattnlay(np.array([x], dtype=np.complex128),
# np.array([m], dtype=np.complex128))[2:4]
# q_sca_analyt, q_abs_analyt = scattnlay(
# np.array([x], dtype=np.complex128),
# np.array([m], dtype=np.complex128)
# )[2:4]
# ```
#
# The numerical values are reported here below:
Expand Down
46 changes: 24 additions & 22 deletions python/demo/demo_biharmonic.py
Original file line number Diff line number Diff line change
Expand Up @@ -33,23 +33,23 @@
# \nabla^{4} u = f \quad {\rm in} \ \Omega,
# $$
#
# where $\nabla^{4} \equiv \nabla^{2} \nabla^{2}$ is the biharmonic operator
# and $f$ is a prescribed source term.
# where $\nabla^{4} \equiv \nabla^{2} \nabla^{2}$ is the biharmonic
# operator and $f$ is a prescribed source term.
# To formulate a complete boundary value problem, the biharmonic equation
# must be complemented by suitable boundary conditions.
#
# ### Weak formulation
#
# Multiplying the biharmonic equation by a test function and integrating
# by parts twice leads to a problem of second-order derivatives, which would
# require $H^{2}$ conforming (roughly $C^{1}$ continuous) basis functions.
# To solve the biharmonic equation using Lagrange finite element basis
# functions, the biharmonic equation can be split into two second-order
# equations (see the Mixed Poisson demo for a mixed method for the Poisson
# equation), or a variational formulation can be constructed that imposes
# weak continuity of normal derivatives between finite element cells.
# This demo uses a discontinuous Galerkin approach to impose continuity
# of the normal derivative weakly.
# by parts twice leads to a problem of second-order derivatives, which
# would require $H^{2}$ conforming (roughly $C^{1}$ continuous) basis
# functions. To solve the biharmonic equation using Lagrange finite element
# basis functions, the biharmonic equation can be split into two second-
# order equations (see the Mixed Poisson demo for a mixed method for the
# Poisson equation), or a variational formulation can be constructed that
# imposes weak continuity of normal derivatives between finite element
# cells. This demo uses a discontinuous Galerkin approach to impose
# continuity of the normal derivative weakly.
#
# Consider a triangulation $\mathcal{T}$ of the domain $\Omega$, where
# the set of interior facets is denoted by $\mathcal{E}_h^{\rm int}$.
Expand All @@ -71,7 +71,8 @@
# \end{align}
# $$
#
# a weak formulation of the biharmonic problem reads: find $u \in V$ such that
# a weak formulation of the biharmonic problem reads: find $u \in V$ such
# that
#
# $$
# a(u,v)=L(v) \quad \forall \ v \in V,
Expand Down Expand Up @@ -158,8 +159,9 @@
# Next, we locate the mesh facets that lie on the boundary
# $\Gamma_D = \partial\Omega$.
# We do this using using {py:func}`exterior_facet_indices
# <dolfinx.mesh.exterior_facet_indices>` which returns all mesh boundary facets
# (Note: if we are only interested in a subset of those, consider {py:func}`locate_entities_boundary
# <dolfinx.mesh.exterior_facet_indices>` which returns all mesh boundary
# facets (Note: if we are only interested in a subset of those, consider
# {py:func}`locate_entities_boundary
# <dolfinx.mesh.locate_entities_boundary>`).

tdim = msh.topology.dim
Expand All @@ -182,8 +184,8 @@

# Next, we express the variational problem using UFL.
#
# First, the penalty parameter $\alpha$ is defined. In addition, we define a
# variable `h` for the cell diameter $h_E$, a variable `n`for the
# First, the penalty parameter $\alpha$ is defined. In addition, we define
# a variable `h` for the cell diameter $h_E$, a variable `n`for the
# outward-facing normal vector $n$ and a variable `h_avg` for the
# average size of cells sharing a facet
# $\left< h \right> = \frac{1}{2} (h_{+} + h_{-})$. Here, the UFL syntax
Expand All @@ -195,12 +197,12 @@
n = ufl.FacetNormal(msh)
h_avg = (h("+") + h("-")) / 2.0

# After that, we can define the variational problem consisting of the bilinear
# form $a$ and the linear form $L$. The source term is prescribed as
# $f = 4.0 \pi^4\sin(\pi x)\sin(\pi y)$. Note that with `dS`, integration is
# carried out over all the interior facets $\mathcal{E}_h^{\rm int}$, whereas
# with `ds` it would be only the facets on the boundary of the domain, i.e.
# $\partial\Omega$. The jump operator
# After that, we can define the variational problem consisting of the
# bilinear form $a$ and the linear form $L$. The source term is prescribed
# as $f = 4.0 \pi^4\sin(\pi x)\sin(\pi y)$. Note that with `dS`,
# integration is carried out over all the interior facets
# $\mathcal{E}_h^{\rm int}$, whereas with `ds` it would be only the facets
# on the boundary of the domain, i.e. $\partial\Omega$. The jump operator
# $[\!\![ w ]\!\!] = w_{+} \cdot n_{+} + w_{-} \cdot n_{-}$ w.r.t. the
# outward-facing normal vector $n$ is in UFL available as `jump(w, n)`.

Expand Down
32 changes: 17 additions & 15 deletions python/demo/demo_cahn-hilliard.py
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,8 @@
# \begin{align}
# \frac{\partial c}{\partial t} - \nabla \cdot M \nabla\mu
# &= 0 \quad {\rm in} \ \Omega, \\
# \mu - \frac{d f}{d c} + \lambda \nabla^{2}c &= 0 \quad {\rm in} \ \Omega.
# \mu - \frac{d f}{d c} + \lambda \nabla^{2}c &= 0 \quad {\rm in}
# \ \Omega.
# \end{align}
# $$
#
Expand All @@ -71,8 +72,8 @@
# \int_{\Omega} \frac{\partial c}{\partial t} q \, {\rm d} x +
# \int_{\Omega} M \nabla\mu \cdot \nabla q \, {\rm d} x
# &= 0 \quad \forall \ q \in V, \\
# \int_{\Omega} \mu v \, {\rm d} x - \int_{\Omega} \frac{d f}{d c} v \, {\rm d} x
# - \int_{\Omega} \lambda \nabla c \cdot \nabla v \, {\rm d} x
# \int_{\Omega} \mu v \, {\rm d} x - \int_{\Omega} \frac{d f}{d c} v \,
# {\rm d} x - \int_{\Omega} \lambda \nabla c \cdot \nabla v \, {\rm d} x
# &= 0 \quad \forall \ v \in V.
# \end{align}
# $$
Expand All @@ -86,15 +87,16 @@
# $$
# \begin{align}
# \int_{\Omega} \frac{c_{n+1} - c_{n}}{dt} q \, {\rm d} x
# + \int_{\Omega} M \nabla \mu_{n+\theta} \cdot \nabla q \, {\rm d} x
# &= 0 \quad \forall \ q \in V \\
# \int_{\Omega} \mu_{n+1} v \, {\rm d} x - \int_{\Omega} \frac{d f_{n+1}}{d c} v \, {\rm d} x
# - \int_{\Omega} \lambda \nabla c_{n+1} \cdot \nabla v \, {\rm d} x
# &= 0 \quad \forall \ v \in V
# + \int_{\Omega} M \nabla \mu_{n+\theta} \cdot \nabla q \, {\rm d} x
# &= 0 \quad \forall \ q \in V \\
# \int_{\Omega} \mu_{n+1} v \, {\rm d} x - \int_{\Omega}
# \frac{d f_{n+1}}{d c} v \, {\rm d} x - \int_{\Omega} \lambda \nabla
# c_{n+1} \cdot \nabla v \, {\rm d} x &= 0 \quad \forall \ v \in V
# \end{align}
# $$
#
# where $dt = t_{n+1} - t_{n}$ and $\mu_{n+\theta} = (1-\theta) \mu_{n} + \theta \mu_{n+1}$.
# where $dt = t_{n+1} - t_{n}$ and $\mu_{n+\theta} = (1-\theta) \mu_{n} +
# \theta \mu_{n+1}$.
# The task is: given $c_{n}$ and $\mu_{n}$, solve the above equation to
# find $c_{n+1}$ and $\mu_{n+1}$.
#
Expand Down Expand Up @@ -263,14 +265,14 @@
# ```
#
# To solve the nonlinear system of equations,
# {py:class}`NonlinearProblem<dolfinx.fem.petsc.NonlinearProblem>` object to
# solve a system of nonlinear equations
# {py:class}`NonlinearProblem<dolfinx.fem.petsc.NonlinearProblem>` object
# to solve a system of nonlinear equations

# +
# For the factorisation of the underlying linearized problems, prefer MUMPS,
# then superlu_dist, then default.
# We measure convergence by looking at the norm of the increment of the solution
# between two iterations, called `stol` in PETSc, see:
# For the factorisation of the underlying linearized problems, prefer
# MUMPS, then superlu_dist, then default.
# We measure convergence by looking at the norm of the increment of the
# solution between two iterations, called `stol` in PETSc, see:
# [`SNES convegence tests`](https://petsc.org/release/manual/snes/#convergence-tests)
# for further details.

Expand Down
3 changes: 2 additions & 1 deletion python/demo/demo_gmsh.py
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,8 @@ def gmsh_sphere(model: gmsh.model, name: str) -> gmsh.model:
# Add physical tag for sphere
model.add_physical_group(dim=3, tags=[sphere], tag=1)

# Embed all sub-entities from the GMSH model into the sphere and tag them
# Embed all sub-entities from the GMSH model into the sphere and tag
# them
for dim in [0, 1, 2]:
entities = model.getEntities(dim)
entity_ids = [entity[1] for entity in entities]
Expand Down
4 changes: 2 additions & 2 deletions python/demo/demo_hdg.py
Original file line number Diff line number Diff line change
Expand Up @@ -103,8 +103,8 @@ def u_e(x):
facets = np.arange(num_facets, dtype=np.int32)

# Create the sub-mesh
# NOTE Despite all facets being present in the submesh, the entity map isn't
# necessarily the identity in parallel
# NOTE Despite all facets being present in the submesh, the entity map
# isn't necessarily the identity in parallel
facet_mesh, facet_mesh_to_mesh = mesh.create_submesh(msh, fdim, facets)[:2]

# Define function spaces
Expand Down
8 changes: 5 additions & 3 deletions python/demo/demo_lagrange_variants.py
Original file line number Diff line number Diff line change
Expand Up @@ -37,8 +37,9 @@
# influence of interpolation point position, we create a degree 10
# element on an interval using equally spaced points, and plot the basis
# functions. We create this element using `basix.ufl`'s
# `element` function. The function `element.tabulate` returns a 3-dimensional
# array with shape (derivatives, points, (value size) * (basis functions)).
# `element` function. The function `element.tabulate` returns a 3-
# dimensional array with shape (derivatives, points, (value size) *
# (basis functions)).
# In this example, we only tabulate the 0th derivative and the value
# size is 1, so we take the slice `[0, :, :]` to get a 2-dimensional
# array.
Expand Down Expand Up @@ -158,7 +159,8 @@ def saw_tooth(x):
# Lagrange compared to the GLL variant. To quantify the error, we
# compute the interpolation error in the $L_2$ norm,
#
# $$\left\|u - u_h\right\|_2 = \left(\int_0^1 (u - u_h)^2\right)^{\frac{1}{2}},$$
# $$\left\|u - u_h\right\|_2 =
# \left(\int_0^1 (u - u_h)^2\right)^{\frac{1}{2}},$$
#
# where $u$ is the function and $u_h$ is its interpolation in the finite
# element space. The following code uses UFL to compute the $L_2$ error
Expand Down
2 changes: 1 addition & 1 deletion python/demo/demo_mixed-poisson.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
# jupytext_version: 1.14.1
# ---

# # Mixed formulation of the Poisson equation with a block-preconditioner/solver
# # Mixed formulation of the Poisson equation with a block-preconditioner/solver # noqa
#
# This demo illustrates how to solve the Poisson equation using a mixed
# (two-field) formulation and a block-preconditioned iterative solver.
Expand Down
4 changes: 2 additions & 2 deletions python/demo/demo_mixed-topology.py
Original file line number Diff line number Diff line change
Expand Up @@ -110,8 +110,8 @@
V_cpp = _cpp.fem.FunctionSpace_float64(mesh, elements_cpp, dofmaps)

# Create forms for each cell type.
# FIXME This hack is required at the moment because UFL does not yet know about
# mixed topology meshes.
# FIXME This hack is required at the moment because UFL does not yet know
# about mixed topology meshes.
a = []
L = []
for i, cell_name in enumerate(["hexahedron", "prism"]):
Expand Down
18 changes: 11 additions & 7 deletions python/demo/demo_navier-stokes.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
# jupytext_version: 1.13.6
# ---

# # Divergence conforming discontinuous Galerkin method for the Navier--Stokes equations
# # Divergence conforming discontinuous Galerkin method for the Navier--Stokes equations # noqa
#
# This demo ({download}`demo_navier-stokes.py`) illustrates how to
# implement a divergence conforming discontinuous Galerkin method for
Expand All @@ -29,7 +29,8 @@
#
# $$
# \begin{align}
# \partial_t u - \nu \Delta u + (u \cdot \nabla)u + \nabla p &= f \text{ in } \Omega_t, \\
# \partial_t u - \nu \Delta u + (u \cdot \nabla)u + \nabla p &= f
# \text{ in } \Omega_t, \\
# \nabla \cdot u &= 0 \text{ in } \Omega_t,
# \end{align}
# $$
Expand Down Expand Up @@ -124,16 +125,17 @@
# \end{cases}
# $$
#
# where $\Gamma^\psi = \left\{x \in \Gamma; \; \psi(x) \cdot n(x) < 0\right\}$.
# where $\Gamma^\psi = \left\{x \in \Gamma; \; \psi(x) \cdot n(x) < 0
# \right\}$.
#
# The semi-discrete version problem (in dimensionless form) is: find
# $(u_h, p_h) \in V_h^{u_D} \times Q_h$ such that
#
# $$
# \begin{align}
# \int_\Omega \partial_t u_h \cdot v + a_h(u_h, v_h) + c_h(u_h; u_h, v_h)
# + b_h(v_h, p_h) &= \int_\Omega f \cdot v_h + L_{a_h}(v_h) + L_{c_h}(v_h)
# \quad \forall v_h \in V_h^0, \\
# + b_h(v_h, p_h) &= \int_\Omega f \cdot v_h + L_{a_h}(v_h) +
# L_{c_h}(v_h) \quad \forall v_h \in V_h^0, \\
# b_h(u_h, q_h) &= 0 \quad \forall q_h \in Q_h,
# \end{align}
# $$
Expand All @@ -151,8 +153,10 @@
# c_h(w; u, v) &= - \sumK \int_K u \cdot \nabla \cdot (v \otimes w)
# + \sumK \int_{\partial_K} w \cdot n \hat{u}^{w} \cdot v, \\
# L_{a_h}(v_h) &= Re^{-1} \left(- \int_{\partial \Omega} u_D \otimes n :
# \nabla_h v_h + \frac{\alpha}{h} u_D \otimes n : v_h \otimes n \right), \\
# L_{c_h}(v_h) &= - \int_{\partial \Omega} u_D \cdot n \hat{u}_D \cdot v_h, \\
# \nabla_h v_h + \frac{\alpha}{h} u_D \otimes n : v_h \otimes n \right),
# \\
# L_{c_h}(v_h) &= - \int_{\partial \Omega} u_D \cdot n \hat{u}_D \cdot
# v_h, \\
# b_h(v, q) &= - \int_K \nabla \cdot v q.
# \end{align}
# $$
Expand Down
16 changes: 9 additions & 7 deletions python/demo/demo_pml.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
# # Electromagnetic scattering from a wire with perfectly matched layer condition
# # Electromagnetic scattering from a wire with perfectly matched layer condition # noqa
#
# Copyright (C) 2022 Michele Castriotta, Igor Baratta, Jørgen S. Dokken
#
Expand Down Expand Up @@ -221,7 +221,8 @@ def generate_mesh_wire(
# being hit normally by a TM-polarized electromagnetic wave.
#
# The formula are taken from:
# Milton Kerker, "The Scattering of Light and Other Electromagnetic Radiation",
# Milton Kerker, "The Scattering of Light and Other Electromagnetic
# Radiation",
# Chapter 6, Elsevier, 1969.
#
# ## Implementation
Expand Down Expand Up @@ -489,11 +490,12 @@ def pml_coordinates(x: ufl.indexed.Indexed, alpha: float, k0: complex, l_dom: fl
#
# $$
# \begin{align}
# \text{PML}_\text{corners} \rightarrow \mathbf{r}^\prime &= (x^\prime, y^\prime) \\
# \text{PML}_\text{corners} \rightarrow \mathbf{r}^\prime &= (x^\prime,
# y^\prime) \\
# \text{PML}_\text{rectangles along x} \rightarrow
# \mathbf{r}^\prime &= (x^\prime, y) \\
# \mathbf{r}^\prime &= (x^\prime, y) \\
# \text{PML}_\text{rectangles along y} \rightarrow
# \mathbf{r}^\prime &= (x, y^\prime).
# \mathbf{r}^\prime &= (x, y^\prime).
# \end{align}
# $$
#
Expand Down Expand Up @@ -662,8 +664,8 @@ def create_eps_mu(
# The final weak form in the PML region is:
#
# $$
# \int_{\Omega_{pml}}\left[\boldsymbol{\mu}^{-1}_{pml} \nabla \times \mathbf{E}
# \right]\cdot \nabla \times \bar{\mathbf{v}}-k_{0}^{2}
# \int_{\Omega_{pml}}\left[\boldsymbol{\mu}^{-1}_{pml} \nabla \times
# \mathbf{E} \right]\cdot \nabla \times \bar{\mathbf{v}}-k_{0}^{2}
# \left[\boldsymbol{\varepsilon}_{pml} \mathbf{E} \right]\cdot
# \bar{\mathbf{v}}~ d x=0,
# $$
Expand Down
3 changes: 2 additions & 1 deletion python/demo/demo_poisson.py
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,8 @@
# $$
# \begin{align}
# a(u, v) &:= \int_{\Omega} \nabla u \cdot \nabla v \, {\rm d} x, \\
# L(v) &:= \int_{\Omega} f v \, {\rm d} x + \int_{\Gamma_{N}} g v \, {\rm d} s.
# L(v) &:= \int_{\Omega} f v \, {\rm d} x + \int_{\Gamma_{N}} g v \,
# {\rm d} s.
# \end{align}
# $$
#
Expand Down
16 changes: 9 additions & 7 deletions python/demo/demo_poisson_matrix_free.py
Original file line number Diff line number Diff line change
Expand Up @@ -107,9 +107,10 @@
#
# Next, we locate the mesh facets that lie on the domain boundary
# $\partial\Omega$. We do this by first calling
# {py:func}`create_connectivity <dolfinx.mesh.topology.create_connectivity>`
# and then retrieving all facets on the boundary using
# {py:func}`exterior_facet_indices <dolfinx.mesh.exterior_facet_indices>`.
# {py:func}`create_connectivity
# <dolfinx.mesh.topology.create_connectivity>` and then retrieving all
# facets on the boundary using {py:func}`exterior_facet_indices
# <dolfinx.mesh.exterior_facet_indices>`.

tdim = mesh.topology.dim
mesh.topology.create_connectivity(tdim - 1, tdim)
Expand Down Expand Up @@ -156,10 +157,11 @@

# ### Matrix-free conjugate gradient solver
#
# The right hand side vector $b - A x_{\rm bc}$ is the assembly of the linear
# form $L$ where the essential Dirichlet boundary conditions are implemented
# using lifting. Since we want to avoid assembling the matrix `A`, we compute
# the necessary matrix-vector product using the linear form `M` explicitly.
# The right hand side vector $b - A x_{\rm bc}$ is the assembly of the
# linear form $L$ where the essential Dirichlet boundary conditions are
# implemented using lifting. Since we want to avoid assembling the matrix
# `A`, we compute the necessary matrix-vector product using the linear form
# `M` explicitly.

# Apply lifting: b <- b - A * x_bc
b = fem.assemble_vector(L_fem)
Expand Down
Loading
Loading