Start documenting demo

Author: jorgensd

_config.yml

@@ -1,7 +1,7 @@
 # Book settings
 # Learn more at https://jupyterbook.org/customize/config.html
 
-title: scifem
+title: fenicsx_ii
 author: Jørgen S. Dokken
 copyright: "2025"
 only_build_toc_files: true

_toc.yml

@@ -5,6 +5,10 @@ format: jb-book
 root: README
 parts:
 
+  - caption: Examples
+    chapters:
+      - file: "demos/coupled_poisson_solver"
+
   - caption: Reference
     chapters:
       - file: "docs/api.rst"

demos/coupled_poisson_solver.py

@@ -1,36 +1,66 @@
 # # A coupled 3D-1D manufactured solution
 #
-# This example is based on the manufactured solution from {cite}`masri2024coupled3d1d`.
-
+# This example is based on the strong formulation and manufactured solution
+# from {cite}`3D1Dmasri-masri2024coupled3d1d`.
+# We consider a 3D domain $\Omega$ with a 1D inclusion $\Lambda$.
+# We assume that $\Lambda$ can be parameterised as a curve $\lambda(s)$, $s\in(0,L)$.
+# Next, we define a generalized cylinder $B_{\Lambda,R(s)}$ with centerline $\Lambda$ and radius $R$.
+# A cross section of the cylinder at $s$ is defined as the circle $\Theta_{R}(s)$ with radius
+# $R: [0, L] \to \mathbb{R}^+$, area $A(s)$ and perimeter $P(s)$.
+# The boundary of $B_\Lambda$ is denoted $\Gamma$.
+# For a function $u\in V(\Omega)$, we define the restriction operator
+# $\Pi_R(s): V(\Omega) \to L^2(\Lambda)$ as
+#
+# $$
+# \Pi_R(u)(s) = \frac{1}{|P(s)|} \int_{\partial \Theta_R(s)} u~\mathrm{d}\gamma.
+# $$
+#
+# ## Mesh generation
+# We start by generating the two domains that we will consider in this demo.
+# ```{admonition} Strict inclusion
+# :class: dropdown
+# Note that $\Gamma$ should be strictly included in $\Omega$.
+# Additionally, we assume that $B_\Lambda \subset \Omega$.
+# ```
+# We start by defining a cube that spans $[-0.5,-0.5,-0.5]\times[0.5,0.5,0.5]$
+
+# +
 from mpi4py import MPI
-
-import basix.ufl
 import dolfinx
 import numpy as np
+import basix.ufl
 import ufl
 
