Documentation of the Config Files for Shape Optimization Problems#

Let us take a detailed look at the config files for shape optimization problems and discusss the corresponding parameters. The corresponding config file used for this discussion is config.ini, which is the config file used for Shape Optimization with a Poisson Problem.

For shape optimization problems, the config file is a lot larger compared to the config files for optimal control. However, several important parameters are shared between both _typing of optimization problems.

A general config file for shape optimization has the following sections Mesh, StateSystem, OptimizationRoutine, AlgoLBFGS, AlgoCG, ShapeGradient, Regularization, MeshQuality and Output. We go over these sections and each parameter in them in the following.

As in Documentation of the Config Files for Optimal Control Problems, we refer to the documentation of the configparser module for a detailed description of how these config files can be structured. Moreover, we remark that cashocs has a default behavior for almost all of these parameters, which is triggered when they are NOT specified in the config file, and we will discuss this behavior for each parameter in this tutorial. For a summary over all parameters and their default values look at the end of this page.

Section Mesh#

The mesh section is, in contrast to the corresponding one for optimal control problems, more important and also has more parameters. However, its primary importance is for remeshing, which we cover in Remeshing with cashocs.

As first parameter, we have

mesh_file = ./mesh/mesh.xdmf

This specifies a path to a .xdmf file containing the discretized geometry. For all purposes, cashocs assumes that this .xdmf file was generated via conversion from a GMSH file using the command line command cashocs.convert().

Note, that the corresponding files for the boundaries and subdomains are generated automatically with cashocs.convert(), and they will also be read by import_mesh if they are present.

The second parameter in the Mesh section, gmsh_file, is defined via

gmsh_file = ./mesh/mesh.msh

This defines the path to the GMSH .msh file which was used to create the .xdmf file specified in mesh_file. As before, this parameter is only relevant for remeshing purposes, and not needed otherwise.

The next parameter is geo_file, which is the final file we need for remeshing ( and only there). It is also given by a path to a file, in this case to the GMSH .geo file used to generate the gmsh_file. It is specified, .e.g., as

geo_file = ./mesh/mesh.geo


For a detailed discussion of how to use these parameters we refer to Remeshing with cashocs.

Next up is a boolean flag that is used to indicate whether remeshing shall be performed

remesh = False

As the remeshing feature is experimental, we do advise to always try without remeshing. Note, that by default this flag is set to remehs = False so that remeshing is disabled.

Finally, we have the boolean flag show_gmsh_output, specified via

show_gmsh_output = False

This is used to toggle on / off the terminal output of GMSH when it performs a remeshing operation. This can be helpful for debugging purposes. By default, this is set to show_gmsh_output = False.

As stated throughout the Mesh section, these parameters are optional most of the time, and are only really required for remeshing. You can safely leave them out of your config file, and you should not need them, unless you want to perform remeshing.

Section StateSystem#

The StateSystem section is in complete analogy to the corresponding one for optimal control problems. For the sake of completeness, we briefly recall the parameters here, anyway.

The first parameter is is_linear, and can be set as

is_linear = True

This is a boolean flag that indicates whether the state system is linear or not. The default value for this parameter is is_linear = False, as every linear problem can also be interpreted as a nonlinear one.

The next parameters are used to define the tolerances of the Newton solver, in case a nonlinear state system has to be solved

newton_rtol = 1e-11
newton_atol = 1e-13

Here, newton_rtol sets the relative, and newton_atol the absolute tolerance for Newton’s method. Their default values are newton_rtol = 1e-11 and newton_atol = 1e-13.

The next parameter for the Newton iteration is the maximum number of iterations it is allowed to perform before the iteration is cancelled. This is controlled via

newton_iter = 50

which defaults to newton_iter = 50.

The parameter newton_damped, which is set via

newton_damped = True

is a boolean flag, indicating whether a damping strategy should be performed for the Newton method, or whether the classical Newton-Raphson iteration shall be used. This defaults to newton_damped = False (as this is faster), but for some problems it might be beneficial to use damping in order to enhance the convergence of the nonlinear solver.

Additionally, we have the boolean parameter newton_inexact, defined via

newton_inexact = False

which sets up an inexact Newton method for solving nonlinear problems in case this is newton_inexact = True. The default is newton_inexact = False.

Next, we have the parameter

newton_verbose = False

This is used to toggle the verbose output of the Newton method for the state system. By default this is set to newton_verbose = False so that there is not too much noise in the terminal.

The upcoming parameters are used to define the behavior of a Picard iteration, that may be used if we have multiple variables.


For a detailed discussion of how to use the Picard iteration to solve a coupled state system, we refer to Coupled Problems - Picard Iteration. Note, that this demo is written for optimal control problems, but the definition of the state system can be transferred analogously to shape optimization problems, too.

First, we have a boolean flag, set via

picard_iteration = False

which determines whether the Picard iteration is enabled or not. This defaults to picard_iteration = False, so that the Picard solver is disabled by default. The following two parameters determine, analogously to above, the tolerances for the Picard iteration

picard_rtol = 1e-10
picard_atol = 1e-12

The default values for these parameters are picard_rtol = 1e-10 and picard_atol = 1e-12. Moreover, note that the tolerances of the Newton solver are adjusted automatically in case a Picard iteration is performedm, so that an inexact Picard iteration is used.

