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

Note

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.

Note

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]

Note

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.

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 = [].

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.

Warning

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.

Note

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.

Note

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.

Moreover, we have the parameter degree_estimation which is specified via

degree_estimation = True

This parameter enables cashocs’ default estimation of the quadrature degree for the shape derivative. If this is set to False, an error related to FEniCS may occur - so this should be always enabled.

Next, we have the parameter global_deformation which is set via the line

global_deformation = False

If this is set to True, cashocs computes the deformation from the initial to the optimized mesh (even when remeshing has been performed). This can, however, lead to some unexpected errors with PETSc, so this should be used with care.

We have the parameter test_for_intersections, which is specified via

test_for_intersections = True

If this parameter is set to True, cashocs will check the deformed meshes for (self) intersections, which would generate non-physical geometries and reject them - so that all generated designs are physically meaningful. This should not be set to False.

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.

Hint

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 MeshQualityConstraints

The parameter min_angle is used to define the threshold angle, i.e. the minimum (dihedral) angle which is feasible for the mesh. The default is given by

min_angle = 0.0

which ensures that the constraints are not active by default. Note that the initial mesh has to be feasible for the method to work, so if the minimum angle in the mesh is smaller than the min_angle specified in the configuration, cashocs will raise an exception. Note that the angle is specified in degree and not radians.

To circumvent this problem for meshes with small angles (which could be used, e.g., to resolve boundary layers, the next parameter feasible_angle_reduction_factor is used. This parameter specifies, how much smaller the (dihedral) angles of the mesh are allowed to become relative to the value in the initial mesh. That means a value of feasible_angle_reduction_factor = 0.25 ensures that no (dihedral) angle in a mesh element will become smaller than one quarter of the smallest angle of the element in the initial mesh. The default is given by

feasible_angle_reduction_factor = 0.0

which ensures that the constraints are not active by default.

Note

If both the feasible_angle_reduction_factor and min_angle are given, cashocs uses the element-wise minimum of the two. In particular, this means that a strategy of using feasible_angle_reduction_factor = 0.9999 and some value for min_angle can be used to constrain the (dihedral) angle to a specific value, wherever this is possible (and leave the angles that are below this threshold as they are).

The parameter tol is used to define a tolerance for which constraints are treated as active or not. As we treat the constraints numerically, they can only be satisfied up to a certain tolerance, which the user can specify here. The default value of

tol = 1e-2

should work well in most situations. In some situations, the optimization could be faster with a tolerance of tol = 1e-1 (but should never be larger) or more accurate when using, e.g., tol = 1e-3 (lower values should most of the time not be necessary).

The parameter mode can only be set to

mode = approximate

at the moment, which is also the default value. In the future, other options might be possible.

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 mesh as a GMSH file in each iteration of the optimization. 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

Summary

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

[Mesh]

Parameter = Default value	Remarks
mesh_file	Only needed for remeshing
gmsh_file	Only needed for remeshing
geo_file	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

[StateSystem]

Parameter = Default value	Remarks
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

[OptimizationRoutine]

Parameter = Default value	Remarks
algorithm	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

[LineSearch]

Parameter = Default value	Remarks
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

[AlgoLBFGS]

Parameter = Default value	Remarks
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

[AlgoCG]

Parameter = Default value	Remarks
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

[ShapeGradient]

Parameter = Default value	Remarks
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.
global_deformation = False	Computes the global deformation from initial to optimized mesh. This can lead to unexpected errors, use with care.
degree_estimation = True	Estimate the required degree for quadrature of the shape derivative. This should be True, otherwise unexpected errors can happen.
test_for_intersections = True	If enabled, the mesh is tested for intersections which would create non-physical meshes. This should always be enabled, otherwise the obtained results might be incorrect.

[Regularization]

Parameter = Default value	Remarks
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.

[MeshQuality]

Parameter = Default value	Remarks
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
remesh_iter	When set to a value > 0, remeshing is performed every remesh_iter iterations.

[MeshQualityConstraints]

Parameter = Default value	Remarks
min_angle = 0.0	The minimum feasible triangle / dihedral angle of the mesh cells in degrees. This is constant for all cells. If this is positive, the constraints are used. If this is 0, no constraints are used.
tol = 1e-2	The tolerance for the mesh quality constraints. If abs(g(x)) < tol, then the constraint is considered active
mode = approximate	The mode for calculating the (shape) derivatives of the constraint functions. At the moment, only “approximate” is supported.
feasible_angle_reduction_factor = 0.0	A factor in the interval [0,1) which sets the feasible reduction of the triangle / dihedral angles. This means, that each cell is only allowed to have angles larger than this times the initial minimum angle.

[Output]

Parameter = Default value	Remarks
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 in each iteration of the optimization; 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