Re-factor shape optimization demos

demos/documented/cashocs_as_solver/control_solver/demo_control_solver.py

@@ -8,6 +8,10 @@
 #       jupytext_version: 1.14.4
 # ---
 
+# ```{eval-rst}
+# .. include:: ../../../global.rst
+# ```
+#
 # (demo_control_solver)=
 # # cashocs as Solver for Optimal Control Problems
 #
@@ -123,8 +127,8 @@
 # For both objects, one has to define them as a single UFL form for cashocs, as with the
 # state system and cost functional. In particular, the adjoint weak form has to be in
 # the form of a nonlinear variational problem, so that
-# `fenics.solve(adjoint_form == 0, p, adjoint_bcs)` could be used to solve it. In
-# particular, both forms have to include {py:class}`fenics.TestFunction` objects from
+# {python}`fenics.solve(adjoint_form == 0, p, adjoint_bcs)` could be used to solve it.
+# In particular, both forms have to include {py:class}`fenics.TestFunction` objects from
 # the control space and adjoint space, respectively, and must not contain
 # {py:class}`fenics.TrialFunction` objects.
 #

demos/documented/cashocs_as_solver/shape_solver/demo_shape_solver.py

@@ -8,6 +8,10 @@
 #       jupytext_version: 1.14.4
 # ---
 
+# ```{eval-rst}
+# .. include:: ../../../global.rst
+# ```
+#
 # (demo_shape_solver)=
 # # cashocs as Solver for Shape Optimization Problems
 #
@@ -132,7 +136,7 @@ def eps(u):
 # at its place.
 #
 # ::::{hint}
-# Alternatively, one could define the variable `vector_field` as follows
+# Alternatively, one could define the variable {python}`vector_field` as follows
 # :::python
 # space = VectorFunctionSpace(mesh, 'CG', 1)
 # vector_field = TestFunction(space)

demos/documented/misc/xdmf_io/demo_xdmf_io.py

@@ -1,20 +1,16 @@
-# Copyright (C) 2020-2022 Sebastian Blauth
-#
-# This file is part of cashocs.
-#
-# cashocs is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# cashocs is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with cashocs.  If not, see <https://www.gnu.org/licenses/>.
-
-"""For the documentation of this demo see https://cashocs.readthedocs.io/en/latest/demos/shape_optimization/doc_shape_poisson.html.
+# ---
+# jupyter:
+#   jupytext:
+#     text_representation:
+#       extension: .py
+#       format_name: light
+#       format_version: '1.5'
+#       jupytext_version: 1.14.4
+# ---
 
-"""
+# ```{eval-rst}
+# .. include:: ../../../global.rst
+# ```
+#
+# (demo_xdmf_io)=
+# # Writing and Reading XDMF Files

demos/documented/optimal_control/box_constraints/demo_box_constraints.py

@@ -8,6 +8,10 @@
 #       jupytext_version: 1.14.4
 # ---
 
+# ```{eval-rst}
+# .. include:: ../../../global.rst
+# ```
+#
 # (demo_box_constraints)=
 # # Control Constraints
 #
@@ -88,8 +92,8 @@
 u_b = interpolate(Expression("50*x[0]", degree=1), V)
 
 # which just corresponds to two functions, generated from {py:class}`fenics.Expression`
-# objects via {py:func}`fenics.interpolate`. These are then put into the list `cc`,
-# which models the control constraints, i.e.,
+# objects via {py:func}`fenics.interpolate`. These are then put into the list
+# {python}`cc`, which models the control constraints, i.e.,
 
 cc = [u_a, u_b]
 
@@ -105,18 +109,19 @@
 # cc = [u_a, u_b]
 # :::
 #