The maximum amout of iterations for the Picard iteration are set with

picard_iter = 10

The default value for this is given by picard_iter = 50.

Finally, we can enable verbose output of the Picard iteration with the following boolean flag

picard_verbose = False

which is set to picard_verbose = False by default.

Section OptimizationRoutine#

The section OptimizationRoutine also closely resembles the one for optimal control problems. Again, we will take a brief look at all parameters here

The first parameter that can be controlled via the config file is algorithm, which is set via

algorithm = lbfgs

There are three possible choices for this parameter for shape optimization problems, namely

  • algorithm = gd or algorithm = gradient_descent : A gradient descent method

  • algorithm = cg, algorithm = conjugate_gradient, algorithm = ncg, algorithm = nonlinear_cg : Nonlinear CG methods

  • algorithm = lbfgs or algorithm = bfgs : limited memory BFGS method.

Thereafter, we specify the tolerances for the optimization algorithm with the parameters

rtol = 5e-3
atol = 0.0

Again, rtol denotes the relative, and atol the absolute tolerance, and the defaults for these parameters are given by rtol = 1e-3, and atol = 0.0.

The next parameter is used to control the maximum number of iterations performed by the optimization algorithm. It is set via

max_iter = 50

and defaults to max_iter = 100.

Next up, we have the initial guess for the step size, which can be determined via

initial_stepsize = 1.0

The default behavior is given by initial_stepsize = 1.0.

The next parameter is given by

safeguard_stepsize = True

This parameter can be used to activate safeguarding of the initial stepsize for line search methods. This helps to choose an apropriate stepsize for the initial iteration even if the problem is poorly scaled.

The upcoming parameters are used for the Armijo rule

epsilon_armijo = 1e-4
beta_armijo = 2

They are used to verify that the condition

\[J((I + t \mathcal{V})\Omega) \leq J(\Omega) + \varepsilon_{\text{Armijo}}\ t\ dJ(\Omega)[\mathcal{V}]\]

holds, and if this is not satisfied, the stepsize is updated via \(t = \frac{t}{\beta_{\text{Armijo}}}\). As default values for these parameters we use epsilon_armijo = 1e-4 as well as beta_armijo = 2.

Next, we have a set of two parameters which detail the methods used for computing gradients in cashocs. These parameters are

gradient_method = direct

as well as

gradient_tol = 1e-9

The first parameter, gradient_method can be either gradient_method = direct or gradient_method = iterative. In the former case, a direct solver is used to compute the gradient (using a Riesz projection) and in the latter case, an iterative solver is used to do so. In case we have gradient_method = iterative, the parameter gradient_tol is used to specify the (relative) tolerance for the iterative solver, in the other case the parameter is not used.

The following parameter, soft_exit, is a boolean flag which determines how the optimization algorithm is terminated in case it does not converge. If soft_exit = True, then an error message is printed, but code after the solve call of the optimization problem will still be executed. However, when soft_exit = False, cashocs raises an exception and terminates. This is set via

soft_exit = False

and is set to soft_exit = False by default.

Section LineSearch#

In this section, parameters regarding the line search can be specified. The type of the line search can be chosen via the parameter

method = armijo

Possible options are method = armijo, which performs a simple backtracking line search based on the armijo rule with fixed steps (think of halving the stepsize in each iteration), and method = polynomial, which uses polynomial models of the cost functional restricted to the line to generate “better” guesses for the stepsize. The default is method = armijo.

The next parameter, polynomial_model, specifies, which type of polynomials are used to generate new trial stepsizes. It is set via

polynomial_model = cubic

The parameter can either be polynomial_model = quadratic or polynomial_model = cubic. If this is polynomial_model = quadratic, a quadratic interpolation polynomial along the search direction is generated and this is minimized analytically to generate a new trial stepsize. Here, only the current function value, the direction derivative of the cost functional in direction of the search direction, and the most recent trial stepsize are used to generate the polynomial. In case that polynomial_model = cubic, the last two trial stepsizes (when available) are used in addition to the current cost functional value and the directional derivative, to generate a cubic model of the one-dimensional cost functional, which is then minimized to compute a new trial stepsize.

For the polynomial models, we also have a safeguarding procedure, which ensures that trial stepsizes cannot be chosen too large or too small, and which can be configured with the following two parameters. The trial stepsizes generate by the polynomial models are projected to the interval \([\beta_{low} \alpha, \beta_{high} \alpha]\), where \(\alpha\) is the previous trial stepsize and \(\beta_{low}, \beta_{high}\) are factors which can be set via the parameters factor_low and factor_high. In the config file, this can look like this

factor_high = 0.5
factor_low = 0.1

and the values specified here are also the default values for these parameters.

Finally, we have the parameter

fail_if_not_converged = False

which determines, whether the line search is terminated if the state system cannot be solved at the current iterate. If this is fail_if_not_converged = True, then an exception is raised. Otherwise, the iterate is counted as having too high of a function value and the stepsize is “halved” and a new iterate is formed.

Section AlgoLBFGS#

Next, we discuss the parameters relevant for the limited memory BFGS method. For details regarding this method, we refer to Schulz, Siebenborn, and Welker, Efficient PDE Constrained Shape Optimization Based on Steklov-Poincaré-Type Metrics, where the methods are introduced.

