Add the parameter "inhomogeneous_exponent" to the Section ShapeGradient

CHANGELOG.rst

@@ -32,6 +32,10 @@ of the maintenance releases, please take a look at
 
     * :ini:`fail_if_not_converged` determines, whether the line search is cancelled once the state system cannot be solved or if a new iterate is tried instead.
 
+  * Section ShapeGradient
+
+    * :ini:`inhomogeneous_exponent` is a float, which specifies an exponent for the inhomogeneous mesh stiffness
+
   * Section MeshQuality
 
     * :ini:`remesh_iter` is used to perform a remeshing after a certain amount of iterations.

cashocs/_forms/shape_form_handler.py

@@ -567,10 +567,15 @@ def _compute_shape_gradient_forms(self) -> ufl.Form:
                 self.volumes.vector().vec().scale(1 / vol_max)
                 self.volumes.vector().apply("")
 
+                self.inhomogeneous_exponent = fenics.Constant(
+                    self.config.getfloat("ShapeGradient", "inhomogeneous_exponent")
+                )
             else:
                 self.volumes.vector().vec().set(1.0)
                 self.volumes.vector().apply("")
 
+                self.inhomogeneous_exponent = fenics.Constant(1.0)
+
             def eps(u: fenics.Function) -> ufl_expr.Expr:
                 """Computes the symmetric gradient of a vector field ``u``.
 
@@ -589,16 +594,16 @@ def eps(u: fenics.Function) -> ufl_expr.Expr:
             riesz_scalar_product = (
                 fenics.Constant(2)
                 * self.mu_lame
-                / self.volumes
+                / pow(self.volumes, self.inhomogeneous_exponent)
                 * fenics.inner(eps(trial), eps(test))
                 * self.dx
                 + fenics.Constant(lambda_lame)
-                / self.volumes
+                / pow(self.volumes, self.inhomogeneous_exponent)
                 * fenics.div(trial)
                 * fenics.div(test)
                 * self.dx
                 + fenics.Constant(damping_factor)
-                / self.volumes
+                / pow(self.volumes, self.inhomogeneous_exponent)
                 * fenics.inner(trial, test)
                 * self.dx
             )

cashocs/io/config.py

@@ -325,6 +325,10 @@ def __init__(self, config_file: Optional[str] = None) -> None:
                 "update_inhomogeneous": {
                     "type": "bool",
                 },
+                "inhomogeneous_exponent": {
+                    "type": "float",
+                    "attributes": ["non_negative"],
+                },
                 "use_distance_mu": {
                     "type": "bool",
                 },
@@ -578,7 +582,6 @@ def __init__(self, config_file: Optional[str] = None) -> None:
 [ShapeGradient]
 lambda_lame = 0.0
 damping_factor = 0.0
-inhomogeneous = False
 mu_def = 1.0
 mu_fix = 1.0
 use_sqrt_mu = False
@@ -594,7 +597,9 @@ def __init__(self, config_file: Optional[str] = None) -> None:
 dist_max = 1.0
 boundaries_dist = []
 smooth_mu = False
+inhomogeneous = False
 update_inhomogeneous = False
+inhomogeneous_exponent = 1.0
 fixed_dimensions = []
 shape_bdry_def = []
 shape_bdry_fix = []

docs/source/user/demos/shape_optimization/doc_config.rst

@@ -673,6 +673,16 @@ Moreover, the parameter
 
 can be used to update the local mesh size after each mesh deformation, in case this is :ini:`update_inhomogeneous = True`, so that elements which become smaller also obtain a higher stiffness and vice versa. The default is :ini:`update_inhomogeneous = False`.
 
+For the inhomogeneous mesh stiffness, we also have the parameter :ini:`inhomogeneous_exponent`, which is specified via
+
+.. code-block:: ini
+
+    inhomogeneous_exponent = 1.0
+
+This parameter can be used to specify an exponent for the inhomogeneous mesh stiffness, so that the parameters
+:math:`\mu, \lambda` and :math:`\delta` are scaled by :math:`\left( \frac{1}{\text{vol} \right)^p`, where
+:math:`p` is specified in :ini:`inhomogeneous_exponent`. The default for this parameter is :ini:`inhomogeneous_exponent = 1.0`.
+
 There is also a different possibility to define the stiffness parameter :math:`\mu`
 using cashocs, namely to define :math:`\mu` in terms of how close a point of the
 computational domain is to a boundary. In the following we will explain this
@@ -1271,6 +1281,8 @@ in the following.
       - if :ini:`inhomogeneous = True`, uses inhomogeneous elasticity equations, weighted by the local mesh size
     * - :ini:`update_inhomogeneous = False`
       - if :ini:`update_inhomogeneous = True` and :ini:`inhomogeneous=True`, then the weighting with the local mesh size is updated as the mesh is deformed.
+    * - :ini:`inhomogeneous_exponent = 1.0`
+      - The exponent for the inhomogeneous mesh stiffness
     * - :ini:`use_distance_mu = False`
       - if :ini:`use_distance_mu = True`, the value of the second Lamé parameter is computed via the distance to the boundary
     * - :ini:`dist_min = 1.0`