-from fenicsx_ii import Average, Circle, LinearProblem
-
-# We create a 3D mesh (a cube)
-M = 64
-N = M
-volume = dolfinx.mesh.create_box(
-    MPI.COMM_WORLD,
+M = 32 # Number of elements in each spatial direction in the box
+N = M # Number of elements in the line
+comm = MPI.COMM_WORLD
+omega = dolfinx.mesh.create_box(
+    comm,
     [[-0.5, -0.5, -0.5], [0.5, 0.5, 0.5]],
     [M, M, M],
     cell_type=dolfinx.mesh.CellType.tetrahedron,
 )
-volume.name = "volume"
-
-
-# We create a 1D mesh (a line), which is not embedded in the 3D mesh
+omega.name = "volume"
+# -
+
+# Next we create a line that spans $[0,0,-0.5]\times[0,0,0.5]$.
+# Note that the discretized lines does not have to conform with the
+# 3D domain.
+# ```{admonition} Line embedded in 3D
+# Note that the geometrical dimension of the line should be 3. Therefore,
+# the built in mesh generators {py:func}`dolfinx.mesh.create_interval`
+# and {py:func}`dolfinx.mesh.create_unit_interval` does not work as
+# they use geometrical dimension 1.
+# ```
+
+# +
 l_min = -0.5
 l_max = 0.5
-if MPI.COMM_WORLD.rank == 0:
+if comm.rank == 0:
     nodes = np.zeros((N, 3), dtype=np.float64)
     nodes[:, 0] = np.full(nodes.shape[0], 0)
     nodes[:, 1] = np.full(nodes.shape[0], 0)
-    nodes[:, 2] = np.linspace(l_min, l_max, nodes.shape[0])[::-1]
+    nodes[:, 2] = np.linspace(l_min, l_max, nodes.shape[0])
     connectivity = np.repeat(np.arange(nodes.shape[0]), 2)[1:-1].reshape(
         nodes.shape[0] - 1, 2
     )
@@ -40,22 +70,55 @@
 c_el = ufl.Mesh(
     basix.ufl.element("Lagrange", basix.CellType.interval, 1, shape=(nodes.shape[1],))
 )
-line_mesh = dolfinx.mesh.create_mesh(
-    MPI.COMM_WORLD,
+lmbda = dolfinx.mesh.create_mesh(
+    comm,
     x=nodes,
     cells=connectivity,
     e=c_el,
     partitioner=dolfinx.mesh.create_cell_partitioner(
         dolfinx.mesh.GhostMode.shared_facet
     ),
 )
+# -
+
+# ## Variational formulation
+#
+# We consider the following coupled 3D-1D problem, which is used in
+# {cite}`3D1Dmasri-masri2024coupled3d1d` and stems from
+# {cite}`3D1Dmasri-laurino20193d1d`:
+# Find $u\in V(\Omega)$, $p\in Q(\Lambda)$ such that
+#
+# $$
+# \begin{align}
+#   - \nabla \cdot (\nabla u) + \xi (\Pi_R(u) - p))\delta_\Gamma &= f && \text{in } \Omega, \\
+#   - d_s(A d_s p) + P \xi (p - \Pi(u)) &= A \hat f  &&  \text{in } \Lambda, \\
+#   u&=g &&\text{on } \partial\Omega,\\
+#   A d_s p &=0 && \text{at } s\in\{0, 1\}.
+# \end{align}
+# $$
+#
+# with the variational formulation
+#
+# $$
+# \begin{align*}
+#   \int_\Omega \nabla u \cdot \nabla v~\mathrm{d}x
+#   + \int_\Gamma P\xi (\Pi_R(u) - p)\Pi_R(v)~\mathrm{d}s
+#   &= \int_\Omega f\cdot v~\mathrm{d}x\\
+#   \int_\Gamma d_s p \cdot d_s q~\mathrm{d}s
+#   + \int_\Gamma P\xi (p - \Pi_R(u))q~\mathrm{d}s
+#   &= \int_\Gamma A\hat f\cdot q~\mathrm{d}s
+# \end{align*}
+# $$
+
+from fenicsx_ii import Average, Circle, LinearProblem
+
 
 
 # We define the appropriate function spaces on each mesh
 
 degree = 1
-V = dolfinx.fem.functionspace(volume, ("Lagrange", degree))
-Q = dolfinx.fem.functionspace(line_mesh, ("Lagrange", degree))
+V = dolfinx.fem.functionspace(omega, ("Lagrange", degree))
+Q = dolfinx.fem.functionspace(lmbda, ("Lagrange", degree))
 
 # We define a mixed function space, to automate block extraction of the system
 
@@ -69,21 +132,21 @@
 
 R = 0.05
 q_degree = 20
-restriction_trial = Circle(line_mesh, R, degree=q_degree)
-restriction_test = Circle(line_mesh, R, degree=q_degree)
+restriction_trial = Circle(lmbda, R, degree=q_degree)
+restriction_test = Circle(lmbda, R, degree=q_degree)
 
 # Next, we define the integration measures on each mesh
 
-dx_3D = ufl.Measure("dx", domain=volume)
-dx_1D = ufl.Measure("dx", domain=line_mesh)
+dx_3D = ufl.Measure("dx", domain=omega)
+dx_1D = ufl.Measure("dx", domain=lmbda)
 
 # We can define the intermediate space that we interpolate the
 # 3D arguments into
 
 q_el = basix.ufl.quadrature_element(
-    line_mesh.basix_cell(), value_shape=(), degree=q_degree
+    lmbda.basix_cell(), value_shape=(), degree=q_degree
 )
-Rs = dolfinx.fem.functionspace(line_mesh, q_el)
+Rs = dolfinx.fem.functionspace(lmbda, q_el)
 
 
 # Next, we define the averaging operator of each of the test and trial functions
@@ -93,7 +156,7 @@
 
 # We define the various constants used in the variational formulation
 
-Alpha1 = dolfinx.fem.Constant(volume, 1.0)
+Alpha1 = dolfinx.fem.Constant(omega, 1.0)
 A = ufl.pi * R**2
 
 # and the variational form itself
@@ -119,38 +182,38 @@ def x(mesh: dolfinx.mesh.Mesh):
 P = 2 * ufl.pi * R
 
 a = Alpha1 * ufl.inner(ufl.grad(u), ufl.grad(v)) * dx_3D
-a += P * xi(line_mesh) * ufl.inner(avg_u - p, avg_v) * dx_1D
+a += P * xi(lmbda) * ufl.inner(avg_u - p, avg_v) * dx_1D
 a += A * ufl.inner(ufl.grad(p), ufl.grad(q)) * dx_1D
-a += P * xi(line_mesh) * ufl.inner(p - avg_u, q) * dx_1D
+a += P * xi(lmbda) * ufl.inner(p - avg_u, q) * dx_1D
 
-vol_rate = xi(volume) / (xi(volume) + 1)
+vol_rate = xi(omega) / (xi(omega) + 1)
 
-u_inside = vol_rate * u_line(x(volume))
-r = ufl.sqrt(x(volume)[0] ** 2 + x(volume)[1] ** 2)
-u_outside = vol_rate * (1 - R * ufl.ln(r / R)) * u_line(x(volume))
+u_inside = vol_rate * u_line(x(omega))
+r = ufl.sqrt(x(omega)[0] ** 2 + x(omega)[1] ** 2)
+u_outside = vol_rate * (1 - R * ufl.ln(r / R)) * u_line(x(omega))
 u_ex = ufl.conditional(r < R, u_inside, u_outside)
 
-p_ex = ufl.sin(ufl.pi * x(line_mesh)[2]) + 2
+p_ex = ufl.sin(ufl.pi * x(lmbda)[2]) + 2
 
-f_vol_in = vol_rate * ufl.pi**2 * ufl.sin(ufl.pi * x(volume)[2])
+f_vol_in = vol_rate * ufl.pi**2 * ufl.sin(ufl.pi * x(omega)[2])
 f_vol_out = (
-    vol_rate * (1 - R * ufl.ln(r / R)) * ufl.pi**2 * ufl.sin(ufl.pi * x(volume)[2])
+    vol_rate * (1 - R * ufl.ln(r / R)) * ufl.pi**2 * ufl.sin(ufl.pi * x(omega)[2])
 )
 f_vol = ufl.conditional(r < R, f_vol_in, f_vol_out)
 # - ufl.div(ufl.grad(u_ex))
 
-line_rate = xi(line_mesh) / (xi(line_mesh) + 1)
+line_rate = xi(lmbda) / (xi(lmbda) + 1)
 # f_line = - A*ufl.div(ufl.grad(p_ex)) + P * line_rate * p_ex
-f_line = A * ufl.sin(ufl.pi * x(line_mesh)[2]) * ufl.pi**2 + P * line_rate * p_ex
+f_line = A * ufl.sin(ufl.pi * x(lmbda)[2]) * ufl.pi**2 + P * line_rate * p_ex
 
 L = f_vol * v * dx_3D
 L += f_line * q * dx_1D
 
 # Exerior boundary conditions
-volume.topology.create_connectivity(volume.topology.dim - 1, volume.topology.dim)
-exterior_facets = dolfinx.mesh.exterior_facet_indices(volume.topology)
+omega.topology.create_connectivity(omega.topology.dim - 1, omega.topology.dim)
+exterior_facets = dolfinx.mesh.exterior_facet_indices(omega.topology)
 exterior_dofs = dolfinx.fem.locate_dofs_topological(
-    V, volume.topology.dim - 1, exterior_facets
+    V, omega.topology.dim - 1, exterior_facets
 )
 bc_expr = dolfinx.fem.Expression(u_ex, V.element.interpolation_points)
 u_bc = dolfinx.fem.Function(V)
@@ -180,12 +243,12 @@ def x(mesh: dolfinx.mesh.Mesh):
 # We move the solution to a Function and save to file
 
 
-with dolfinx.io.VTXWriter(volume.comm, "u_3D_solver.bp", [uh]) as vtx:
+with dolfinx.io.VTXWriter(omega.comm, "u_3D_solver.bp", [uh]) as vtx:
     vtx.write(0.0)
-with dolfinx.io.VTXWriter(line_mesh.comm, "p_1D_solver.bp", [ph]) as vtx:
+with dolfinx.io.VTXWriter(lmbda.comm, "p_1D_solver.bp", [ph]) as vtx:
     vtx.write(0.0)
 
-with dolfinx.io.VTXWriter(volume.comm, "u_bc.bp", [u_bc]) as vtx:
+with dolfinx.io.VTXWriter(omega.comm, "u_bc.bp", [u_bc]) as vtx:
     vtx.write(0.0)
 
 
@@ -217,24 +280,24 @@ def H1(f: ufl.core.expr.Expr, dx: ufl.Measure):
 H1_error_u = H1(uh - u_ex, dx_3D)
 H1_error_p = H1(ph - p_ex, dx_1D)
 
-h_vol = volume.comm.allreduce(
+h_vol = omega.comm.allreduce(
     np.max(
-        volume.h(
-            volume.topology.dim,
+        omega.h(
+            omega.topology.dim,
             np.arange(
-                volume.topology.index_map(volume.topology.dim).size_local,
+                omega.topology.index_map(omega.topology.dim).size_local,
                 dtype=np.int32,
             ),
         )
     ),
     op=MPI.MAX,
 )
-h_line = line_mesh.comm.allreduce(
+h_line = lmbda.comm.allreduce(
     np.max(
-        line_mesh.h(
-            line_mesh.topology.dim,
+        lmbda.h(
+            lmbda.topology.dim,
             np.arange(
-                line_mesh.topology.index_map(line_mesh.topology.dim).size_local,
+                lmbda.topology.index_map(lmbda.topology.dim).size_local,
                 dtype=np.int32,
             ),
         )
@@ -248,3 +311,12 @@ def H1(f: ufl.core.expr.Expr, dx: ufl.Measure):
     print(f"L2-error p: {L2_error_p:.6e}")
     print(f"H1-error u: {H1_error_u:.6e}")
     print(f"H1-error p: {H1_error_p:.6e}")
+
+
+## References
+
+# ```{bibliography}
+# :filter: cited
+# :labelprefix:
+# :keyprefix: 3D1Dmasri-
+# ```

references.bib

@@ -12,7 +12,7 @@ @article{masri2024coupled3d1d
 
 
 @article{dangelo20083d1d,
-author = {D'ANGELO, C. and QUARTERONI, A.},
+author = {D'Angelo, C. and Quarteroni, A.},
 title = {ON THE COUPLING OF 1D AND 3D DIFFUSION-REACTION EQUATIONS: APPLICATION TO TISSUE PERFUSION PROBLEMS},
 journal = {Mathematical Models and Methods in Applied Sciences},
 volume = {18},
@@ -39,5 +39,16 @@ @InProceedings{kutcha2021trace
 }
 
 
+@article{laurino20193d1d,
+	author = {Laurino, Federica and Zunino, Paolo},
+	title = {Derivation and analysis of coupled {PDEs} on manifolds with high dimensionality gap arising from topological model reduction},
+	DOI= "10.1051/m2an/2019042",
+	journal = {ESAIM: M2AN},
+	year = 2019,
+	volume = 53,
+	number = 6,
+	pages = "2047-2080",
+}
+