The first parameter, bfgs_memory_size, determines how large the storage of the BFGS method is. It is set via

bfgs_memory_size = 3

Usually, a higher storage leads to a better Hessian approximation, and thus to faster convergence. However, this also leads to an increased memory usage. Typically, values below 5 already work very well. The default is bfgs_memory_size = 5.

The other parameter for the BFGS method is

use_bfgs_scaling = True

This determines, whether one should use a scaling of the initial Hessian approximation (see Nocedal and Wright, Numerical Optimization). This is usually very beneficial and should be kept enabled (which is the default).

Third, we have the parameter bfgs_periodic_restart, which is set in the line

bfgs_periodic_restart = 0

This is a non-negative integer value, which indicates the number of BFGS iterations, before a reinitialization takes place. In case that this is bfgs_periodic_restart = 0 (which is the default), no restarts are performed.

Finally, we have the parameter damped, which can be set with

damped = False

This parameter is a boolean flag, which indicates whether Powell’s damping (on H) should be used or not. This is useful, when the curvature condition is not satisfied and (without damping) a restart would be required. The default is damped = False.

Section AlgoCG#

The following parameters are used to define the behavior of the nonlinear conjugate gradient methods for shape optimization. For more details on this, we refer to the preprint Blauth, Nonlinear Conjugate Gradient Methods for PDE Constrained Shape Optimization Based on Steklov-Poincaré-Type Metrics.

First, we define which nonlinear CG method is used by

cg_method = DY

Available options are

  • cg_method = FR : The Fletcher-Reeves method

  • cg_method = PR : The Polak-Ribiere method

  • cg_method = HS : The Hestenes-Stiefel method

  • cg_method = DY : The Dai-Yuan method

  • cg_method = HZ : The Hager-Zhang method

The default value is cg_method = FR. As for optimal control problems, the subsequent parameters are used to define the restart behavior of the nonlinear CG methods. First, we have

cg_periodic_restart = False

This boolean flag en- or disables that the NCG methods are restarted after a fixed amount of iterations, which is specified via

cg_periodic_its = 5

i.e., if cg_periodic_restart = True and cg_periodic_its = n, then the NCG method is restarted every \(n\) iterations. The default behavior is given by cg_periodic_restart = False and cg_periodic_its = 10.

Alternatively, there also exists a relative restart criterion (see Nocedal and Wright, Numerical Optimization), which can be enabled via the boolean flag cg_relative_restart, which is defined in the line

cg_relative_restart = False

and the corresponding restart tolerance is set in

cg_restart_tol = 0.5

Note, that cg_restart_tol should be in \((0, 1)\). If two subsequent gradients generated by the nonlinear CG method are not “sufficiently orthogonal”, the method is restarted with a gradient step. The default behavior is given by cg_relative_restart = False and cg_restart_tol = 0.25.

Section ShapeGradient#

After we have specified the behavior of the solution algorithm, this section is used to specify parameters relevant to the computation of the shape gradient. Note, that by shape gradient we refer to the following object.

Let \(\mathcal{S} \subset \{ \Omega \;\vert\; \Omega \subset \mathbb{R}^d \}\) be a subset of the power set of \(\mathbb{R}^d\). Let \(J\) be a shape differentiable functional \(J \colon \mathcal{S} \to \mathbb{R}\) with shape derivative \(dJ(\Omega)[\mathcal{V}]\). Moreover, let \(a \colon H \times H \to \mathbb{R}\) be a symmetric, continuous, and coercive bilinear form on the Hilbert space \(H\). Then, the shape gradient \(\mathcal{G}\) of \(J\) (w.r.t. \(a\)) is defined as the solution of the problem

\[\begin{split}\text{Find } \mathcal{G} \in H \text{ such that } \\ \quad a(\mathcal{G}, \mathcal{V}) = dJ(\Omega)[\mathcal{V}].\end{split}\]

For PDE constrained shape optimization, it is common to use a bilinear form based on the linear elasticity equations, which enables smooth mesh deformations. This bilinear form is given as follows, in a general form, that is also implemented in cashocs

\[a \colon H \times H; \quad a(\mathcal{W}, \mathcal{V}) = \int_\Omega 2 \mu \left( \varepsilon(\mathcal{W}) : \varepsilon(\mathcal{V}) \right) + \lambda \left( \text{div}(\mathcal{W}) \text{div}(\mathcal{V}) \right) + \delta \left( V \cdot W \right) \text{ d}x,\]

where \(H\) is some suitable subspace of \(H^1(\Omega)^d\) and \(\varepsilon(\mathcal{V}) = \frac{1}{2}(D\mathcal{V} + D\mathcal{V}^\top)\) is the symmetric part of the Jacobian. The subspace property is needed to include certain geometrical constraints of the shape optimization problem, which fix certain boundaries, into the shape gradient. For a detailed description of this setting we refer to the preprint Blauth, Nonlinear Conjugate Gradient Methods for PDE Constrained Shape Optimization Based on Steklov-Poincaré-Type Metrics. Moreover, we note that for the second Lamé parameter \(\mu\), cashocs implements an idea from Schulz and Siebenborn, Computational Comparison of Surface Metric for PDE Constrained Shape Optimization: There, it is proposed to compute \(\mu\) as the solution of the Laplace problem

