From 3406739288cf46694cfc8cfec0888d31e02887cc Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Sun, 23 May 2021 22:34:13 +0200 Subject: [PATCH] Some minor changes to the step-79 documentation --- examples/step-79/doc/intro.dox | 54 ++++++++++++++++++---------------- examples/step-79/step-79.cc | 8 ++--- 2 files changed, 32 insertions(+), 30 deletions(-) diff --git a/examples/step-79/doc/intro.dox b/examples/step-79/doc/intro.dox index df186bfbe0..2d737e5bf3 100644 --- a/examples/step-79/doc/intro.dox +++ b/examples/step-79/doc/intro.dox @@ -6,7 +6,7 @@ structure that is bearing some load. Ideally, we would like to minimize the maximum stress placed on a structure by selecting a region $E$ where material is placed. In other words, @f[ - \text{Minimize}\| \boldsymbol{\sigma} (\mathbf{u}) \|_\infty + \text{minimize}\| \boldsymbol{\sigma} (\mathbf{u}) \|_\infty @f] @f[ \text{subject to } |E|\leq V_{\max}, @@ -35,7 +35,7 @@ measure of the potential energy stored in an object due to its deformation, but also works as a measure of total deformation over the structure. @f[ - \text{Minimize } \int_E \frac{1}{2}\boldsymbol{\sigma} : \boldsymbol{\varepsilon} dV + \text{minimize } \int_E \frac{1}{2}\boldsymbol{\sigma} : \boldsymbol{\varepsilon} dV @f] @f[ \text{subject to } \|E\| \leq V_{\max} @@ -91,10 +91,10 @@ optimization problem. Use of SIMP then turns the optimization problem into the following: @f[ - \text{Minimize } \int_\Omega \frac{1}{2}\boldsymbol{\sigma}(\rho) : \boldsymbol{\varepsilon}(\rho) d_\Omega + \text{minimize } \int_\Omega \frac{1}{2}\boldsymbol{\sigma}(\rho) : \boldsymbol{\varepsilon}(\rho) d\Omega @f] @f[ - \text{subject to } \int_\Omega \rho(x) d_\Omega= V_{\max}, + \text{subject to } \int_\Omega \rho(x) d\Omega= V_{\max}, @f] @f[ 0<\rho_{\min}\leq \rho(x) \leq 1, @@ -103,7 +103,7 @@ following: \nabla \cdot \boldsymbol{\sigma}(\rho) + \mathbf{F} = 0 \quad \text{on } \Omega @f] -The final constraint, the balance of linear momentum (sometimes referred to as the elasticity equation), +The final constraint, the balance of linear momentum (which we will refer to as the elasticity equation), gives a method for finding $\boldsymbol{\sigma}$ and $\boldsymbol{\varepsilon}$ given the density $\rho$.

Elasticity Equation

@@ -123,7 +123,7 @@ $\mathbf{u}$ and $\boldsymbol{\sigma}$ as solution variables (as is done in mixe formulations). Furthermore, we will make the assumption that the material is linear isotropic, -in which case the stress-strain tensor can be expressed in terms of the Lam\'{e} +in which case the stress-strain tensor can be expressed in terms of the Lamé parameters $\lambda,\mu$ such that @f{align} \boldsymbol{\sigma} &= \rho^p (\lambda \text{tr}(\boldsymbol{\varepsilon}) \mathbf{I} + 2 \mu \boldsymbol{\varepsilon}) , \\ @@ -242,21 +242,21 @@ the objective function and the constraints by multiplying the constraints by Lagrange multipliers. Specifically, we will use the following symbols for the Lagrange multipliers for the various constraints:
    -
  1. $\mathbf{y}_1 $ - a Lagrange multiplier corresponding to the - elasticity constraint
    2. -
  2. $y_2$ - a Lagrange multiplier corresponding to the convolution - filter constraint
    3. -
  3. $z_1$ - a Lagrange multiplier corresponding to the lower slack variable
    4. -
  4. $z_2$ - a Lagrange multiplier corresponding to the upper slack variable.
    5. +
  5. $\mathbf{y}_1 $: a Lagrange multiplier corresponding to the + elasticity constraint,
    6. +
  6. $y_2$: a Lagrange multiplier corresponding to the convolution + filter constraint,
    7. +
  7. $z_1$: a Lagrange multiplier corresponding to the lower slack variable, and
    8. +
  8. $z_2$: a Lagrange multiplier corresponding to the upper slack variable.
With these variables, the Lagrangian function reads as follows: @f{align}{ \mathcal{L} =& \int_{\partial\Omega} \mathbf{u} \cdot \mathbf{t} d\partial\Omega - - \alpha \int_\Omega \left(\log(s_1) + \log(s_2)\right) d\Omega- \left(\int_\Omega + - \alpha \int_\Omega \left(\log(s_1) + \log(s_2)\right) d\Omega- \int_\Omega \rho^p \left(\frac{\mu}{2}\left(\boldsymbol{\varepsilon}(\mathbf{y}_1):\boldsymbol{\varepsilon}(\mathbf{u})) \right) + \lambda \left( \nabla \cdot \mathbf{u} \nabla \cdot \mathbf{y}_1 - \right) d\Omega \right)- \int_{\partial \Omega} \mathbf{y}_1 \cdot \mathbf{t} d\partial\Omega \right) \\ + \right)\right) d\Omega - \int_{\partial \Omega} \mathbf{y}_1 \cdot \mathbf{t} d\partial\Omega \\ & -\int_\Omega y_2 (\rho - H(\varrho)) d\Omega - \int_\Omega z_1 (\varrho-s_1) d\Omega - \int_\Omega z_2 (1 - s_2 -\varrho) d\Omega @f} @@ -345,13 +345,14 @@ relationships are also inequalities. We will address the nonlinearity using a Newton method to compute search directions, and come back to how to deal with the inequalities below when talking about step length procedures. -Newton's method applied to the equations above results in the following system. -Here, variational derivatives with respect to the $\{\bullet\}$ variable are taken - in the $c_{\{\bullet\}}$ direction. This gives +Newton's method applied to the equations above results in the system of equations +listed below. +Therein, variational derivatives with respect to the $\{\bullet\}$ variable are +taken in the $c_{\{\bullet\}}$ direction.
    -
  1. Stationarity - these equations ensure we are at a critical point of the -objective function when constrained +
  2. Stationarity: These equations ensure we are at a critical point of the +objective function when constrained. Equation 1 @f{align}{ @@ -387,7 +388,7 @@ Equation 3 @f]
    3. -
  3. Primal Feasibility - these equations ensure the equality constraints +
  4. Primal Feasibility: These equations ensure the equality constraints are met. Equation 4 @@ -421,8 +422,8 @@ Equation 7 @f]
    5. -
  5. Complementary Slackness - these equations essentially ensure the barrier -is met - in the final solution, we need $s^T z = 0$ +
  6. Complementary Slackness: These equations essentially ensure the barrier +is met - in the final solution, we need $s^T z = 0$. Equation 8 @f[ @@ -437,9 +438,9 @@ Equation 9 @f]
    7. -
  7. Dual Feasibility - Multiplier on slacks and slack variables must be kept +
  8. Dual Feasibility: The Lagrange multiplier on slacks and slack variables must be kept greater than 0. (This is the only part not implemented in the -SANDTopOpt::assemble_system() function) +`SANDTopOpt::assemble_system()` function.) @f[ s,z \geq 0 @f] @@ -564,7 +565,7 @@ penalty parameter should be in the next iteration, rather than using reduction factors that are independent of the current solution. That said, such methods are also more complicated and we will not do this here. -

    Merit Function

    +

    Merit %Function

    The algorithm outlined above makes use of a "merit function". Merit functions are used to determine whether a step from $x_k$ to a proposed point $x_{k+1}$ is @@ -596,7 +597,8 @@ there exists some $p_0$ so that for any $p > p_0$, the merit function has its minima at the same location as the original problem. This penalty parameter is updated (by recommendation of Nocedal and Wright @cite Benson2002) as follows: @f[ - p > \frac{\frac{1}{2} \mathbf{x}^T \cdot \mathbf{H} \cdot \mathbf{x} - \mathbf{x}^T \cdot \nabla f}{\|c_i\|_{l_\infty}, i \in \mathcal{E}} , + p > \frac{\frac{1}{2} \mathbf{x}^T \cdot \mathbf{H} \cdot \mathbf{x} - \mathbf{x}^T \cdot \nabla f}{\|c_i\|_{l_\infty}} + \quad , i \in \mathcal{E}, @f] where $\mathbf{H}$ is the Hessian of the objective function, $\mathbf{x}$ is a vector of our decision (primal) variables, $f$ is the objective function, and $c_i$ is the error on a diff --git a/examples/step-79/step-79.cc b/examples/step-79/step-79.cc index d89e0ab170..4a2f50236e 100644 --- a/examples/step-79/step-79.cc +++ b/examples/step-79/step-79.cc @@ -464,7 +464,7 @@ namespace SAND // The use of the symbolic names defined in namespace // `SolutionComponents` helps understand what each of the // following terms corresponds to, but it also makes the - // expressions lengthy and unwieldy: An term such as + // expressions lengthy and unwieldy: A term such as // `coupling[SolutionComponents::density_upper_slack_multiplier][SolutionComponents::density]` // just doesn't read very well, and would either have to be // split over several lines or run off the right edge of @@ -792,9 +792,9 @@ namespace SAND // us how far off our filtered density is from the filter // applied to the unfiltered density. That is because while at // the solution of the nonlinear problem, we have - // $\rho=H\sigma$, but at intermediate iterations, we in - // general have $\rho^k\neq H\sigma^k$ and the "residual" - // $\rho^k-H\sigma^k$ will then appear as the right hand side + // $\rho=H\varrho$, but at intermediate iterations, we in + // general have $\rho^k\neq H\varrho^k$ and the "residual" + // $\rho^k-H\varrho^k$ will then appear as the right hand side // of one of the Newton update equations that we compute // below. BlockVector filtered_unfiltered_density_solution = -- 2.39.5