-# and completely analogous with `float('-inf')` for no constraint on the lower bound.
-# Moreover, note that the specification of using either constant `float` values and
-# {py:class}`fenics.Function` objects can be mixed arbitrarily, so that one can, e.g.,
-# specify a constant value for the upper boundary and use a {py:class}`fenics.Function`
-# on the lower one.
+# and completely analogous with {python}`float('-inf')` for no constraint on the lower
+# bound. Moreover, note that the specification of using either constant {python}`float`
+# values and {py:class}`fenics.Function` objects can be mixed arbitrarily, so that one
+# can, e.g., specify a constant value for the upper boundary and use a
+# {py:class}`fenics.Function` on the lower one.
 # ::::
 #
 # ### Setup of the optimization problem and its solution
 #
 # Now, we can set up the optimal control problem as we did before, using the additional
-# keyword argument `control_constraints` into which we put the list `cc`, and then solve
-# it via the {py:meth}`ocp.solve() <cashocs.OptimalControlProblem.solve>` method
+# keyword argument {python}`control_constraints` into which we put the list
+# {python}`cc`, and then solve it via the {py:meth}`ocp.solve()
+# <cashocs.OptimalControlProblem.solve>` method
 
 ocp = cashocs.OptimalControlProblem(
     e, bcs, J, y, u, p, config=config, control_constraints=cc

demos/documented/optimal_control/constraints/demo_constraints.py

@@ -8,6 +8,10 @@
 #       jupytext_version: 1.14.4
 # ---
 
+# ```{eval-rst}
+# .. include:: ../../../global.rst
+# ```
+#
 # (demo_constraints)=
 # # Treatment of additional constraints
 #
@@ -103,10 +107,10 @@
 # indicator function for the lower left quarter multiplied by the state variable y.
 # The next argument is the right-hand side of the equality constraint, i.e.,
 # $c_{b,l}$ which we choose as 0 in this example. Finally, the keyword argument
-# `measure` is used to specify the integration measure that should be used to define
-# where the constraint is given. Typical examples are a volume measure (`dx`, as it is
-# the case here) or surface measure (`ds`, which could be used if we wanted to pose the
-# constraint only on the boundary).
+# {python}`measure` is used to specify the integration measure that should be used to
+# define where the constraint is given. Typical examples are a volume measure (`dx`, as
+# it is the case here) or surface measure (`ds`, which could be used if we wanted to
+# pose the constraint only on the boundary).
 #
 # Let's move on to the next constraint. Again, we have an equality constraint, but now
 # it is a scalar value which is constrained, and its given by the integral over some
@@ -122,7 +126,7 @@
 # case $y^2$ multiplied by the indicator function of the bottom right quarter, i.e.,
 # the left-hand side of the constraint. The second and final argument for this
 # constraint is right-hand side of the constraint, i.e., $c_{b,r}$, which we choose as
-# `0.01` in this example.
+# {python}`0.01` in this example.
 #
 # Let's move on to the interesting case of inequality constraints. Let us first consider
 # a setting similar to before, where the constraint's left-hand side is given by an
@@ -135,9 +139,9 @@
 
 # Here, as before, the first argument is the left-hand side of the constraint, i.e., the
 # UFL form of the integrand, in this case $y$ times the indicator function of the top
-# left quarter, which is to be integrated over the measure `dx`. The second argument
-# `lower_bound = -0.025` specifies the lower bound for this inequality constraint, that
-# means, that $c_{t,l} = -0.025$ in our case.
+# left quarter, which is to be integrated over the measure {python}`dx`. The second
+# argument {python}`lower_bound = -0.025` specifies the lower bound for this inequality
+# constraint, that means, that $c_{t,l} = -0.025$ in our case.
 #
 # Finally, let us take a look at the case of pointwise inequality constraint. This is,
 # as before, implemented via the {py:class}`cashocs.InequalityConstraint` class
@@ -148,14 +152,16 @@
 
 # Here, again the first argument is the function on the left-hand side of the
 # constraint, i.e., $y$ times the indicator function of the top right quarter. The
-# second argument, `upper_bound=0.25`, defines the right-hand side of the constraint,
-# i.e., we choose $c_{t,r} = 0.25$. Finally, as for the pointwise equality constraint,
-# we specify the integration measure for which the constraint is posed, in our case
-# `measure=dx`, as we consider the constraint pointwise in the domain $\Omega$.
+# second argument, {python}`upper_bound=0.25`, defines the right-hand side of the
+# constraint, i.e., we choose $c_{t,r} = 0.25$. Finally, as for the pointwise equality
+# constraint, we specify the integration measure for which the constraint is posed, in
+# our case {python}`measure=dx`, as we consider the constraint pointwise in the domain
+# $\Omega$.
 #
 # :::{note}
-# For bilateral inequality constraints we can use both keyword arguments `upper_bound`
-# and `lower_bound` to define both bounds for the constraint.
+# For bilateral inequality constraints we can use both keyword arguments
+# {python}`upper_bound` and {python}`lower_bound` to define both bounds for the
+# constraint.
 # :::
 #
 # As is usual in cashocs, once we have defined multiple constraints, we gather them into
@@ -182,7 +188,7 @@
 # :::{note}
 # To be able to treat (nearly) arbitrary types of constraints, cashocs regularizes
 # these using either an augmented Lagrangian method or a quadratic penalty method.
-# Which method is used can be specified via the keyword argument `method`, which
+# Which method is used can be specified via the keyword argument {python}`method`, which
 # is chosen to be an augmented Lagrangian method (`'AL'`) in this demo.
 # :::
 #

demos/documented/optimal_control/control_boundary_conditions/demo_control_boundary_conditions.py

@@ -8,6 +8,10 @@
 #       jupytext_version: 1.14.4
 # ---
 
+# ```{eval-rst}
+# .. include:: ../../../global.rst
+# ```
+#
 # (demo_control_boundary_conditions)=
 # # Boundary conditions for control variables
 #
@@ -102,8 +106,9 @@
 ocp.solve()
 
 # where the only additional parts in comparison to {ref}`demo_poisson` are the keyword
-# arguments `riesz_scalar_products`, which was already covered in
-# {ref}`demo_neumann_control`, and `control_bcs_list`, which we have defined previously.
+# arguments {python}`riesz_scalar_products`, which was already covered in
+# {ref}`demo_neumann_control`, and {python}`control_bcs_list`, which we have defined
+# previously.
 #
 # After solving this problem with cashocs, we visualize the solution with the code
 

demos/documented/optimal_control/dirichlet_control/demo_dirichlet_control.py

@@ -8,6 +8,10 @@
 #       jupytext_version: 1.14.4
 # ---
 
+# ```{eval-rst}
+# .. include:: ../../../global.rst
+# ```
+#
 # (demo_dirichlet_control)=
 # # Dirichlet Boundary Control
 #
@@ -127,9 +131,9 @@
 u = Function(V)
 # -
 
-# The only difference is, that we now also define `n`, which is the outer unit
-# normal vector on $\Gamma$, and `h`, which is the maximum length of an edge of the
-# respective finite element (during assemly).
+# The only difference is, that we now also define {python}`n`, which is the outer unit
+# normal vector on $\Gamma$, and {python}`h`, which is the maximum length of an edge of
+# the respective finite element (during assemly).
 #
 # Then, we define the Dirichlet boundary conditions, which are enforced strongly.
 # As we use Nitsche's method to implement the boundary conditions on the entire
@@ -195,7 +199,7 @@
 [bc.apply(bdry_idx.vector()) for bc in bcs]
 mask = np.where(bdry_idx.vector()[:] == 1)[0]
 
-# Then, we restrict both `y` and `u` to the boundary by
+# Then, we restrict both {python}`y` and {python}`u` to the boundary by
 
 y_bdry = Function(V)
 u_bdry = Function(V)
@@ -217,7 +221,7 @@
 print("Error L^\infty: " + format(error_inf, ".3e") + " %")
 print("Error L^2: " + format(error_l2, ".3e") + " %")
 
-# We see, that with `eta = 1e4` we get a relative error of under 5e-3 % in the
+# We see, that with {python}`eta = 1e4` we get a relative error of under 5e-3 % in the
 # $L^\infty(\Omega)$ norm, and under 5e-4 in the $L^2(\Omega)$ norm, which is
 # sufficient for applications.
 #

demos/documented/optimal_control/heat_equation/demo_heat_equation.py

@@ -8,9 +8,13 @@
 #       jupytext_version: 1.14.4
 # ---
 
+# ```{eval-rst}
+# .. include:: ../../../global.rst
+# ```
+#
 # (demo_heat_equation)=
 # # Distributed Control for Time Dependent Problems
-
+#
 # ## Problem Formulation
 #
 # In this demo  we take a look at how time dependent problems can be treated with
@@ -117,7 +121,7 @@
 t_end = 1.0
 t_array = np.linspace(t_start, t_end, int(1 / dt))
 
-# Here, `t_array` is a numpy array containing all time steps. Note, that we do **not**
+# Here, {python}`t_array` is a numpy array containing all time steps. Note, that we do **not**
 # include $t=0$ in the array. This is due to the fact, that the initial condition
 # is prescribed and fixed.
 #
@@ -130,8 +134,8 @@
 controls = [Function(V) for i in range(len(t_array))]
 adjoints = [Function(V) for i in range(len(t_array))]
 
-# Note, that `states[k]` corresponds to $y_{k+1}$ since indices start at
-# `0` in most programming languages (as it is the case in python).
+# Note, that {python}`states[k]` corresponds to $y_{k+1}$ since indices start at
+# {python}`0` in most programming languages (as it is the case in python).
 #
 # As the boundary conditions are not time dependent, we can initialize them now, and
 # repeat them in a list, since they are the same for every state
@@ -193,14 +197,14 @@
     )
 
 # ::::{admonition} Description of the for-loop