\[\begin{split}\begin{alignedat}{2} - \Delta \mu &= 0 \quad &&\text{ in } \Omega, \\ \mu &= \mu_\text{def} \quad &&\text{ on } \Gamma^\text{def},\\ \mu &= \mu_\text{fix} \quad &&\text{ on } \Gamma^\text{fix}.\\ \end{alignedat}\end{split}\]

This allows to give the deformable and fixed boundaries a different stiffness, which is then smoothly extended into the interior of the domain. Moreover, they propose to use the solution of this Laplace equation directly for 2D problems, and to use \(\sqrt{\mu}\) for 3D problems.

Moreover, let us take a look at the possible _typing of boundaries that can be used with cashocs. In principle, there exist two _typing: deformable and fixed boundaries. On fixed boundaries, we impose homogeneous Dirichlet boundary conditions for the shape gradient, so that these are not moved under the corresponding deformation. In cashocs, we define what boundaries are fixed and deformable via their markers, which are either defined in the corresponding python script, or in the GMSH file, if such a mesh is imported.

The config file for Shape Optimization with a Poisson Problem defines the deformable boundaries with the command

shape_bdry_def = [1]


Remember, that in Shape Optimization with a Poisson Problem, we defined boundaries with the commands

boundary = CompiledSubDomain('on_boundary')
boundaries = MeshFunction('size_t', mesh, dim=1)
boundary.mark(boundaries, 1)

Hence, we see that the marker 1 corresponds to the entire boundary, so that this is set to being deformable through the config.

As we do not have a fixed boundary for this problem, the corresponding list for the fixed boundaries is empty

shape_bdry_fix = []

Note, that cashocs also gives you the possibility of defining partially constrainted boundaries, where only one axial component is fixed, whereas the other two are not. These are defined in

shape_bdry_fix_x = []
shape_bdry_fix_y = []
shape_bdry_fix_z = []

For these, we have that shape_bdry_fix_x is a list of all markers whose corresponding boundaries should not be deformable in x-direction, but can be deformed in the y- and z-directions. Of course you can constrain a boundary to be only variable in a single direction by adding the markers to the remaining lists.

The next parameter is specified via

use_pull_back = True

This parameter is used to determine, whether the material derivative should be computed for objects that are not state or adjoint variables. This is enabled by default.


This parameter should always be set to use_pull_back = True, otherwise the shape derivative might be wrong. Only disable it when you are sure what you are doing.

Furthermore, note that the material derivative computation is only correct, as long as no differential operators act on objects that are not state or adjoint variables. However, this should only be a minor restriction and not relevant for almost all problems.


See Inverse Problem in Electric Impedance Tomography for a case, where we use use_pull_back = False.

The next parameters determine the coefficients of the bilinear form \(a\). First, we have the first Lamé parameter \(\lambda\), which is set via

lambda_lame = 1.428571428571429

The default value for this is lambda_lame = 0.0.

Next, we specify the damping parameter \(\delta\) with the line

damping_factor = 0.2

The default for this is damping_factor = 0.0.


As the default value for the damping factor is damping_factor = 0.0, this should be set to a positive value in case the entire boundary of a problem is deformable. Otherwise, the Riesz identification problem for the shape gradient is not well-posed.

Finally, we define the values for \(\mu_\text{def}\) and \(\mu_\text{fix}\) via

mu_fix = 0.35714285714285715
mu_def = 0.35714285714285715

The default behavior is given by mu_fix = 1.0 and mu_def = 1.0.

The parameter use_sqrt_mu is a boolean flag, which switches between using \(\mu\) and \(\sqrt{\mu}\) as the stiffness for the linear elasticity equations, as discussed above. This is set via

use_sqrt_mu = False

and the default value is use_sqrt_mu = False.

The next line in the config file is

inhomogeneous = False

This determines, whether an inhomogeneous linear elasticity equation is used to project the shape gradient. This scales the parameters \(\mu, \lambda\) and \(\delta\) by \(\frac{1}{\text{vol}}\), where \(\text{vol}\) is the volume of the current element (during assembly). This means, that smaller elements get a higher stiffness, so that the deformation takes place in the larger elements, which can handle larger deformations without reducing their quality too much. For more details on this approach, we refer to the paper Blauth, Leithäuser, and Pinnau, Model Hierarchy for the Shape Optimization of a Microchannel Cooling System.

Moreover, the parameter

update_inhomogeneous = False

can be used to update the local mesh size after each mesh deformation, in case this is update_inhomogeneous = True, so that elements which become smaller also obtain a higher stiffness and vice versa. The default is update_inhomogeneous = False.

For the inhomogeneous mesh stiffness, we also have the parameter inhomogeneous_exponent, which is specified via

inhomogeneous_exponent = 1.0

This parameter can be used to specify an exponent for the inhomogeneous mesh stiffness, so that the parameters \(\mu, \lambda\) and \(\delta\) are scaled by \(\left( \frac{1}{\text{vol}} \right)^p\), where \(p\) is specified in inhomogeneous_exponent. The default for this parameter is inhomogeneous_exponent = 1.0.

There is also a different possibility to define the stiffness parameter \(\mu\) using cashocs, namely to define \(\mu\) in terms of how close a point of the computational domain is to a boundary. In the following we will explain this alternative way of defining \(\mu\). To do so, we must first set the boolean parameter

use_distance_mu = True

which enables this formulation and deactivates the previous one. Note that by default, the value is use_distance_mu = False. Next, we have the parameters dist_min, dist_max, mu_min and mu_max. These do the following: If the distance to the boundary is smaller than dist_min, the value of \(\mu\) is set to mu_min, and if the distance to the boundary is larger than dist_max, \(\mu\) is set to mu_max. If the distance to the boundary is between dist_min and dist_max, the value of \(\mu\) is interpolated between mu_min and mu_max. The type of this interpolation is determined by the parameter

smooth_mu = True

If this parameter is set to smooth_mu = True, then a smooth, cubic polynomial is used to interplate between mu_min and mu_max, which yields a continuously differentiable \(\mu\). If this is set to smooth_mu = False, then a linear interpolation is used, which only yields a continuous \(\mu\). The default for this parameter is smooth_mu = False.

Finally, we can specify which boundaries we want to incorporate when computing the distance. To do so, we can specify a list of indices which contain the boundary markers in the parameter

boundaries_dist = [1,2,3]

This means, that only boundaries marked with 1, 2, and 3 are considered for computing the distance, and all others are ignored. The default behavior is that all (outer) boundaries are considered.

There is also another possibility to compute the shape gradient in cashocs, namely using the \(p\)-Laplacian, as proposed by Müller, Kühl, Siebenborn, Deckelnick, Hinze, and Rung. In order to do so, we have the following line

use_p_laplacian = False

If this is set to use_p_laplacian = True, the \(p\)-Laplacian is used to compute the shape gradient, as explained in Shape Optimization with the p-Laplacian. However, by default this is disabled. The value of \(p\) which is then used is defined in the next line

p_laplacian_power = 6

which defaults to p_laplacian_power = 2, whenever the parameter is not defined. The higher \(p\) is chosen, the better the numerical are expected to be, but the numerical solution of the problem becomes more involved.

Finally, there is the possibility to use a stabilized weak form for the \(p\)-Laplacian operator, where the stabilization parameter can be defined in the line

p_laplacian_stabilization = 0.0

The default value of this parameter is p_laplacian_stabilization = 0.0. Note, that the parameter should be chosen comparatively small, i.e., significantly smaller than 1.

Furthermore, we have the parameter fixed_dimensions, which enables us to restrict the shape gradient to specific dimensions. It is set via

fixed_dimensions = []

In case fixed_dimensions = [], there is no restriction on the shape gradient. However, if fixed_dimensions = [i], then the i-th component of the shape gradient is set to 0, so that we have no deformation in the i-th coordinate direction. For example, if fixed_dimensions = [0, 2], we only have a deformation in the \(y\)-component of the mesh. The default is fixed_dimensions = [].

Section Regularization#

In this section, the parameters for shape regularizations are specified. For a detailed discussion of their usage, we refer to Regularization for Shape Optimization Problems.

First, we have the parameters factor_volume and target_volume. These are set via the lines

factor_volume = 0.0
target_volume = 3.14

They are used to implement the (target) volume regularization term

\[\frac{\mu_\text{vol}}{2} \left( \int_{\Omega} 1 \text{ d}x - \text{vol}_\text{des} \right)^2\]

Here, \(\mu_\text{vol}\) is specified via factor_volume, and \(\text{vol}_\text{des}\) is the target volume, specified via target_volume. The default behavior is factor_volume = 0.0 and target_volume = 0.0, so that we do not have a volume regularization.

The next line, i.e.,

use_initial_volume = True

determines the boolean flag use_initial_volume. If this is set to use_initial_volume = True, then not the value given in target_volume is used, but instead the volume of the initial geometry is used for \(\text{vol}_\text{des}\).

For the next two _typing of regularization, namely the (target) surface and (target) barycenter regularization, the syntax for specifying the parameters is completely analogous. For the (target) surface regularization we have

factor_surface = 0.0
target_surface = 1.0

These parameter are used to implement the regularization term

\[\frac{\mu_\text{surf}}{2} \left( \int_{\Gamma} 1 \text{ d}s - \text{surf}_\text{des} \right)^2\]

Here, \(\mu_\text{surf}\) is determined via factor_surface, and \(\text{surf}_\text{des}\) is determined via target_surface. The default values are given by factor_surface = 0.0 and target_surface = 0.0.

As for the volume regularization, the parameter

use_initial_surface = True

determines whether the target surface area is specified via target_surface or if the surface area of the initial geometry should be used instead. The default behavior is given by use_initial_surface = False.

Next, we have the curvature regularization, which is controlled by the parameter

factor_curvature = 0.0

This is used to determine the size of \(\mu_\text{curv}\) in the regularization term

\[\frac{\mu_\text{curv}}{2} \int_{\Gamma} \kappa^2 \text{ d}s,\]

where \(\kappa\) denotes the mean curvature. This regularization term can be used to generate more smooth boundaries and to prevent kinks from occurring.

Finally, we have the (target) barycenter regularization. This is specified via the parameters

factor_barycenter = 0.0
target_barycenter = [0.0, 0.0, 0.0]

and implements the term

\[\frac{\mu_\text{bary}}{2} \left\lvert \frac{1}{\text{vol}(\Omega)} \int_\Omega x \text{ d}x - \text{bary}_\text{des} \right\rvert^2\]