-# At the beginning, the 'current' time t is determined from ``t_array``, and the
+# At the beginning, the 'current' time t is determined from {python}`t_array`, and the
 # expression for the desired state is updated to reflect the current time with the code
 # :::python
 # t = t_array[k]
 # y_d_expr.t = t
 # :::
 #
-# The following line sets the object `y` to $y_{k+1}$
+# The following line sets the object {python}`y` to $y_{k+1}$
 # :::python
 # y = states[k]
 # :::
@@ -215,7 +219,7 @@
 # :::
 #
 # which ensures that $y_0 = 0$, which corresponds to the initial condition
-# $y^{(0)} = 0$. Hence, `y_prev` indeed corresponds to $y_{k}$.
+# $y^{(0)} = 0$. Hence, {python}`y_prev` indeed corresponds to $y_{k}$.
 #
 # Moreover, we get the current control and adjoint state via
 # :::python

demos/documented/optimal_control/iterative_solvers/demo_iterative_solvers.py

@@ -8,9 +8,13 @@
 #       jupytext_version: 1.14.4
 # ---
 
+# ```{eval-rst}
+# .. include:: ../../../global.rst
+# ```
+#
 # (demo_iterative_solvers)=
 # # Iterative Solvers for State and Adjoint Systems
-
+#
 # ## Problem Formulation
 #
 # cashocs is also capable of using iterative solvers through the linear algebra
@@ -187,13 +191,13 @@
 # -
 
 # :::{note}
-# Note, that if the `ksp_options` and `adjoint_ksp_options` are not passed
-# to the {py:class}`OptimalControlProblem <cashocs.OptimalControlProblem>` or `None`,
-# which is the default value of these keyword parameters, then the direct solver MUMPS
-# is used. Moreover, if one wants to use identical options for state and adjoint
-# systems, then only the `ksp_options` have to be passed. This is because
-# `adjoint_ksp_options` always mirrors the ksp_options in case that the input is `None`
-# for `adjoint_ksp_options`.
+# Note, that if the {python}`ksp_options` and {python}`adjoint_ksp_options` are not
+# passed to the {py:class}`OptimalControlProblem <cashocs.OptimalControlProblem>` or
+# {python}`None`, which is the default value of these keyword parameters, then the
+# direct solver MUMPS is used. Moreover, if one wants to use identical options for state
+# and adjoint systems, then only the {python}`ksp_options` have to be passed. This is
+# because {python}`adjoint_ksp_options` always mirrors the ksp_options in case that the
+# input is {python}`None` for {python}`adjoint_ksp_options`.
 # :::
 #
 # We visualize the results of the optimization with the lines