The default behavior is given by factor_barycenter = 0.0 and target_barycenter = [0,0,0], so that we do not have a barycenter regularization.

The flag

use_initial_barycenter = True

again determines, whether \(\text{bary}_\text{des}\) is determined via target_barycenter or if the barycenter of the initial geometry should be used instead. The default behavior is given by use_initial_barycenter = False.


The object target_barycenter has to be a list. For 2D problems it is also sufficient, if the list only has two entries, for the \(x\) and \(y\) barycenters.

Finally, we have the parameter use_relative_scaling which is set in the line

use_relative_scaling = False

This boolean flag does the following. For some regularization term \(J_\text{reg}(\Omega)\) with corresponding factor \(\mu\) (as defined above), the default behavior is given by use_relative_scaling = False adds the term \(\mu J_\text{reg}(\Omega)\) to the cost functional, so that the factor specified in the configuration file is actually used as the factor for the regularization term. In case use_relative_scaling = True, the behavior is different, and the following term is added to the cost functional: \(\frac{\mu}{\left\lvert J_\text{reg}(\Omega_0) \right\rvert} J_\text{reg}(\Omega)\), where \(\Omega_0\) is the initial guess for the geometry. In particular, this means that the magnitude of the regularization term is equal to \(\mu\) on the initial geometry. This allows a detailed weighting of multiple regularization terms, which is particularly useful in case the cost functional is also scaled (see Scaling of the Cost Functional).

Section MeshQuality#

This section details the parameters that influence the quality of the computational mesh. First, we have the lines

volume_change = inf
angle_change = inf

These parameters are used to specify how much the volume and the angles, respectively, of the mesh elements are allowed to change in a single transformation. In particular, they implement the following criteria (see Etling, Herzog, Loayza, Wachsmuth, First and Second Order Shape Optimization Based on Restricted Mesh Deformations)

\[\begin{split}\frac{1}{\alpha} &\leq \det\left( \text{id} + D\mathcal{V} \right) \leq \alpha \\ \left\lvert\left\lvert D\mathcal{V} \right\rvert\right\rvert_{F} &\leq \beta.\end{split}\]

Here, \(\alpha\) corresponds to volume_change and \(\beta\) corresponds to angle_change, and \(\mathcal{V}\) is the deformation. The default behavior is given by volume_change = inf and angle_change = inf, so that no restrictions are posed. Note, that, e.g., Etling, Herzog, Loayza, Wachsmuth, First and Second Order Shape Optimization Based on Restricted Mesh Deformations use the values volume_change = 2.0 and angle_change = 0.3.

The next two parameters are given byx

tol_lower = 0.0
tol_upper = 1e-15

These parameters specify a kind of interval for the mesh quality. In particular, we have the following situation (note that the mesh quality is always an element in \([0,1]\)):

  • If the mesh quality is in \([\texttt{tol upper}, 1]\), the mesh is assumed to be “good”, so that finite element solutions of the corresponding PDEs are sensible and not influenced by the mesh quality or discretization artifacts.

  • If the mesh quality is in \([\texttt{tol lower}, \texttt{tol upper}]\), a kind of breaking point is reached. Here, it is assumed that the mesh is sufficiently good so that the solution of the state system is still possible. However, a mesh whose quality is in this interval should not be used anymore to compute the solution of the adjoint system or to compute the shape gradient, as the quality is too poor for this purpose. Usually, this means that the algorithm is terminated, unless remeshing is enabled. In the latter case, remeshing is performed.

  • If the mesh quality is in the interval \([0, \texttt{tol lower}]\), the mesh quality is assumed to be so poor, that even the solution of the state system is not possible anymore. In practice, this can only happen during the Armijo line search. Thanks to our previous considerations, we also know that the mesh, that is to be deformed, has at least a quality of tol_lupper, so that this quality might be reached again, if the step size is just decreased sufficiently often. This way, it is ensured that the state system is only solved when the mesh quality is larger than tol_lower, so that the corresponding cost functional value is reasonable.

The default behavior is given by tol_lower = 0.0 and tol_upper = 1e-15, so that there are basically no requirements on the mesh quality.

Finally, the upcoming two parameters specify how exactly the mesh quality is measured. The first one is

measure = condition_number

and determines one of the four mesh quality criteria, as defined in MeshQuality. Available options are

  • measure = skewness

  • measure = maximum_angle

  • measure = radius_ratios

  • measure = condition_number

(see MeshQuality for a detailed description). The default value is given by measure = skewness.

The parameter type determines, whether the minimum quality over all elements (type = min) or the average quality over all elements (type = avg) shall be used. This is set via

type = min

and defaults to type = min.

Finally, we have the parameter remesh_iter in which the user can specify after how many iterations a remeshing should be performed. It is given by

remesh_iter = 0

where remesh_iter = 0 means that no automatic remeshing is performed (this is the default), and remesh_iter = n means that remeshing is performed after each n iterations. Note that to use this parameter and avoid unexpected results, it might be beneficial to the the lower and upper mesh quality tolerances to a low value, so that the “quality based remeshing” does not interfere with the “iteration based remeshing”, but both can be used in combination.

Section Output#

In this section, the parameters for the output of the algorithm, either in the terminal or as files, are specified. First, we have the parameter verbose. This is used to toggle the output of the optimization algorithm. It defaults to verbose = True and is controlled via

verbose = True

The parameter save_results is a boolean flag, which determines whether a history of the optimization algorithm, including cost functional value, gradient norm, accepted step sizes, and mesh quality, shall be saved to a .json file. This defaults to save_results = True, and can be set with

save_results = False

Moreover, we define the parameter save_txt

save_txt = False

This saves the output of the optimization, which is usually shown in the terminal, to a .txt file, which is human-readable.

The next line in the config file is

save_state = False

Here, the parameter save_state is set. This is a boolean flag, which can be set to save_state = True to enable that cashocs generates .xdmf files for the state variables for each iteration the optimization algorithm performs. These are great for visualizing the steps done by the optimization algorithm, but also need some disc space, so that they are disabled by default. Note, that for visualizing these files, you need Paraview.

The next parameter, save_adjoint works analogously, and is given in the line

save_adjoint = False

If this is set to True, cashocs generates .xdmf files for the adjoint variables in each iteration of the optimization algorithm. Its main purpose is for debugging.

The next parameter is given by save_gradient, which is given in the line

save_gradient = False

This boolean flag ensures that a paraview with the computed shape gradient is saved in result_dir/xdmf. The main purpose of this is for debugging.

Moreover, we also have the parameter save_mesh that is set via

save_mesh = False

This is used to save the optimized geometry to a GMSH file. The default behavior is given by save_mesh = False. Note, that this is only possible if the input mesh was already generated by GMSH, and specified in the Mesh section of the config file. For any other meshes, the underlying mesh is also saved in the .xdmf files, so that you can at least always visualize the optimized geometry.

In the end, we also have, like for optimal control problems, a parameter that specifies where the output is placed, again named result_dir, which is given in the config file in the line

result_dir = ./results

As before, this is either a relative or absolute path to the directory where the results should be placed.

The parameter precision, which is set via

precision = 3

is an integer parameter which determines how many significant digits are printed in the output to the console and / or the result file.

Moreover, we have the parameter time_suffix, which adds a suffix to the result directory based on the current time. It is controlled by the line

time_suffix = False


Finally, an overview over all parameters and their default values can be found in the following.


Parameter = Default value



Only needed for remeshing


Only needed for remeshing


Only needed for remeshing

remesh = False

if remesh = True, remeshing is enabled; this feature is experimental, use with care

show_gmsh_output = False

if show_gmsh_output = True, shows the output of GMSH during remeshing in the console


Parameter = Default value


is_linear = False

using is_linear = True gives an error for nonlinear problems

newton_rtol = 1e-11

relative tolerance for Newton’s method

newton_atol = 1e-13

absolute tolerance for Newton’s method

newton_iter = 50

maximum iterations for Newton’s method

newton_damped = False

if newton_damped = True, damping is enabled

newton_inexact = False

if newton_inexact = True, an inexact Newton’s method is used

newton_verbose = False

newton_verbose = True enables verbose output of Newton’s method

picard_iteration = False

picard_iteration = True enables Picard iteration; only has an effect for multiple variables

picard_rtol = 1e-10

relative tolerance for Picard iteration

picard_atol = 1e-12

absolute tolerance for Picard iteration

picard_iter = 50

maximum iterations for Picard iteration

picard_verbose = False

picard_verbose = True enables verbose output of Picard iteration


Parameter = Default value



has to be specified by the user; see solve

rtol = 1e-3

relative tolerance for the optimization algorithm

atol = 0.0

absolute tolerance for the optimization algorithm

max_iter = 100

maximum iterations for the optimization algorithm

gradient_method = direct

specifies the solver for computing the gradient, can be either gradient_method = direct or gradient_method = iterative

gradient_tol = 1e-9

the relative tolerance in case an iterative solver is used to compute the gradient.

soft_exit = False

if soft_exit = True, the optimization algorithm does not raise an exception if it did not converge


Parameter = Default value


method = armijo

method = armijo is a simple backtracking line search, whereas method = polynomial uses polynomial models to compute trial stepsizes.

initial_stepsize = 1.0

initial stepsize for the first iteration in the Armijo rule

epsilon_armijo = 1e-4

beta_armijo = 2.0

safeguard_stepsize = True

De(-activates) a safeguard against poor scaling

polynomial_model = cubic

This specifies, whether a cubic or quadratic model is used for computing trial stepsizes

factor_high = 0.5

Safeguard for stepsize, upper bound

factor_low = 0.1

Safeguard for stepsize, lower bound

fail_if_not_converged = False

if this is True, then the line search fails if the state system can not be solved at the new iterate


Parameter = Default value


bfgs_memory_size = 5

memory size of the L-BFGS method

use_bfgs_scaling = True

if use_bfgs_scaling = True, uses a scaled identity mapping as initial guess for the inverse Hessian

bfgs_periodic_restart = 0

specifies, after how many iterations the method is restarted. If this is 0, no restarting is done.

damped = False

specifies whether damping for the BFGS method should be used to enforce the curvature condition and prevent restarting


Parameter = Default value


cg_method = FR

specifies which nonlinear CG method is used

cg_periodic_restart = False

if cg_periodic_restart = True, enables periodic restart of NCG method

cg_periodic_its = 10

specifies, after how many iterations the NCG method is restarted, if applicable

cg_relative_restart = False

if cg_relative_restart = True, enables restart of NCG method based on a relative criterion

cg_restart_tol = 0.25

the tolerance of the relative restart criterion, if applicable


Parameter = Default value


shape_bdry_def = []

list of indices for the deformable boundaries

shape_bdry_fix = []

list of indices for the fixed boundaries

shape_bdry_fix_x = []

list of indices for boundaries with fixed x values

shape_bdry_fix_y = []

list of indices for boundaries with fixed y values

shape_bdry_fix_z = []

list of indices for boundaries with fixed z values

fixed_dimensions = []

a list of coordinates which should be fixed during the shape optimization (x=0, y=1, etc.)

use_pull_back = True

if use_pull_back = False, shape derivative might be wrong; no pull-back for the material derivative is performed; only use with caution

lambda_lame = 0.0

value of the first Lamé parameter for the elasticity equations

damping_factor = 0.0

value of the damping parameter for the elasticity equations

mu_def = 1.0

value of the second Lamé parameter on the deformable boundaries

mu_fix = 1.0

value of the second Lamé parameter on the fixed boundaries

use_sqrt_mu = False

if use_sqrt_mu = True, uses the square root of the computed mu_lame; might be good for 3D problems

inhomogeneous = False

if inhomogeneous = True, uses inhomogeneous elasticity equations, weighted by the local mesh size

update_inhomogeneous = False

if update_inhomogeneous = True and inhomogeneous=True, then the weighting with the local mesh size is updated as the mesh is deformed.

inhomogeneous_exponent = 1.0

The exponent for the inhomogeneous mesh stiffness

use_distance_mu = False

if use_distance_mu = True, the value of the second Lamé parameter is computed via the distance to the boundary

dist_min = 1.0

Specifies the distance to the boundary, until which \(\mu\) is given by mu_min

dist_max = 1.0

Specifies the distance to the boundary, until which \(\mu\) is given by mu_max

mu_min = 1.0

The value of \(\mu\) for a boundary distance smaller than dist_min

mu_max = 1.0

The value of \(\mu\) for a boundary distance larger than dist_max

boundaries_dist = []

The indices of the boundaries, which shall be used to compute the distance, boundaries_dist = [] means that all boundaries are considered

smooth_mu = False

If false, a linear (continuous) interpolation between mu_min and mu_max is used, otherwise a cubic \(C^1\) interpolant is used

use_p_laplacian = False

If use_p_laplacian = True, then the \(p\)-Laplacian is used to compute the shape gradient

p_laplacian_power = 2

The parameter \(p\) of the \(p\)-Laplacian

p_laplacian_stabilization = 0.0

The stabilization parameter for the \(p\)-Laplacian problem. No stabilization is used when this is p_laplacian_stabilization = 0.0.


Parameter = Default value


factor_volume = 0.0

value of the regularization parameter for volume regularization; needs to be non-negative

target_volume = 0.0

prescribed volume for the volume regularization

use_initial_volume = False

if use_initial_volume = True uses the volume of the initial geometry as prescribed volume

factor_surface = 0.0

value of the regularization parameter for surface regularization; needs to be non-negative

target_surface = 0.0

prescribed surface for the surface regularization

use_initial_surface = False

if use_initial_surface = True uses the surface area of the initial geometry as prescribed surface

factor_curvature = 0.0

value of the regularization parameter for curvature regularization; needs to be non-negative

factor_barycenter = 0.0

value of the regularization parameter for barycenter regularization; needs to be non-negative

target_barycenter = [0.0, 0.0, 0.0]

prescribed barycenter for the barycenter regularization

use_initial_barycenter = False

if use_initial_barycenter = True uses the barycenter of the initial geometry as prescribed barycenter

use_relative_scaling = False

if use_relative_scaling = True, the regularization terms are scaling so that they have the magnitude specified in the respective factor for the initial iteration.


Parameter = Default value


volume_change = inf

determines by what factor the volume of a cell is allowed to change within a single deformation

angle_change = inf

determines how much the angles of a cell are allowed to change within a single deformation

tol_lower = 0.0

if the mesh quality is lower than this tolerance, the state system is not solved for the Armijo rule, instead step size is decreased

tol_upper = 1e-15

if the mesh quality is between tol_lower and tol_upper, the state system will still be solved for the Armijo rule. If the accepted step yields a quality lower than this, algorithm is terminated (or remeshing is initiated)

measure = skewness

determines which quality measure is used

type = min

determines if minimal or average quality is considered


Parameter = Default value


verbose = True

if verbose = True, the history of the optimization is printed to the console

save_results = True

if save_results = True, the history of the optimization is saved to a .json file

save_txt = True

if save_txt = True, the history of the optimization is saved to a human readable .txt file

save_state = False

if save_state = True, the history of the state variables over the optimization is saved in .xdmf files

save_adjoint = False

if save_adjoint = True, the history of the adjoint variables over the optimization is saved in .xdmf files

save_gradient = False

if save_gradient = True, the history of the shape gradient over the optimization is saved in .xdmf files

save_mesh = False

if save_mesh = True, saves the mesh for the optimized geometry; only available for GMSH input

result_dir = ./results

path to the directory, where the output should be placed

precision = 3

number of significant digits to be printed

time_suffix = False

Boolean flag, which adds a suffix to result_dir based on the current time