diff --git a/shapo/antora.html b/shapo/antora.html deleted file mode 100644 index 016ed80..0000000 --- a/shapo/antora.html +++ /dev/null @@ -1,372 +0,0 @@ - - - - - - Antora environment :: Feel++ ShapO // Docs - - - - - - - - - - - - - - - - - - - - - - -
- -
-
- -
- -
-
- - -
- -
-

Antora environment

-
-

Antora is our static website generator. -We use it to generate the documentation of the project and upload it to docs.feelpp.org[Feel++ docs] website.

-
-
-
To generate the documentation
-
-
cd docs
-antora site.yml
-
-
-
-
To visualizethe documentation
-
-
# from topevel directory
-node-srv docs/public
-# from docs/
-node-srv docs
-
-
-
-
-
-
- - - - - - - - - - - - diff --git a/shapo/cantilever/index.html b/shapo/cantilever/index.html new file mode 100644 index 0000000..ffa72c5 --- /dev/null +++ b/shapo/cantilever/index.html @@ -0,0 +1,1001 @@ + + + + + + Linear-elasticity : Cantilever :: Feel++ ShapO // Docs + + + + + + + + + + + + + + + + + + + + + + + + +
+ +
+
+ +
+ +
+
+ + +
+ +
+

Linear-elasticity : Cantilever

+
+
+
+

The examples are based on Lucas Palazzolo’s course. For more details, see [internship_Palazzolo].

+
+
+
+
+

1. Introduction

+
+
+

We consider a homogeneous isotropic elastic solid which occupies a bounded domain \(\Omega\) in \(\mathbb{R}^2\). We suppose that the boundary is composed of three parts

+
+
+
+\[\begin{equation*} +\partial \Omega = \Gamma \cup \Gamma_N \cup \Gamma_D +\end{equation*}\] +
+
+
+

where \(\Gamma\) is a variable part traction-free (homogeneous Neumann), \(\Gamma_D\) is a fixed part on which the solid is fixed (homogeneous Dirichlet) and \(\Gamma_N\) is a fixed part on which \(f\) forces are applied (inhomogeneous Neumann). The cantilever problem is illustrated in the following image.

+
+
+
+cantilever +
+
Figure 1. Illustration of the 2D cantilever with a hole. The \(\Gamma_D\) and \(\Gamma_N\) boundaries (Dirichlet condition and Neumann condition) are fixed. The \(\Gamma\) boundary (of which the hole is a part) is movable in order to optimize the shape.
+
+
+
+
Cantilever PDE
+
+

The displacement \(u : \mathbb{R}^N \to \mathbb{R}^N\) is the solution of the system

+
+
+
+\[\begin{equation} +\begin{cases} +-\nabla \cdot \sigma = 0 &\qquad \text{in } \Omega,\\ +\sigma = 2\mu e(u) + \lambda tr(e(u))I &\qquad \text{in } \Omega,\\ +u=0 &\qquad \text{on } \Gamma_D,\\ +\sigma n = f &\qquad \text{on } \Gamma_N,\\ +\sigma n = 0 &\qquad \text{on } \Gamma, + \end{cases} +\end{equation}\] +
+
+
+

with \(\sigma(u)\in M_N(\mathbb{R})\) the stress tensor, \(\lambda, \mu >0\) the Lamé coefficients of the material and \(e(u)=\frac{1}{2}\left(\nabla u + \nabla u^T\right)\) the deformation tensor.

+
+
+
+
+
+
Cantilever cost function
+
+

We wish to solve the following minimisation problem

+
+
+
+\[\begin{equation}\label{cantilever:minpb} +\inf_{\Omega \in \Omega_{ad}} \left\{ J(\Omega) = \int_{\Gamma_N} f \cdot u \right\} +\end{equation}\] +
+
+
+

where the set of admissible forms is defined as

+
+
+
+\[\begin{equation}\label{cantilever:thetaad} +\Omega_{ad}=\left\{\Omega \in C(\Omega_0) \mid \Gamma_D \cup \Gamma_N \subset \partial \Omega, ~ g(\Omega)=0\right\} +\end{equation}\] +
+
+
+

with \(g\) defined by [eq:Omegad].

+
+
+
+
+ + + + + +
+ + +
+

It should be noted that the domain of integration of the cost function does not depend on \(\Omega\) because \(\Gamma_N\) is fixed. Thus, the dependence with respect to the domain is ensured only through the solution of the previous model.

+
+
+
+
+
+
+

2. Theoretical formulation of the problem

+
+
+
+
The weak formulation of Cantilever PDE
+
+

The weak formulatio of the cantilever problem is given by

+
+
+
+\[\begin{equation} +\int_{\Omega}2\mu e(u):e(\phi) + \int_{\Omega}\lambda (\nabla \cdot u)(\nabla \cdot \phi) = \int_{\Gamma_N}f \cdot \phi. +\end{equation}\] +
+
+
+
+
+Proof +
+
+

The divergence of \(\sigma\) is given by

+
+
+
+\[ \begin{equation*} +\nabla \cdot \sigma = \left(\sum_{j=1}^N \frac{\partial \sigma_{ij}}{\partial x_j}\right)_{1\leq i\leq N}. +\end{equation*}\] +
+
+
+

Let \(u \in H^1_{\Gamma_N}(\Omega)\) solution of Cantilever PDE. By using the fact that \(tr(e(u))=\nabla \cdot u\), the equation can be written for \(1\leq i\leq N\) as

+
+
+
+\[\begin{equation*} + -\sum_{j=1}^N \frac{\partial}{\partial x_j}\left[\mu\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)+\lambda (\nabla \cdot u)\delta_{ij}\right]=0. +\end{equation*}\] +
+
+
+

By multiplying by a test function \(\phi \in H^1_{\Gamma_D}(\Omega)\) and integrating over the domain \(\Omega\), we have

+
+
+
+\[\begin{equation*} + \int_{\Omega}\sum_{j=1}^N\mu\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_j}\right)\frac{\partial \phi_i}{\partial x_j}+\int_{\Omega}\lambda (\nabla \cdot u) \frac{\partial \phi_i}{\partial x_i}=\int_{\Gamma_N}f_i \phi_i. +\end{equation*}\] +
+
+
+

In addition, we have the following result

+
+
+
+\[\begin{equation}\label{ipp:utils1} +\sum_{i,j=1}^{N}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\frac{\partial \phi_i}{\partial x_j}=\frac{1}{2}\sum_{i,j=1}^{N}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\left(\frac{\partial \phi_i}{\partial x_j}+\frac{\partial \phi_j}{\partial x_j}\right)=2e(u):e(\phi) +\end{equation}\] +
+
+
+

Finally, summing over the index \(i\) in the previous formulation and using this result, we obtain the result.

+
+
+
+
+
+
Definition : Lagrangian of Cantilever PDE and Cantilever cost function
+
+

We denote the Lagrangian of this problem by

+
+
+
+\[\begin{equation*} +\begin{array}{rcl} +\mathcal{L}:\Omega_{ad}\times \left(H^1_{\Gamma_D}(\mathbb{R}^N)\right)^2&\to& \mathbb{R}\\ +(\Omega, v,q) &\mapsto &\int_{\Omega}2\mu e(v):e(q) + \int_{\Omega} \lambda (\nabla \cdot v)(\nabla \cdot q) - \int_{\Gamma_N} f \cdot q + \int_{\Gamma_N} f \cdot v. +\end{array} +\end{equation*}\] +
+
+
+
+
+ + + + + +
+ + +
+

It should be noted that the different variables of \(\mathcal{L}\) are independent because we consider Lagrange multipliers on \(H^1_{\Gamma_D}(\mathbb{R}^N)\) and not \(H^1_{\Gamma_D}(\Omega)\) and especially that \(\Gamma_N\) is independent of \(\Omega\) because it is fixed.

+
+
+
+
+
+
Dual PDE of the cantilever problem
+
+

The dual PDE of the cantilever problem is

+
+
+
+\[\begin{equation*} +\begin{cases} +-\nabla \cdot \sigma = 0 &\qquad \text{on } \Omega\\ +\sigma = 2\mu e(p) + \lambda tr(e(p))I &\qquad \text{on } \Omega\\ +p=0 &\qquad \text{on } \Gamma_D\\ +\sigma n = -f &\qquad \text{on } \Gamma_N\\ +\sigma n = 0 &\qquad \text{on } \Gamma + \end{cases} +\end{equation*}\] +
+
+
+
+
+Proof +
+
+

The partial derivative with respect to \(v\) : let \(\phi \in H^1_{\Gamma_D}(\mathbb{R}^N)\)

+
+
+
+\[\begin{equation*} +\left\langle \frac{\partial \mathcal{L}}{\partial v}(\Omega, v, q), \phi \right\rangle = \int_{\Omega}2\mu e(\phi):e(q) + \int_{\Omega} \lambda (\nabla \cdot \phi) (\nabla \cdot q) + \int_{\Gamma_N} f \cdot \phi, +\end{equation*}\] +
+
+
+

which, when it vanishes, corresponds to the weak formulation of the dual problem.

+
+
+
+
+ + + + + +
+ + +
+

Note that this corresponds exactly to the primal problem with a minus sign. This is due in particular to the choice of \(J\) and the boundary conditions of the primal problem. Indeed, we have \(J\) which depends on \(f\) on \(\Gamma_N\) and the boundary conditions are all null except on \(\Gamma_N\) which is equal to \(f\). Thus, we can avoid solving the dual problem in this case and simply use the solution of the primal problem.

+
+
+
+
+
+
Shape gradient for the cantilever problem
+
+

The shape gradient for the cantiler problem is defined by

+
+
+
+\[\begin{equation}\label{cantilever:gradshape} + G(\Omega)=-2\mu\|e(u)\|^2-\lambda(tr(e(u)))^2. +\end{equation}\] +
+
+
+
+
+Proof +
+
+

Let’s calculate the partial derivative of \(\mathcal{L}\) with respect to the \(\Omega\) domain in the \(\theta\) direction

+
+
+
+\[\begin{equation*} + \frac{\partial \mathcal{L}}{\partial \Omega}(\Omega, v, q)(\theta)= \int_{\partial \Omega}\theta \cdot n \left[ 2\mu e(v):e(q) + \lambda (\nabla \cdot v)(\nabla \cdot q)\right] +\end{equation*}\] +
+
+
+

by using [Prop:diffJOmega]. When we evaluate this derivative with the state \(u(\Omega)\) and the adjoint state \(p(\Omega)=u(\Omega)\), we find exactly the value of the derivative of the cost function

+
+
+
+\[\begin{equation}\label{cantilever:shaped} +\frac{\partial \mathcal{L}}{\partial \Omega}\left(\Omega, u(\Omega), p(\Omega)\right)(\theta)= -\int_{\Gamma}\theta \cdot n \left[ 2\mu \|e(u)\|^2 + \lambda tr(e(u))^2\right] = DJ(\Omega)(\theta), +\end{equation}\] +
+
+
+

as explained here [lagmethod:gradientJ]. Thus by definition we obtain the shape gradient.

+
+
+
+
+
+
+

3. Experimental Evaluation

+
+
+

In this section, we will present various results on the optimization problem concerning the shape of cantilevers. The initial domain considered is a trapezoid with two feet in the 2D case (four feet in the 3D case). The choice of a trapezoid as the initial shape is advantageous due to its simplicity and its ability to approximate the solution of the problem. Throughout this section, unless specified otherwise, all results have been obtained using \(\mathbb{P}_1\) continuous finite elements. All simulations are carried out using the classic gradient descent method. The NSGF method is used in the following example, a rigid body in a Stokes fluid.

+
+
+

3.1. 2D simulations for various types of cantilever

+
+

We present the application of shape optimization to a 2D cantilever with two pillars. The specific parameter values used in the analysis are detailed in the following table.

+
+ + ++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Table 1. Geometric parameters and Lamé coefficient of the initial domain for the 2D case of the cantilever. The Lamé coefficients are \(\lambda\) and \(\mu\). \(H\) represents the height of the trapezoid, \(L_1\) the large base and \(L_2\) the small base. The force applied to the curve \(\Gamma_N\) is defined by \(f\).

Symbol

Value (dimensionless)

\(\lambda\)

\(50/9\)

\(\mu\)

\(350/27\)

\(H\)

\(9\)

\(L_1\)

\(8\)

\(L_2\)

\(2\)

\(f\)

\((0,-1)\)

+
+

No hole :

+
+
+

We begin by considering the simplest case, which is the cantilever without any holes, as depicted in Figures below. To ensure clear visibility of the mesh in print, a discretization parameter of \(h=0.4\) is set. For the initialization of the Lagrange multiplier, we choose \(l=0.5\). Additionally, the values of \(a=b=0.5\) and \(c=10\) are set, along with a descent step of \(t=0.2\). These specific values have been determined through empirical testing to achieve optimal performance.

+
+ + +++++++ + + + + + + + + + + + + + + + + +
Table 2. Comparison of the cost function between the initial domain and the final domain, the volume error between the two domains, the \(H^1\)-norm of the displacement field at the end and the number of iterations for the 2D cantilever without hole.

\(J(\Omega_0)\)

\(J(\Omega_{n_{final}})\)

\(||\Omega_0|-|\Omega_{n_{final}}||\)

\(||\theta_{n_{final}}||_{H^1}\)

\(n_{final}\)

\(4.6834\)

\(3.1346\)

\(1.3630e-2\)

\(1.696e-2\)

\(125\)

+
+
+results nohole 2D l0.5 t0.2 h0.4 a0.5b0.5 c10 n200 +
+
Figure 2. Evolution of the cost function, the volume of the domain \(\Omega_k\) and the norm of \(\theta_k\) of the 2D cantilever without hole with the following parameters: \(h=0.4\), \(t=0.2\), \(l=0.5\), \(a=b=0.5\) and \(c=10\). The iterations range from stem[0] to \(125\).
+
+ + +++ + + + + + +
Table 3. Visualisation of the results obtained during different iterations of the algorithm for the shape optimisation of the 2D cantilever without hole with the following parameters: \(h=0.4\), \(t=0.2\), \(l=0.5\), \(a=b=0.5\) and \(c=10\). The initial domain (\(k=0\)) is shown on the left. The intermediate domain (\(k=60\)) is displayed in the middle. The final domain (\(k=125\)) is displayed on the right.

nohole 0 +nohole 60 +nohole 125

+
+

One hole :

+
+
+

The second test consists of studying the cantilever with a single hole. The results obtained are presented in Figures below. We use a discretization parameter of \(h=0.4\). For the Lagrange multiplier, we initialize with \(l=0.5\), and set \(a=b=0.5\) and \(c=10\) with a descent step of \(t=0.2\).

+
+ + +++++++ + + + + + + + + + + + + + + + + +
Table 4. Comparison of the cost function between the initial domain and the final domain, the volume error between the two domains, the \(H^1\)-norm of the displacement field at the end and the number of iterations for the one hole 2D cantilever.

\(J(\Omega_0)\)

\(J(\Omega_{n_{final}})\)

\(||\Omega_0|-|\Omega_{n_{final}}||\)

\(||\theta_{n_{final}}||_{H^1}\)

\(n_{final}\)

\(4.7035\)

\(3.2162\)

\(2.9980e-3\)

\(1.8906e-2\)

\(125\)

+
+
+results hole 2D l0.5 t0.2 h0.4 a0.5b0.5 c10 n200 +
+
Figure 3. Evolution of the cost function, the volume of the domain \(\Omega_k\) and the norm of \(\theta_k\) of the 2D cantilever with one hole with the following parameters: \(h=0.\), \(t=0.2\), \(l=0.5\), \(a=b=0.5\) and \(c=10\). The iterations range from \(0\) to \(125\).
+
+ + +++ + + + + + +
Table 5. Visualisation of the results obtained during different iterations of the algorithm for the shape optimisation of the 2D cantilever with one hole with the following parameters: \(h=0.4\), \(t=0.2\), \(l=0.5\), \(a=b=0.5\) and \(c=10\). The initial domain (\(k=0\)) is shown on the left. The intermediate domain (\(k=60\)) is displayed in the middle. The final domain (\(k=125\)) is displayed on the right.

hole 0 +hole 60 +hole 125

+
+

Four holes :

+
+
+

The last type of cantilever studied is one with four holes in its domain. The results are displayed in the figures below. We use a discretization parameter of \(h=0.4\). For the Lagrange multiplier, we initialize with \(l=0.5\), and set \(a=b=0.5\) and \(c=10\) with a descent step of \(t=0.2\).

+
+ + +++++++ + + + + + + + + + + + + + + + + +
Table 6. Comparison of the cost function between the initial domain and the final domain, the volume error between the two domains, the \(H^1\)-norm of the displacement field at the end and the number of iterations for the four holes 2D cantilever.

\(J(\Omega_0)\)

\(J(\Omega_{n_{final}})\)

\(||\Omega_0|-|\Omega_{n_{final}}||\)

\(||\theta_{n_{final}}||_{H^1}\)

\(n_{final}\)

\(5.1589\)

\(3.3770\)

\(2.5317e-3\)

\(1.5615e-2\)

\(125\)

+
+
+results 4holes 2D l0.5 t0.2 h0.4 a0.5b0.5 c10 n200 +
+
Figure 4. Evolution of the cost function, the volume of the domain \(\Omega_k\) and the norm of \(\theta_k\) of the 2D cantilever with 4 holes with the following parameters: \(h=0.4\), \(t=0.2\), \(l=0.5\), \(a=b=0.5\) and \(c=10\). The iterations range from \(0\) to \(125\).
+
+ + +++ + + + + + +
Table 7. Visualisation of the results obtained during different iterations of the algorithm for the shape optimisation of the 2D cantilever with 4 holes with the following parameters: \(h=0.4\), \(t=0.2\), \(l=0.5\), \(a=b=0.5\) and \(c=10\). The initial domain (\(k=0\)) is shown on the left. The intermediate domain (\(k=60\)) is displayed in the middle. The final domain (\(k=125\)) is displayed on the right.

4hole 0 +4hole 60 +4hole 125

+
+

The results presented above demonstrate that satisfactory outcomes can be achieved using the proposed method. Notably, the obtained shapes closely resemble those reported in various papers that focus on optimizing the geometrical shape of cantilevers, such as [allaire_conception_2006]. Furthermore, we successfully decrease the cost function while approaching the initial volume of the domain, which aligns precisely with the desired behavior. It is worth noting that potential geometric instabilities, in the form of spikes, may emerge on \(\Gamma_D\) after a certain number of iterations. These spikes are also observed in the 3D case, as illustrated in the subsequent figures.

+
+
+
+

3.2. 3D simulations for various types of cantilever :

+
+

In this section, we address the case of the 3D cantilever, which introduces additional complexities compared to the 2D case and results in longer computation times. Incorporating holes into the structure also increases the risk of encountering mesh superposition issues. The overall geometry of the 3D cantilever remains similar to the 2D case, with the addition of four pillars, and boundary conditions are now applied to surfaces instead of curves. The specific parameter values used in the analysis are provided in detail in following table.

+
+ + ++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Table 8. Geometric parameters and Lamé coefficient of the initial domain for the 2D case of the cantilever. The Lamé coefficients are \(\lambda\) and \(\mu\). \(H\) represents the height of the truncated pyramid, \(C_1\) the side of the largest square base and \(C_2\) the side of the smallest square base. The force applied to the surface \(\Gamma_N\) is defined by \(f\).

Symbol

Value (dimensionless)

\(\lambda\)

\(50/9\)

\(\mu\)

\(350/27\)

\(H\)

\(9\)

\(C_1\)

\(8\)

\(C_2\)

\(2\)

\(f\)

\((0,-1,0)\)

+
+

No hole :

+
+
+

Let’s first consider the case without a hole. We use a discretization parameter of \(h=0.8\). For the coefficients affecting the optimization problem, we set \(l=0.5\), \(a=b=0.5\), \(c=3\), and choose a descent step \(t\) of \(0.1\). The results of this test are presented in Figures below.

+
+ + +++++++ + + + + + + + + + + + + + + + + +
Table 9. Comparison of the cost function between the initial domain and the final domain, the volume error between the two domains, the \(H^1\)-norm of the displacement field at the end and the number of iterations for the no hole 3D cantilever.

\(J(\Omega_0)\)

\(J(\Omega_{n_{final}})\)

\(||\Omega_0|-|\Omega_{n_{final}}||\)

\(||\theta_{n_{final}}||_{H^1}\)

\(n_{final}\)

\(4.8060\)

\(3.0484\)

\(0.1439\)

\(2.1647e-2\)

\(400\)

+
+
+results nohole 3D l0.5 t0.1 h0.8 a0.5b0.5 c2 n500 +
+
Figure 5. Evolution of the cost function, the volume of the domain \(\Omega_k\) and the norm of \(\theta_k\) of the 3D cantilever without hole with the following parameters: \(h=0.8\), \(t=0.1\), \(l=0.5\), \(a=b=0.5\) and \(c=2\). The iterations range from \(0\) to \(400\).
+
+ + +++ + + + + + +
Table 10. Visualisation of the results obtained during different iterations of the algorithm for the shape optimisation of the 3D cantilever without hole with the following parameters: \(h=0.8\), \(t=0.1\), \(l=0.5\), \(a=b=0.5\) and \(c=2\). The initial domain (\(k=0\)) is shown on the left. The intermediate domain (\(k=150\)) is displayed in the middle. The final domain (\(k=400\)) is displayed on the right.

3Dnohole 0 +3Dnohole 120 +3Dnohole 400

+
+

One hole :

+
+
+

In the next test, we add a hole inside the material. The results are presented in Figures below. We use a discretization parameter of \(h=0.8\). The other coefficients used are \(l=0.5\), \(a=b=0.5\), \(c=2\), and \(t=0.08\).

+
+ + +++++++ + + + + + + + + + + + + + + + + +
Table 11. Comparison of the cost function between the initial domain and the final domain, the volume error between the two domains, the \(H^1\)-norm of the displacement field at the end and the number of iterations for the one hole 3D cantilever.

\(J(\Omega_0)\)

\(J(\Omega_{n_{final}})\)

\(||\Omega_0|-|\Omega_{n_{final}}||\)

\(||\theta_{n_{final}}||_{H^1}\)

\(n_{final}\)

\(6.4813\)

\(3.7726\)

\(0.2582\)

\(4.7253e-2\)

\(250\)

+
+
+results hole 3D l0.5 t0.08 h0.8 a0.5b0.5 c2 n500 +
+
Figure 6. Evolution of the cost function, the volume of the domain \(\Omega_k\) and the norm of \(\theta_k\) of the 3D cantilever with holes with the following parameters: \(h=0.8\), \(t=0.08\), \(l=0.5\), \(a=b=0.5\) and \(c=2\). The iterations range from \(0\) to \(250\).
+
+ + +++ + + + + + +
Table 12. Visualisation of the results obtained during different iterations of the algorithm for the shape optimisation of the 3D cantilever with holes with the following parameters: \(h=0.8\), \(t=0.08\), \(l=0.5\), \(a=b=0.5\) and \(c=2\). The initial domain (\(k=0\)) is shown on the left. The intermediate domain (\(k=150\)) is displayed in the middle. The final domain (\(k=250\)) is displayed on the right.

3Dhole 0 +3Dhole 120 +3Dhole 250

+
+

It is important to note that the chosen optimal shape in the simulations effectively reduces the cost function while maintaining the volume. However, due to limitations, we had to prematurely halt the simulations. If allowed to continue, the cost function would have been further minimized, resulting in a more optimal shape. Nonetheless, the obtained shape presents several issues. Specifically, certain areas such as the feet and top part exhibit more prominent outgrowths, which can be attributed to volume conservation. The shape optimization process tends to hollow out the cantilever below the desired volume (between the four legs). As a result, to compensate for this volume loss, protrusions seem to emerge after a certain number of iterations. Another significant issue encountered is mesh collision and overlap. Additionally, the discontinuities observed in the curves of the different simulations correspond to the remeshing that occurs during the iterations.

+
+
+
+
+
+

References on geometric shape optimisation

+
+
+
    +
  • +

    [internship_Palazzolo] Lucas Palazzolo. Shape optimization for rigid objects in Stokes flow. Stage de M2, Univerité de Strasbourg/Inria Sophia-Antipolis, 2023. repo

    +
    • +
  • +

    [allaire_conception_2006] Grégoire Allaire. Conception optimale de structures. volume 58 of Mathématiques & applications. Springer Berlin Heidelberg, 2006

    +
    • +
  • +

    [cea_conception_1986] Jean Céa. Conception optimale ou identification de formes, calcul rapide de la dérivée directionnelle de la fonction coût. M2AN - Modélisation mathématique et analyse numérique, 20(3):371-402, 1986

    +
    • +
  • +

    [moreau_shapes_2022] Clément Moreau, Kenta Ishimoto and Yannick Privat. Shapes optimising grand resistance tensor entries for a rigid body in a Stokes flow. July 2022

    +
    • +
  • +

    [feppon_shape_2019] Florian Feppon. Shape and topology optimization of multiphysics systems. These de doctorat, Université Paris-Saclay (ComUE), Dec. 2019

    +
    • +
+
+
+
+
+
+
+
+ + + + + + + + + + + + diff --git a/shapo/cmake.html b/shapo/cmake.html deleted file mode 100644 index c5d61db..0000000 --- a/shapo/cmake.html +++ /dev/null @@ -1,366 +0,0 @@ - - - - - - Cmake environment :: Feel++ ShapO // Docs - - - - - - - - - - - - - - - - - - - - - - -
- -
-
- -
- -
-
- - -
- -
-

Cmake environment

-
-

The CMakeLists.txt is setup with Feel++. -A sample Feel++ application is available in src.

-
-
-
    -
  • -

    cpack is configured

    -
    • -
  • -

    ctest is configured

    -
    • -
-
-
-
-
-
- - - - - - - - - - - - diff --git a/shapo/geometricshapeopti.html b/shapo/geometricshapeopti.html deleted file mode 100644 index 545d85e..0000000 --- a/shapo/geometricshapeopti.html +++ /dev/null @@ -1,3217 +0,0 @@ - - - - - - Geometric Shape Optimisation :: Feel++ ShapO // Docs - - - - - - - - - - - - - - - - - - - - - - -
- -
-
- -
- -
-
- - -
- -
-

Geometric Shape Optimisation

- -
-
-
-

This introduction to geometric shape optimisation is based on Lucas Palazzolo’s course. For more details, see [internship_Palazzolo].

-
-
-

Shape optimization is a field of study that focuses on finding the best possible shape for a given object or system in order to optimize certain performance criteria or objectives which are typically a function of differential equations defined on the shape itself. The optimisation problem can be written as

-
-
-
-\[\begin{equation*} - \inf_{\Omega \in \Omega_{ad}}J(\Omega, u(\Omega)) -\end{equation*}\] -
-
-
-

with \(u\) solution of a PDE defined on the domain \(\Omega\).

-
-
-

Geometric shape optimization focuses on directly modifying the shape of an object to achieve the desired objectives via a non-parametric deformation field \(\theta\). It involves manipulating the geometry of the object to improve its performance. We consider a reference domain \(\Omega_0\), which we assume to be a regular bounded open subset of \(\mathbb{R}^N\). Let \(\partial \Omega_0\) the boundary of this domain. We denote with \(\theta\) the deformation applied to the initial domain to get the new domain \(\Omega\), i.e.

-
-
-
-\[\begin{equation*} - \Omega = (\textrm{Id} + \theta)(\Omega_0). -\end{equation*}\] -
-
-
-

This deformation is illustrated in the following figure where a deformation field is applied to a reference domain.

-
-
-
-shape opti principle -
-
Figure 1. Shape optimisation principle : the \(\Omega_0\) domain is deformed according to a deformation field \(\theta\) such that the new \(\Omega\) domain is given by \(\Omega = (\textrm{Id} + \theta)(\Omega_0)\).
-
-
-

In the pursuit of finding shape derivative of the cost function, a well-known method called Cea’s method, developed by Cea in [cea_conception_1986], proves to be both simple and effective. However, it should be noted that Cea’s method is considered a formal approach as it relies on certain assumptions regarding the regularity of the problem’s data. We look at the numerical solution of these two problems using two numerical methods : the gradient descent [allaire_conception_2006] and the null space gradient flow [feppon_shape_2019].

-
-
-

The gradient descent method is a classic method for solving geometric shape optimisation problems. In this method, an initial shape is iteratively modified by moving in the direction of the negative gradient of the objective function. The process continues until convergence to an optimal shape is achieved. The null space gradient flow method involves finding the shape that minimizes the objective function while satisfying a set of constraints. It utilizes the notion of null space and range space directions to guide the shape optimization process.

-
-
-
-
-

1. Introduction

-
-
-

1.1. Definitions

-
-

As presented in [allaire_conception_2006], in order to describe the set of admissible forms, we need to introduce the following set of diffeomorphisms.

-
-
-
-
Definition : Deformation of the identity
-
-

We define a space of diffeomorphisms on \(\mathbb{R}^N\) (which can be seen as a deformation of the identity) by

-
-
-
-\[\begin{equation*} -\mathcal{T} = \left\{ T \mid (T-\textrm{Id})\in W^{1,\infty}(\Omega, \mathbb{R}^N) \text{ and } (T^{-1}-\textrm{Id})\in W^{1,\infty}(\Omega, \mathbb{R}^N) \right\}. -\end{equation*}\] -
-
-
-
-
-

Now we can clarify the admissible shapes obtained by deformation of \(\Omega_0\), that are defined by the following set.

-
-
-
-
Definition : The set of admissible shapes
-
-

The set of admissible shapes obtained by deformation of \(\Omega_0\) is given by

-
-
-
-\[\begin{equation*} - C(\Omega_0)=\left\{\Omega \subset \mathbb{R}^N \mid \exists T \in \mathcal{T}, \quad \Omega=T(\Omega_0) \right\}, -\end{equation*}\] -
-
-
-

where \(\mathcal{T}\) is defined by Definition : Deformation of the identity.

-
-
-
-
-

In view of the above definitions, it is natural to consider the vector field \(\theta\) such as

-
-
-
-\[\begin{equation*} - T = \textrm{Id} + \theta \qquad \text{with} \qquad \theta \in W^{1,\infty}(\Omega, \mathbb{R}^N) ~ \text{and} ~ T \in \mathcal{T}. -\end{equation*}\] -
-
-
-

In order to be able to define a differentiability notion in \(\Omega_0\) with respect to \(\theta\), we need the following lemma which guarantees that if \(\theta\) is small enough then \(T= \textrm{Id} + \theta\) is indeed a diffeomorphism which belongs to \(\mathcal{T}\).

-
-
-
-
Lemma 1
-
-

For all \(\theta \in W^{1,\infty}(\Omega, \mathbb{R}^N)\) that verify \(\|\theta\|_{W^{1,\infty}(\Omega, \mathbb{R}^N)}<1\), the map \(T=\textrm{Id} +\theta\) is a bijection of \(\mathbb{R}^N\) that belongs to \(\mathcal{T}\) defined by Definition : Deformation of the identity.

-
-
-
-
-Proof -
-
-

See [[allaire_conception_2006], Lemma 6.13]

-
-
-
-
-

All that remains is to define the notion of differentiability in relation to the domain that we call shape derivation or differentiability with respect to the domain. The following definition applies to differentiability in the Fréchet sense as well as in the Gâteaux sense.

-
-
-
-
Definition : Differentiability with respect to the domain
-
-

Let \(J(\Omega)\) a map such that \(J : C(\Omega_0) \to \mathbb{R}\). \(J\) is said to be differentiable with respect to the domain on \(\Omega_0\) if

-
-
-
-\[\begin{equation*} - \theta \mapsto J\left((\textrm{Id} +\theta)(\Omega_0)\right) -\end{equation*}\] -
-
-
-

is differentiable in 0 in the Banach space \(W^{1,\infty}(\Omega, \mathbb{R}^N).\)

-
-
-
-
-

Let us introduce some notation to define the matrix scalar product.

-
-
-
-
Definition : Matrix scalar product
-
-

Let \(A\) and \(B\) be two real matrices of dimension \(N\). We define the scalar product of two matrices as

-
-
-
-\[\begin{equation*} - A : B = tr(A^TB). -\end{equation*}\] -
-
-
-

and we will note thereafter \(\|\cdot \|\) the associated norm.

-
-
-
-
-
-
Definition
-
-

We denote

-
-
-
-\[\begin{equation*} - H^{k}_{\Gamma}(\Omega) = \left\{v \in H^{k}(\Omega) \mid v|_{\Gamma} = 0 \right\}, -\end{equation*}\] -
-
-
-

the set of functions belonging to \(H^{k}(\Omega)\) and which vanish on \(\Gamma\).

-
-
-
-
-
-

1.2. Derivation of integrals

-
-

In this subsection, Definition : Differentiability with respect to the domain of differentiability with respect to the domain will be applied to volume or surface integrals. The paragraph will be brief, for more details we point the reader to [allaire_conception_2006]. We denote \(n\) as the unit normal pointing outwards to the domain.

-
-
-
-
Proposition 1
-
-

Let \(\Omega_0\) be a regular bounded open subset of \(\mathbb{R}^N\). Let \(f \in W^{1,1}(\mathbb{R}^N)\) and \(J\) be the map of \(C(\Omega_0)\) in \(\mathbb{R}\) defined by

-
-
-
-\[\begin{equation*} - J(\Omega)=\int_{\Omega}f. -\end{equation*}\] -
-
-
-

Then \(J\) is differentiable in \(\Omega\) and for all \(\theta \in W^{1,\infty}(\Omega,\mathbb{R}^N)\) we have

-
-
-
-\[\begin{equation*} - DJ(\Omega)(\theta)=\int_{\Omega}\nabla \cdot \left(\theta f\right) = \int_{\partial \Omega}\theta\cdot n f. -\end{equation*}\] -
-
-
-
-
-Proof -
-
-

See [[allaire_conception_2006], Proposition 6.22]

-
-
-
-
-
-
Propostion 2
-
-

Let \(\Omega_0\) be a regular bounded open subset of \(\mathbb{R}^N\). Let \(f \in W^{2,1}(\mathbb{R}^N)\) and \(J\) be the map of \(C(\Omega_0)\) in \(\mathbb{R}\) defined by

-
-
-
-\[\begin{equation*} - J(\Omega)=\int_{\partial \Omega}f. -\end{equation*}\] -
-
-
-

Then \(J\) is differentiable in \(\Omega\) and for all \(\theta \in C^1(\mathbb{R}^N,\mathbb{R}^N)\) we have

-
-
-
-\[\begin{equation*} -DJ(\Omega)(\theta)=\int_{\partial \Omega}\left(\nabla f \cdot \theta + f(\nabla \cdot \theta)-f(\nabla\theta n \cdot n)\right) = \int_{\partial \Omega}\theta \cdot n \left(\frac{\partial f}{\partial n}+(\nabla \cdot n)f\right). -\end{equation*}\] -
-
-
-
-
-Proof -
-
-

See [[allaire_conception_2006], Proposition 6.24 & Lemma 6.25]

-
-
-
-
- - - - - -
- - -
-

Note that we need \(\theta \in C^1(\mathbb{R}^N,\mathbb{R}^N)\). Indeed, we need to define the trace of \(\nabla \theta\) on \(\partial \Omega\), and therefore the hypothesis \(\theta \in W^{1,\infty}(\Omega,\mathbb{R}^N)\) is not sufficient.

-
-
-
-
-
-

1.3. Derivation of domain-dependent functions and equations: Cea’s method

-
-

In the following, we consider a cost function \(J\) (depending on a domain \(\Omega\)) which we wish to minimize (or maximize) by finding the shape of the optimal domain. Additionally, we assume that \(J\) relies on a variable \(u\), which is a solution to a specific differential equation defined on the domain \(\Omega\). Consequently, it becomes necessary to differentiate the cost function and the differential equation with respect to the domain.

-
-
-

A simple and efficient method for calculating values is the Cea’s method developed in [cea_conception_1986]. This method also allows us to "guess" the definitions of the adjoint state \(p\). However, it is important to note that this method is formal, as it relies on certain assumptions about the regularities of the problem’s data, particularly concerning the solution \(u\). For more rigorous calculations, the use of Eulerian and Lagrangian derivatives, as described in [allaire_conception_2006], becomes necessary. This method incorporates the concepts of the primal problem or state problem, the dual problem or adjoint problem, and the Lagrangian. The problem is presented as follows.

-
-
-

Cost Function : Let \(J\) be a (domain-dependent) cost function that we seek to minimize (or maximize) such that

-
-
-
-\[\begin{equation*} - \begin{array}{rcl} -J:C(\Omega_0)&\to& \mathbb{R}\\ -\Omega &\mapsto &J(\Omega, u(\Omega)), -\end{array} -\end{equation*}\] -
-
-
-

where \(J\) depends of the solution \(u\) of a differential equation. For ease of reading, we omit the \(u\) in the notation of the cost function. We seek to solve the following problem

-
-
-
-\[\begin{equation*} - \inf_{\Omega \in \Omega_{ad}} J(\Omega), -\end{equation*}\] -
-
-
-

where \(\Omega_{ad}\) corresponds to the admissible domains.

-
-
-

Primal Equation : Let \(u\) be the solution of a problem, called a state problem or primal problem. Let \(E\) be the map governing the variational formulation associated with the primal problem such that

-
-
-
-\[\begin{equation*} -\begin{array}{rcl} -E:V\times V&\to& \mathbb{R}\\ -(v,q) &\mapsto &E(v,q), -\end{array} -\end{equation*}\] -
-
-
-

where \(V\) is, in our case, a Sobolev space. Then, by definition of the weak formulation, for all \(q \in V\) we have

-
-
-
-\[\begin{equation*} - E(u,q)=0, -\end{equation*}\] -
-
-
-

with \(u\) the solution of the primal problem.

-
-
-

The Lagrangian : The Cea’s method is based on duality. For this purpose, we consider the primal equation as a constraint and introduce the following Lagrangian.

-
-
-
-
Definition : Lagrangian without Dirichlet condition
-
-
-\[\begin{equation*} -\begin{array}{rcl} -\mathcal{L}:C(\Omega_0)\times V\times V&\to& \mathbb{R}\\ -(\Omega, v,q) &\mapsto &J(\Omega) + E(v,q), -\end{array} -\end{equation*}\] -
-
-
-

which is the sum of the objective function and the variational formulation of the primal equation.

-
-
-
-
- - - - - -
- - -
-

We must not forget that \(J\) also depends on \(v\) ! It is absolutely important for the following to consider that the three variables are independent of each other. Thus, the space \(V\) must be independent of the domain \(\Omega\).

-
-
-
-
-

When the state problem has a Dirichlet condition, we must add another Lagrange multiplier, as explained in [allaire_conception_2006].

-
-
-
-
Deifnition : Lagrangian with Dirichlet condition
-
-
-\[\begin{equation*} -\begin{array}{rcl} -\mathcal{L}:C(\Omega_0)\times V\times V\times V&\to& \mathbb{R}\\ -(\Omega, v,q, \psi) &\mapsto &J(\Omega) + E(v,q) + F(v,\psi), -\end{array} -\end{equation*}\] -
-
-
-

where \(F\) represents the constraint associated with the Dirichlet condition.

-
-
-
-
-

To simplify the following, we will assume that \(E\) and \(F\) are bilinear. We then quickly obtain the different equations and the gradient of \(J\) (assuming all the necessary regularities).

-
-
-
-
Proposition : Primal problem thanks to the Lagrangian
-
-

For all \((q,\phi) \in V^2\)

-
-
-
-\[\begin{equation} - \left\langle\frac{\partial \mathcal{L}}{\partial q}(\Omega, u, q),\phi \right\rangle = 0, -\end{equation}\] -
-
-
-

with \(u\) solution of the primal problem, i.e. for all \(\phi \in V\)

-
-
-
-\[\begin{equation} - E(u,\phi) = 0, -\end{equation}\] -
-
-
-

by bilinearity of \(E\). This results in \(u\) being the solution to the variational formulation of the primal problem.

-
-
-
-
-
-
Proposition : Dual problem thanks to the Lagrangian
-
-

For all \((v,\phi)\in V^2\), we have

-
-
-
-\[\begin{equation} - \left\langle \frac{\partial \mathcal{L}}{\partial v}(\Omega, v, p),\phi \right\rangle = 0, -\end{equation}\] -
-
-
-

with \(p\) solution of the dual problem, i.e. for all \(\phi \in V\)

-
-
-
-\[\begin{equation} - E(\phi, p)+ \frac{\partial J}{\partial v}(\Omega)(\phi) = 0, -\end{equation}\] -
-
-
-

by bilinearity of \(E\) and by the fact that \(J\) also depends on \(v\). This results in \(p\) being solution to the weak formulation of the dual problem.

-
-
-
-
-
-
Proposition : Differential of the cost fucntion thanks to the Lagrangian
-
-

For all \(\theta\) in V, we have

-
-
-
-\[\begin{equation} - DJ(\Omega)(\theta)=\frac{\partial \mathcal{L}}{\partial \Omega}(\Omega, u(\Omega), p(\Omega))(\theta). -\end{equation}\] -
-
-
-

with \(u\) and \(p\) solution of the primal and dual problem.

-
-
-
-
-Proof -
-
-

Since \(u(\Omega)\) is a solution of the primal problem and thus vanishes the weak formulation, we have

-
-
-
-\[\begin{equation*} - \mathcal{L}(\Omega, u(\Omega), q) = J(\Omega) \qquad \forall q\in V. -\end{equation*}\] -
-
-
-

The key point is that this is valid for all \(q\) in \(V\). Thus, we can take \(q\) in particular in the following as the solution of the dual problem (similarly for Lagrange multipliers in the case of Dirichlet conditions). To emphasize the dependence of the solution \(u\) on the domain, we write \(u(\Omega)\). By independence of the variable \(q\) with respect to \(\Omega\) and deriving this relation using chain rule (always assuming that we have the necessary regularities to do so), we get

-
-
-
-\[\begin{equation*} - DJ(\Omega)(\theta) = \frac{\partial \mathcal{L}}{\partial \Omega}\left(\Omega, u(\Omega), q\right)(\theta) + \left\langle \frac{\partial \mathcal{L}}{\partial v}\left(\Omega, u(\Omega), q\right), u'(\Omega)(\theta) \right\rangle. -\end{equation*}\] -
-
-
-

By taking \(q=p(\Omega)\) solution of the dual problem, the last term vanishes. Finally, we come to

-
-
-
-\[\begin{equation} - DJ(\Omega)(\theta)=\frac{\partial \mathcal{L}}{\partial \Omega}(\Omega, u(\Omega), p(\Omega))(\theta). -\end{equation}\] -
-
-
-
-
-
-
-
-

2. Solution methods

-
-
-

To minimize the cost function \(J(\Omega)\), we differentiate according to the variable \(\theta\) which parametrizes the form \(\Omega=(\textrm{Id}+\theta)(\Omega_0)\). In the various problems we study, certain boundaries will be assumed to be fixed, i.e. non-deformable, and noted \(\Gamma_{fixed}\). We denote \(\Gamma\) the deformable part of \(\partial \Omega\), thus

-
-
-
-\[\begin{equation*} - \partial \Omega = \Gamma \cup \Gamma_{fixed}. -\end{equation*}\] -
-
-
-

In [allaire_conception_2006] a method is quickly defined for all cost functions whose derivative can be written as follows

-
-
-
-\[\begin{equation}\label{eq:DJ} - DJ(\Omega)(\theta)=\int_{\Gamma} \theta \cdot n G(\Omega), -\end{equation}\] -
-
-
-

where \(G(\Omega)\) is a function (which depends of the state and of the adjoint) which is called the shape gradient. In the following, we will assume that we are studying an optimisation problem with an equality constraint in order to preserve the volume of the domain.

-
-
-
-
Definition : Set of admissible domains
-
-

The set of admissible domains, \(\Omega_{ad}\), is given by

-
-
-
-\[\begin{equation}\label{eq:Omegad} -\Omega_{ad}=\left\{\Omega \in C(\Omega_0) \mid \Gamma_{fixed} \subset \partial \Omega, ~ g(\Omega)= 0\right\} -\end{equation}\] -
-
-
-

with -\(g\) the equality constraint for the volume conservation, written as

-
-
-
-\[\begin{equation*} - \begin{array}{rcl} - g:V&\to& \mathbb{R}\\ - \theta &\mapsto &|(\textrm{Id} + \theta)(\Omega)|-|\Omega_0|. - \end{array} -\end{equation*}\] -
-
-
-
-
- - - - - -
- - -
-

By abuse of notation, the most of the time we write \(g(\Omega)\) instead of \(g(\theta)\).

-
-
-
-
-

The non-deformation constraints the boundaries \(\Gamma_{fixed}\) are taken into account by simply imposing

-
-
-
-\[\begin{equation*} -\theta=0 \qquad \text{on} \qquad \Gamma_{fixed}, -\end{equation*}\] -
-
-
-

i.e. we impose a null displacement field at these boundaries.

-
-
-
-
Proposition : Differential of the volume constraint
-
-

The differential of \(g\), definined in Definition : Set of admissible domains, is given by

-
-
-
-\[\begin{equation*} - \begin{array}{rcl} - Dg(\Omega):V&\to& \mathbb{R}\\ - \theta &\mapsto &\int_{\Gamma}\theta \cdot n. - \end{array} -\end{equation*}\] -
-
-
-
-
-Proof -
-
-

By using Proposition 1 on \(g\).

-
-
-
-
-

Two methods will be presented to solve this type of problems: gradient descent and Null Space Gradient Flow.

-
-
-

2.1. Solution by a gradient descent method

-
-

The gradient descent (GD) method is an optimization technique used to minimize a function iteratively. By taking a new form, \(\Omega_t\), such as

-
-
-
-\[\begin{equation*} -\Omega_t=(\textrm{Id}+\theta_t)(\Omega_0) \qquad \text{with} \qquad \theta_t = -t G(\Omega_0)n , -\end{equation*}\] -
-
-
-

where \(t>0\), we have

-
-
-
-\[\begin{equation*} - DJ(\Omega_0)(\theta_t)=-t\int_{\Gamma_0}G(\Omega_0)^2 <0. -\end{equation*}\] -
-
-
-

Thus, for a sufficiently small step size \(t\), we have

-
-
-
-\[\begin{equation*} - J(\Omega_{t}) < J(\Omega_0) \qquad \text{if} \qquad G(\Omega_0)\ne 0. -\end{equation*}\] -
-
-
-

It can be noted that to have \(DJ(\Omega_0)(\theta)<0\), we need to choose \(\theta \cdot n >0\). Let \(l \in \mathbb{R}\) a Lagrange multiplier for the volume constraint. The domain optimality condition, by using the Lagrange multiplier theorem, can be seen as

-
-
-
-\[\begin{equation}\label{eq:Lm} -DJ(\Omega)(\theta)+l Dg(\Omega)(\theta) = \int_{\Gamma} \theta \cdot n \left[l +G(\Omega)\right] = 0. -\end{equation}\] -
-
-
-

Numerical implementation : For the numerical implementation, a domain sequence, \(\Omega_k\), is computed which verifies the following constraints

-
-
-
-\[\begin{equation*} - \partial \Omega_k = \Gamma_k \cup \Gamma_{fixed}. -\end{equation*}\] -
-
-
-

Let \(t>0\) a fixed descent step, as explained in [allaire_conception_2006] the following iterations are performed until convergence, for \(k\geq 0\) :

-
-
-
-\[\begin{equation*} -\Omega_{k+1} = (\textrm{Id} + \theta_k)(\Omega_k) -\end{equation*}\] -
-
-
-

with

-
-
-
-\[\begin{equation*} -\theta_k = \begin{cases} --t\left[l_k+G(\Omega_k)\right]n_k &\qquad \text{on } \Gamma_k \\ -0 &\qquad \text{on } \Gamma_{fixed} - \end{cases} -\end{equation*}\] -
-
-
-

where \(n_k\) is the normal vector of \(\partial \Omega_k\) and \(l_k \in \mathbb{R}\) a Lagrange multiplier. To extend the trace of \(\theta_k\) on \(\partial \Omega_k\) inside \(\Omega_k\), we can solve the following system

-
-
-
-
Expansion PDE for GD
-
-
-\[\begin{equation*} -\begin{cases} --\Delta\theta_k = 0 &\qquad \text{on } \Omega_k \\ -\theta_k = 0 &\qquad \text{on } \Gamma_{fixed} \\ -\theta_k = -t\left[l_k+G(\Omega_k)\right]n_k &\qquad \text{on } \Gamma_k. - \end{cases} -\end{equation*}\] -
-
-
-
-
-

Once we know \(\theta_k\) about \(\Omega_k\), we can deform the whole mesh to obtain a new one of the domain \(\Omega_{k+1}\). It will be advisable to remesh the domain from time to time when the mesh becomes of poor quality measure in \(\texttt{Feel++}\) (which leads to numerical errors).

-
-
-

For a higher regularity, we can solve the following system

-
-
-
-
Regularised expansion PDE for GD
-
-
-\[\begin{equation*} -\begin{cases} --\Delta\theta_k = 0 &\qquad \text{on } \Omega_k \\ -\theta_k = 0 &\qquad \text{on } \Gamma_{fixed} \\ -\frac{\partial \theta_k}{\partial n} = -t\left[l_k+G(\Omega_k)\right]n_k &\qquad \text{on } \Gamma_k. - \end{cases} -\end{equation*}\] -
-
-
-
-
-
-
Proposition : Descent direction of Regularised expansion PDE for GD
-
-

The Regularised expansion PDE for GD still provides a descent direction.

-
-
-
-
-Proof -
-
-

We have

-
-
-
-\[\begin{align*} - DJ(\Omega_k)(\theta_k) + l_k Dg(\Omega_k)(\theta_k) &= \int_{\Gamma_k}\theta_k \cdot n_k \left[l_k + G(\Omega_k)\right]\\ - &=t^{-1}\int_{\Omega_k}\Delta\theta_k \cdot \theta_k - t^{-1}\int_{\Gamma_k}\theta_k \cdot \nabla \theta_k n_k, -\end{align*}\] -
-
-
-

because of the boundary conditions and the fact that \(\Delta \theta_k =0\). By integrating, we finally obtain that

-
-
-
-\[\begin{equation*} -DJ(\Omega_k)(\theta_k) + l_k Dg(\Omega_k)(\theta_k) =-t^{-1}\int_{\Omega_k}\|\nabla \theta_k\|^2 \leq 0. -\end{equation*}\] -
-
-
-
-
-

Concerning the choice of the implementation of the Lagrange multiplier, the same model as G.Allaire in his code will be taken into account (a kind of augmented Lagrangian). We thus have

-
-
-
-\[\begin{equation*} - l_k = al_{k-1}+b\frac{\int_{\Gamma_k}G(\Omega_k)}{|\Gamma_k|}+c\frac{|\Omega_k|-|\Omega_0|}{|\Omega_0|}, -\end{equation*}\] -
-
-
-

with \(a\), \(b\) and \(c\) parameters to be chosen according to the problem studied.

-
-
-
-

2.2. Null Space Gradient Flow

-
-

Null Space Gradient Flow (NSGF), presented in [feppon_shape_2019], is a new approach to solving constrained optimization problems. Unlike traditional methods, this technique does not require necessarily the adjustment of non-physical parameters, making it highly practical and reliable. Its applicability to shape optimization applications further enhances its usefulness.

-
-
-

The key concept behind NSGF is to modify the gradient flow algorithm to accommodate the presence of constraints. This is achieved by solving an Ordinary Differential Equation (ODE) that governs the optimization trajectories, denoted as \(x(t)\). The ODE takes the form:

-
-
-
-\[\begin{equation*} -\dot x = -\alpha_J \xi_J(x(t)) - \alpha_C \xi_C(x(t)) -\end{equation*}\] -
-
-
-

Here, \(\alpha_J\) and \(\alpha_C\) positive parameters, control the influence of the objective function and constraint violation, respectively. \(\xi_J(x)\) and \(\xi_C(x)\) represent the null space and range space directions, respectively. The null space direction, \(\xi_J(x)\), is obtained by projecting the gradient \(\nabla J(x)\) onto the cone of feasible directions, ensuring descent while respecting the constraints. The range space direction, \(\xi_C(x)\), guides the optimization path smoothly towards the feasible region.

-
-
-

We consider the case where the optimization takes place on a Hilbert space \(V\) with inner product \(\langle ., . \rangle_{V}\). The transpose of the differential of a function in infinite dimension is defined as follows.

-
-
-
-
Definition : Transpose
-
-

Let \(g : V \to \mathbb{R}^p\) a differentiable function and \(x\) a point of \(V\). Then, for any \(\mu \in \mathbb{R}^p\) it exists a unique vector \(Dg^T(x)(\mu) \in V\) such as

-
-
-
-\[\begin{equation*} - \forall \mu \in \mathbb{R}^p \quad \forall \phi \in V \qquad \left\langle Dg^T(x)(\mu), \phi \right\rangle_{V} = \mu^T Dg(x)(\phi). -\end{equation*}\] -
-
-
-

The linear operator \(Dg^T(x) : \mathbb{R}^p \to V\) is called the transpose of \(Dg(x)\).

-
-
-
-
-

For equality constrained problems (see Definition 3.3 in [feppon_shape_2019]) , null space step \(\xi_J\) and range space step \(\xi_C\) are defined as the following.

-
-
-
-
Definition : Null space and range space directions
-
-

For any domain \(\Omega\), the null space and range space directions \(\xi_C(\Omega) : \mathbb{R}^N \to \mathbb{R}^N\) and \(\xi_J(\Omega) : \mathbb{R}^N \to \mathbb{R}^N\) are defined by

-
-
-
-\[\begin{align*} - \xi_C(\Omega)(X)= &Dg^T(\Omega)\left[(Dg(\Omega) Dg^T(\Omega))^{-1}g(\Omega)\right](X), \\ - \xi_J(\Omega)(X)= &\nabla J(\Omega)(X) - Dg^T(\Omega)\left[(Dg(\Omega)Dg^T(\Omega))^{-1}Dg(\Omega)\nabla J(\Omega)\right](X). -\end{align*}\] -
-
-
-
-
-

In order to control the step size \(\|\theta_n\|_{L^{\infty}(\mathbb{R}^N, \mathbb{R}^N)}\) and keep all values of the displacement of the order of the mesh size, the parameters \(\alpha_{J,n}\) and \(\alpha_{C,n}\) are updated dynamically.

-
-
-
-
Definition : Parameters of directions
-
-

We consider \(A_J>0\) and \(A_C>0\) two parameters and \(n_0>0\) the \(n_0\)-th iteration. The coefficients are updated at every iteration according to the following rules :

-
-
-
-\[ \begin{align*} - \alpha_{J,n} &= \begin{cases} - \frac{A_J \texttt{hmin}}{\|\xi_J(\Omega_n)\|_{L^{\infty}}} \qquad \qquad \qquad \qquad \quad ~ n<n_0,\\ - \frac{A_J \texttt{hmin}}{\max\left\{\|\xi_J(\Omega_n)\|_{L^{\infty}},\|\xi_J(\Omega_{n_0})\|_{L^{\infty}}\right\}} \qquad n\geq n_0, - \end{cases} \\ - \alpha_{C,n} &= \min\left\{0.9, \frac{A_C \texttt{hmin}}{\max\{\texttt{1e-9}, \|\xi_C(\Omega_n)\|_{L^{\infty}}\}}\right\}. - \end{align*}\] -
-
-
-
-
-

As explained in [feppon_shape_2019], \(A_J\) and \(A_C\) don’t need fine tuning. Thus, to obtain a more general algorithm with the least possible parameters, we take \(A_J=A_C=1\).

-
-
-
-
Proposition : Bounded directions
-
-

These normalizations ensure that the infinite norm of the various terms is less than the smallest mesh size, i.e.

-
-
-
-\[\begin{equation*} - \forall n\geq 0 \quad \|\alpha_{J,n}\xi_J(\Omega_n)\|_{L^{\infty}}\leq A_J\texttt{hmin} \quad \text{and} \quad \|\alpha_{C,n}\xi_C(\Omega_n)\|_{L^{\infty}}\leq \min\{\texttt{0.9}, A_C\texttt{hmin}\}. -\end{equation*}\] -
-
-
-
-
-Proof -
-
-

Trivial.

-
-
-
-
-

The key point of this method is to be able to determine the transpose of the differential of the equality constraint. To do this, consider a Hilbert space \(V=H^1(\mathbb{R}^N)\) such that \(V\subset W^{1,\infty}(\Omega, \mathbb{R}^N)\) with the following scalar product

-
-
-
-\[\begin{equation*} - \forall (\theta,\theta') \in V^2, \quad \left\langle \theta, \theta'\right\rangle_{V}=\int_{\Omega}\gamma^2 \nabla \theta:\nabla \theta' + \theta \cdot \theta' -\end{equation*}\] -
-
-
-

where the parameter \(\gamma\) is set proportional to the minimum mesh element size.

-
-
-

The differential of the cost function is assumed to be of the previous form with \(G(\Omega)\) the shape gradient. Thus, as described in [feppon_shape_2019] Chapter 1 on page 54, the gradient of the cost function, \(\nabla J\), is a regularised extension of \(G(\Omega)n\), i.e

-
-
-
-\[\begin{equation}\label{eq:gradJ} - \nabla J(\Omega) = G(\Omega)n \quad \text{on } \Gamma, -\end{equation}\] -
-
-
-

with \(n\) the external unit normal.

-
-
-
-
Proposition : Transposition of the differential for volume conservation
-
-

The transposition of Proposition : Differential of the volume constraint is given by

-
-
-
-\[ \begin{equation}\label{eq:DgTpb} - \begin{cases} --\gamma^2\Delta(Dg^T(\Omega)(e_1)) + Dg^T(\Omega)(e_1) = 0 &\qquad \text{in } \Omega, \\ -\frac{\partial Dg^T(\Omega)(e_1)}{\partial n} = \frac{1}{\gamma^2}n &\qquad \text{on } \Gamma,\\ -\frac{\partial Dg^T(\Omega)(e_1)}{\partial n} = 0 &\qquad \text{on } \Gamma_{fixed}. - \end{cases} - \end{equation}\] -
-
-
-
-
-Proof -
-
-

Let’s determine \(Dg^T\). First at all, let \(e={e_1}\) be the canonical basis of \(\mathbb{R}\). Then, by linearity of differential, we have

-
-
-
-\[\begin{equation*} - \forall \mu \in \mathbb{R} \quad Dg^T(\Omega)(\mu) =\mu_1 Dg^T(\Omega)(e_1). -\end{equation*}\] -
-
-
-

with \(\mu=\mu_1 e_1\). Then, by taking \(\mu=e_1\), for all \(\phi \in V\) and by using Definition : Transpose, we obtain

-
-
-
-\[\begin{equation*} - \left\langle Dg^T(\Omega)(e_1), \phi \right\rangle_{V} = e_1 Dg(\Omega)(\phi), -\end{equation*}\] -
-
-
-

in other words,

-
-
-
-\[\begin{equation*} - \int_{\Omega}\gamma^2 \nabla \left(Dg^T(\Omega)(e_1)\right):\nabla \phi + Dg^T(\Omega)(e_1)\cdot \phi = \int_{\Gamma}\phi \cdot n. -\end{equation*}\] -
-
-
-

We want to solve

-
-
-
-\[\begin{equation*} - \int_{\Omega}\gamma^2 \nabla U : \nabla \phi + U \cdot \phi = \int_{\Gamma}\phi \cdot n, -\end{equation*}\] -
-
-
-

for all \(\phi\) in \(V\) with \(U=Dg^T(\Omega)(e_1)\). By using a integration of parts formula, it follows that

-
-
-
-\[\begin{equation*} - \int_{\Omega}\left(-\gamma^2 \Delta U + U \right)\cdot \phi = \int_{\Gamma}\left(I-\gamma^2 \nabla U\right)n \cdot \phi -\int_{\Gamma_{fixed}}\gamma^2\nabla U n \cdot \phi. -\end{equation*}\] -
-
-
-

By taking \(\nabla U n = I\frac{n}{\gamma^2}\) on \(\Gamma\) and \(\nabla U n = 0\) on \(\Gamma_{fixed}\), we therefore need to solve the following problem

-
-
-
-\[\begin{equation} - \begin{cases} --\gamma^2\Delta(Dg^T(\Omega)(e_1)) + Dg^T(\Omega)(e_1) = 0 &\qquad \text{in } \Omega, \\ -\frac{\partial Dg^T(\Omega)(e_1)}{\partial n} = \frac{1}{\gamma^2}n &\qquad \text{on } \Gamma,\\ -\frac{\partial Dg^T(\Omega)(e_1)}{\partial n} = 0 &\qquad \text{on } \Gamma_{fixed}. - \end{cases} -\end{equation}\] -
-
-
-
-
-
-
Proposition : Null space and range space directions for volume conservation
-
-

The null space and range space directions, in the case of Definition : Set of admissible domains, can be written as follow

-
-
-
-\[\begin{equation} -\xi_J(\Omega)(X)=\nabla J(\Omega)(X) -\frac{\left(\int_{\partial \Omega} \nabla J(\Omega) \cdot n\right)Dg^T(e_1)(X)}{\int_{\Omega}Dg^T(\Omega)\left(e_1\right)\cdot n}, -\end{equation}\] -
-
-
-

and

-
-
-
-\[\begin{equation} -\xi_C(\Omega)(X)=\left(\frac{g(\Omega)}{\int_{\partial \Omega}Dg^T(\Omega)(e_1)\cdot n}\right)Dg^T(\Omega)\left(e_1\right)(X). -\end{equation}\] -
-
-
-
-
-Proof -
-
-

By definition of \(g\) and \(Dg\), we have

-
-
-
-\[\begin{equation*} (Dg(\Omega)Dg^T(\Omega))^{-1}g(\Omega)=\frac{g(\Omega)}{\int_{\partial \Omega}Dg^T(\Omega)(e_1)\cdot n} -\end{equation*}\] -
-
-
-

and

-
-
-
-\[\begin{equation*} -(Dg(\Omega)Dg^T(\Omega))^{-1}Dg(\Omega)\nabla J(\Omega)=\frac{\left(\int_{\partial \Omega} \nabla J(\Omega) \cdot n\right)}{\int_{\Omega}Dg^T(\Omega)\left(e_1\right)\cdot n}. -\end{equation*}\] -
-
-
-

Thus, the result can be deduced.

-
-
-
-
- - - - - -
- - -
-

All the steps above are valid when there are several equality constraints, i.e. \(g : V \to \mathbb{R}^p\) with \(p\geq 1\). Furthermore, it is possible to extend this method to inequality constraints as explained in [feppon_shape_2019].

-
-
-
-
-

Numerical Implementation : For the numerical implementation, a domain sequence, \(\Omega_k\), is computed which verifies the following constraints

-
-
-
-\[\begin{equation*} - \partial \Omega_k = \Gamma_k \cup \Gamma_{fixed}, -\end{equation*}\] -
-
-
-

and

-
-
-
-\[\begin{equation*} -\Omega_{k+1} = (\textrm{Id} + \theta_k)(\Omega_k). -\end{equation*}\] -
-
-
-
-
Expansion problem for the NSGF method
-
-

Using the NSGF method, we define the displacement field \(\theta_k\) over all \(\Omega_k\) such that

-
-
-
-\[\begin{equation*} -\begin{cases} --\Delta\theta_k = 0 &\qquad \text{in } \Omega_k \\ -\theta_k = 0 &\qquad \text{on } \Gamma_{fixed} \\ -\frac{\partial \theta_k}{\partial n} = -\alpha_C \xi_C(\Omega_k) - \alpha_J \xi_J(\Omega_k) &\qquad \text{on } \Gamma_k. - \end{cases} -\end{equation*}\] -
-
-
-
-
-
-
-
-

3. Results

-
-
- - - - - -
- - -
-

It is important to note that any other shape could have been selected for the initial domain as long as it shared the same fixed boundaries.

-
-
-
-
-

3.1. Linear-elasticity : the cantilever

-
-

We consider a homogeneous isotropic elastic solid which occupies a bounded domain \(\Omega\) in \(\mathbb{R}^2\). We suppose that the boundary is composed of three parts

-
-
-
-\[\begin{equation*} -\partial \Omega = \Gamma \cup \Gamma_N \cup \Gamma_D -\end{equation*}\] -
-
-
-

where \(\Gamma\) is a variable part traction-free (homogeneous Neumann), \(\Gamma_D\) is a fixed part on which the solid is fixed (homogeneous Dirichlet) and \(\Gamma_N\) is a fixed part on which \(f\) forces are applied (inhomogeneous Neumann). The cantilever problem is illustrated in the following image.

-
-
-
-cantilever -
-
Figure 2. Illustration of the 2D cantilever with a hole. The \(\Gamma_D\) and \(\Gamma_N\) boundaries (Dirichlet condition and Neumann condition) are fixed. The \(\Gamma\) boundary (of which the hole is a part) is movable in order to optimize the shape.
-
-
-
-
Cantilever PDE
-
-

The displacement \(u : \mathbb{R}^N \to \mathbb{R}^N\) is the solution of the system

-
-
-
-\[\begin{equation} -\begin{cases} --\nabla \cdot \sigma = 0 &\qquad \text{in } \Omega,\\ -\sigma = 2\mu e(u) + \lambda tr(e(u))I &\qquad \text{in } \Omega,\\ -u=0 &\qquad \text{on } \Gamma_D,\\ -\sigma n = f &\qquad \text{on } \Gamma_N,\\ -\sigma n = 0 &\qquad \text{on } \Gamma, - \end{cases} -\end{equation}\] -
-
-
-

with \(\sigma(u)\in M_N(\mathbb{R})\) the stress tensor, \(\lambda, \mu >0\) the Lamé coefficients of the material and \(e(u)=\frac{1}{2}\left(\nabla u + \nabla u^T\right)\) the deformation tensor.

-
-
-
-
-
-
Cantilever cost function
-
-

We wish to solve the following minimisation problem

-
-
-
-\[\begin{equation}\label{cantilever:minpb} -\inf_{\Omega \in \Omega_{ad}} \left\{ J(\Omega) = \int_{\Gamma_N} f \cdot u \right\} -\end{equation}\] -
-
-
-

where the set of admissible forms is defined as

-
-
-
-\[\begin{equation}\label{cantilever:thetaad} -\Omega_{ad}=\left\{\Omega \in C(\Omega_0) \mid \Gamma_D \cup \Gamma_N \subset \partial \Omega, ~ g(\Omega)=0\right\} -\end{equation}\] -
-
-
-

with \(g\) defined by Definition : Set of admissible domains.

-
-
-
-
- - - - - -
- - -
-

It should be noted that the domain of integration of the cost function does not depend on \(\Omega\) because \(\Gamma_N\) is fixed. Thus, the dependence with respect to the domain is ensured only through the solution of the previous model.

-
-
-
-
-

3.1.1. Theoretical formulation of the problem

-
-
-
The weak formulation of Cantilever PDE
-
-

The weak formulatio of the cantilever problem is given by

-
-
-
-\[\begin{equation} -\int_{\Omega}2\mu e(u):e(\phi) + \int_{\Omega}\lambda (\nabla \cdot u)(\nabla \cdot \phi) = \int_{\Gamma_N}f \cdot \phi. -\end{equation}\] -
-
-
-
-
-Proof -
-
-

The divergence of \(\sigma\) is given by

-
-
-
-\[ \begin{equation*} -\nabla \cdot \sigma = \left(\sum_{j=1}^N \frac{\partial \sigma_{ij}}{\partial x_j}\right)_{1\leq i\leq N}. -\end{equation*}\] -
-
-
-

Let \(u \in H^1_{\Gamma_N}(\Omega)\) solution of Cantilever PDE. By using the fact that \(tr(e(u))=\nabla \cdot u\), the equation can be written for \(1\leq i\leq N\) as

-
-
-
-\[\begin{equation*} - -\sum_{j=1}^N \frac{\partial}{\partial x_j}\left[\mu\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)+\lambda (\nabla \cdot u)\delta_{ij}\right]=0. -\end{equation*}\] -
-
-
-

By multiplying by a test function \(\phi \in H^1_{\Gamma_D}(\Omega)\) and integrating over the domain \(\Omega\), we have

-
-
-
-\[\begin{equation*} - \int_{\Omega}\sum_{j=1}^N\mu\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_j}\right)\frac{\partial \phi_i}{\partial x_j}+\int_{\Omega}\lambda (\nabla \cdot u) \frac{\partial \phi_i}{\partial x_i}=\int_{\Gamma_N}f_i \phi_i. -\end{equation*}\] -
-
-
-

In addition, we have the following result

-
-
-
-\[\begin{equation}\label{ipp:utils1} -\sum_{i,j=1}^{N}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\frac{\partial \phi_i}{\partial x_j}=\frac{1}{2}\sum_{i,j=1}^{N}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\left(\frac{\partial \phi_i}{\partial x_j}+\frac{\partial \phi_j}{\partial x_j}\right)=2e(u):e(\phi) -\end{equation}\] -
-
-
-

Finally, summing over the index \(i\) in the previous formulation and using this result, we obtain the result.

-
-
-
-
-
-
Definition : Lagrangian of Cantilever PDE and Cantilever cost function
-
-

We denote the Lagrangian of this problem by

-
-
-
-\[\begin{equation*} -\begin{array}{rcl} -\mathcal{L}:\Omega_{ad}\times \left(H^1_{\Gamma_D}(\mathbb{R}^N)\right)^2&\to& \mathbb{R}\\ -(\Omega, v,q) &\mapsto &\int_{\Omega}2\mu e(v):e(q) + \int_{\Omega} \lambda (\nabla \cdot v)(\nabla \cdot q) - \int_{\Gamma_N} f \cdot q + \int_{\Gamma_N} f \cdot v. -\end{array} -\end{equation*}\] -
-
-
-
-
- - - - - -
- - -
-

It should be noted that the different variables of \(\mathcal{L}\) are independent because we consider Lagrange multipliers on \(H^1_{\Gamma_D}(\mathbb{R}^N)\) and not \(H^1_{\Gamma_D}(\Omega)\) and especially that \(\Gamma_N\) is independent of \(\Omega\) because it is fixed.

-
-
-
-
-
-
Dual PDE of the cantilever problem
-
-

The dual PDE of the cantilever problem is

-
-
-
-\[\begin{equation*} -\begin{cases} --\nabla \cdot \sigma = 0 &\qquad \text{on } \Omega\\ -\sigma = 2\mu e(p) + \lambda tr(e(p))I &\qquad \text{on } \Omega\\ -p=0 &\qquad \text{on } \Gamma_D\\ -\sigma n = -f &\qquad \text{on } \Gamma_N\\ -\sigma n = 0 &\qquad \text{on } \Gamma - \end{cases} -\end{equation*}\] -
-
-
-
-
-Proof -
-
-

The partial derivative with respect to \(v\) : let \(\phi \in H^1_{\Gamma_D}(\mathbb{R}^N)\)

-
-
-
-\[\begin{equation*} -\left\langle \frac{\partial \mathcal{L}}{\partial v}(\Omega, v, q), \phi \right\rangle = \int_{\Omega}2\mu e(\phi):e(q) + \int_{\Omega} \lambda (\nabla \cdot \phi) (\nabla \cdot q) + \int_{\Gamma_N} f \cdot \phi, -\end{equation*}\] -
-
-
-

which, when it vanishes, corresponds to the weak formulation of the dual problem.

-
-
-
-
- - - - - -
- - -
-

Note that this corresponds exactly to the primal problem with a minus sign. This is due in particular to the choice of \(J\) and the boundary conditions of the primal problem. Indeed, we have \(J\) which depends on \(f\) on \(\Gamma_N\) and the boundary conditions are all null except on \(\Gamma_N\) which is equal to \(f\). Thus, we can avoid solving the dual problem in this case and simply use the solution of the primal problem.

-
-
-
-
-
-
Shape gradient for the cantilever problem
-
-

The shape gradient for the cantiler problem is defined by

-
-
-
-\[\begin{equation}\label{cantilever:gradshape} - G(\Omega)=-2\mu\|e(u)\|^2-\lambda(tr(e(u)))^2. -\end{equation}\] -
-
-
-
-
-Proof -
-
-

Let’s calculate the partial derivative of \(\mathcal{L}\) with respect to the \(\Omega\) domain in the \(\theta\) direction

-
-
-
-\[\begin{equation*} - \frac{\partial \mathcal{L}}{\partial \Omega}(\Omega, v, q)(\theta)= \int_{\partial \Omega}\theta \cdot n \left[ 2\mu e(v):e(q) + \lambda (\nabla \cdot v)(\nabla \cdot q)\right] -\end{equation*}\] -
-
-
-

by using Proposition 1. When we evaluate this derivative with the state \(u(\Omega)\) and the adjoint state \(p(\Omega)=u(\Omega)\), we find exactly the value of the derivative of the cost function

-
-
-
-\[\begin{equation}\label{cantilever:shaped} -\frac{\partial \mathcal{L}}{\partial \Omega}\left(\Omega, u(\Omega), p(\Omega)\right)(\theta)= -\int_{\Gamma}\theta \cdot n \left[ 2\mu \|e(u)\|^2 + \lambda tr(e(u))^2\right] = DJ(\Omega)(\theta), -\end{equation}\] -
-
-
-

as explained here Proposition : Differential of the cost fucntion thanks to the Lagrangian. Thus by definition we obtain the shape gradient.

-
-
-
-
-
-

3.1.2. Experimental Evaluation

-
-

In this section, we will present various results on the optimization problem concerning the shape of cantilevers. The initial domain considered is a trapezoid with two feet in the 2D case (four feet in the 3D case). The choice of a trapezoid as the initial shape is advantageous due to its simplicity and its ability to approximate the solution of the problem. Throughout this section, unless specified otherwise, all results have been obtained using \(\mathbb{P}_1\) continuous finite elements. All simulations are carried out using the classic gradient descent method. The NSGF method is used in the following example, a rigid body in a Stokes fluid.

-
-
-

2D simulations for various types of cantilever :

-
-
-

We present the application of shape optimization to a 2D cantilever with two pillars. The specific parameter values used in the analysis are detailed in the following table.

-
- - ---- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Table 1. Geometric parameters and Lamé coefficient of the initial domain for the 2D case of the cantilever. The Lamé coefficients are \(\lambda\) and \(\mu\). \(H\) represents the height of the trapezoid, \(L_1\) the large base and \(L_2\) the small base. The force applied to the curve \(\Gamma_N\) is defined by \(f\).

Symbol

Value (dimensionless)

\(\lambda\)

\(50/9\)

\(\mu\)

\(350/27\)

\(H\)

\(9\)

\(L_1\)

\(8\)

\(L_2\)

\(2\)

\(f\)

\((0,-1)\)

-
-

No hole :

-
-
-

We begin by considering the simplest case, which is the cantilever without any holes, as depicted in Figures below. To ensure clear visibility of the mesh in print, a discretization parameter of \(h=0.4\) is set. For the initialization of the Lagrange multiplier, we choose \(l=0.5\). Additionally, the values of \(a=b=0.5\) and \(c=10\) are set, along with a descent step of \(t=0.2\). These specific values have been determined through empirical testing to achieve optimal performance.

-
- - ------- - - - - - - - - - - - - - - - - -
Table 2. Comparison of the cost function between the initial domain and the final domain, the volume error between the two domains, the \(H^1\)-norm of the displacement field at the end and the number of iterations for the 2D cantilever without hole.

\(J(\Omega_0)\)

\(J(\Omega_{n_{final}})\)

\(||\Omega_0|-|\Omega_{n_{final}}||\)

\(||\theta_{n_{final}}||_{H^1}\)

\(n_{final}\)

\(4.6834\)

\(3.1346\)

\(1.3630e-2\)

\(1.696e-2\)

\(125\)

-
-
-results nohole 2D l0.5 t0.2 h0.4 a0.5b0.5 c10 n200 -
-
Figure 3. Evolution of the cost function, the volume of the domain \(\Omega_k\) and the norm of \(\theta_k\) of the 2D cantilever without hole with the following parameters: \(h=0.4\), \(t=0.2\), \(l=0.5\), \(a=b=0.5\) and \(c=10\). The iterations range from stem[0] to \(125\).
-
- - --- - - - - - -
Table 3. Visualisation of the results obtained during different iterations of the algorithm for the shape optimisation of the 2D cantilever without hole with the following parameters: \(h=0.4\), \(t=0.2\), \(l=0.5\), \(a=b=0.5\) and \(c=10\). The initial domain (\(k=0\)) is shown on the left. The intermediate domain (\(k=60\)) is displayed in the middle. The final domain (\(k=125\)) is displayed on the right.

nohole 0 -nohole 60 -nohole 125

-
-

One hole :

-
-
-

The second test consists of studying the cantilever with a single hole. The results obtained are presented in Figures below. We use a discretization parameter of \(h=0.4\). For the Lagrange multiplier, we initialize with \(l=0.5\), and set \(a=b=0.5\) and \(c=10\) with a descent step of \(t=0.2\).

-
- - ------- - - - - - - - - - - - - - - - - -
Table 4. Comparison of the cost function between the initial domain and the final domain, the volume error between the two domains, the \(H^1\)-norm of the displacement field at the end and the number of iterations for the one hole 2D cantilever.

\(J(\Omega_0)\)

\(J(\Omega_{n_{final}})\)

\(||\Omega_0|-|\Omega_{n_{final}}||\)

\(||\theta_{n_{final}}||_{H^1}\)

\(n_{final}\)

\(4.7035\)

\(3.2162\)

\(2.9980e-3\)

\(1.8906e-2\)

\(125\)

-
-
-results hole 2D l0.5 t0.2 h0.4 a0.5b0.5 c10 n200 -
-
Figure 4. Evolution of the cost function, the volume of the domain \(\Omega_k\) and the norm of \(\theta_k\) of the 2D cantilever with one hole with the following parameters: \(h=0.\), \(t=0.2\), \(l=0.5\), \(a=b=0.5\) and \(c=10\). The iterations range from \(0\) to \(125\).
-
- - --- - - - - - -
Table 5. Visualisation of the results obtained during different iterations of the algorithm for the shape optimisation of the 2D cantilever with one hole with the following parameters: \(h=0.4\), \(t=0.2\), \(l=0.5\), \(a=b=0.5\) and \(c=10\). The initial domain (\(k=0\)) is shown on the left. The intermediate domain (\(k=60\)) is displayed in the middle. The final domain (\(k=125\)) is displayed on the right.

hole 0 -hole 60 -hole 125

-
-

Four holes :

-
-
-

The last type of cantilever studied is one with four holes in its domain. The results are displayed in the figures below. We use a discretization parameter of \(h=0.4\). For the Lagrange multiplier, we initialize with \(l=0.5\), and set \(a=b=0.5\) and \(c=10\) with a descent step of \(t=0.2\).

-
- - ------- - - - - - - - - - - - - - - - - -
Table 6. Comparison of the cost function between the initial domain and the final domain, the volume error between the two domains, the \(H^1\)-norm of the displacement field at the end and the number of iterations for the four holes 2D cantilever.

\(J(\Omega_0)\)

\(J(\Omega_{n_{final}})\)

\(||\Omega_0|-|\Omega_{n_{final}}||\)

\(||\theta_{n_{final}}||_{H^1}\)

\(n_{final}\)

\(5.1589\)

\(3.3770\)

\(2.5317e-3\)

\(1.5615e-2\)

\(125\)

-
-
-results 4holes 2D l0.5 t0.2 h0.4 a0.5b0.5 c10 n200 -
-
Figure 5. Evolution of the cost function, the volume of the domain \(\Omega_k\) and the norm of \(\theta_k\) of the 2D cantilever with 4 holes with the following parameters: \(h=0.4\), \(t=0.2\), \(l=0.5\), \(a=b=0.5\) and \(c=10\). The iterations range from \(0\) to \(125\).
-
- - --- - - - - - -
Table 7. Visualisation of the results obtained during different iterations of the algorithm for the shape optimisation of the 2D cantilever with 4 holes with the following parameters: \(h=0.4\), \(t=0.2\), \(l=0.5\), \(a=b=0.5\) and \(c=10\). The initial domain (\(k=0\)) is shown on the left. The intermediate domain (\(k=60\)) is displayed in the middle. The final domain (\(k=125\)) is displayed on the right.

4hole 0 -4hole 60 -4hole 125

-
-

The results presented above demonstrate that satisfactory outcomes can be achieved using the proposed method. Notably, the obtained shapes closely resemble those reported in various papers that focus on optimizing the geometrical shape of cantilevers, such as [allaire_conception_2006]. Furthermore, we successfully decrease the cost function while approaching the initial volume of the domain, which aligns precisely with the desired behavior. It is worth noting that potential geometric instabilities, in the form of spikes, may emerge on \(\Gamma_D\) after a certain number of iterations. These spikes are also observed in the 3D case, as illustrated in the subsequent figures.

-
-
-

3D simulations for various types of cantilever :

-
-
-

In this section, we address the case of the 3D cantilever, which introduces additional complexities compared to the 2D case and results in longer computation times. Incorporating holes into the structure also increases the risk of encountering mesh superposition issues. The overall geometry of the 3D cantilever remains similar to the 2D case, with the addition of four pillars, and boundary conditions are now applied to surfaces instead of curves. The specific parameter values used in the analysis are provided in detail in following table.

-
- - ---- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Table 8. Geometric parameters and Lamé coefficient of the initial domain for the 2D case of the cantilever. The Lamé coefficients are \(\lambda\) and \(\mu\). \(H\) represents the height of the truncated pyramid, \(C_1\) the side of the largest square base and \(C_2\) the side of the smallest square base. The force applied to the surface \(\Gamma_N\) is defined by \(f\).

Symbol

Value (dimensionless)

\(\lambda\)

\(50/9\)

\(\mu\)

\(350/27\)

\(H\)

\(9\)

\(C_1\)

\(8\)

\(C_2\)

\(2\)

\(f\)

\((0,-1,0)\)

-
-

No hole :

-
-
-

Let’s first consider the case without a hole. We use a discretization parameter of \(h=0.8\). For the coefficients affecting the optimization problem, we set \(l=0.5\), \(a=b=0.5\), \(c=3\), and choose a descent step \(t\) of \(0.1\). The results of this test are presented in Figures below.

-
- - ------- - - - - - - - - - - - - - - - - -
Table 9. Comparison of the cost function between the initial domain and the final domain, the volume error between the two domains, the \(H^1\)-norm of the displacement field at the end and the number of iterations for the no hole 3D cantilever.

\(J(\Omega_0)\)

\(J(\Omega_{n_{final}})\)

\(||\Omega_0|-|\Omega_{n_{final}}||\)

\(||\theta_{n_{final}}||_{H^1}\)

\(n_{final}\)

\(4.8060\)

\(3.0484\)

\(0.1439\)

\(2.1647e-2\)

\(400\)

-
-
-results nohole 3D l0.5 t0.1 h0.8 a0.5b0.5 c2 n500 -
-
Figure 6. Evolution of the cost function, the volume of the domain \(\Omega_k\) and the norm of \(\theta_k\) of the 3D cantilever without hole with the following parameters: \(h=0.8\), \(t=0.1\), \(l=0.5\), \(a=b=0.5\) and \(c=2\). The iterations range from \(0\) to \(400\).
-
- - --- - - - - - -
Table 10. Visualisation of the results obtained during different iterations of the algorithm for the shape optimisation of the 3D cantilever without hole with the following parameters: \(h=0.8\), \(t=0.1\), \(l=0.5\), \(a=b=0.5\) and \(c=2\). The initial domain (\(k=0\)) is shown on the left. The intermediate domain (\(k=150\)) is displayed in the middle. The final domain (\(k=400\)) is displayed on the right.

3Dnohole 0 -3Dnohole 120 -3Dnohole 400

-
-

One hole : -In the next test, we add a hole inside the material. The results are presented in Figures below. We use a discretization parameter of \(h=0.8\). The other coefficients used are \(l=0.5\), \(a=b=0.5\), \(c=2\), and \(t=0.08\).

-
- - ------- - - - - - - - - - - - - - - - - -
Table 11. Comparison of the cost function between the initial domain and the final domain, the volume error between the two domains, the \(H^1\)-norm of the displacement field at the end and the number of iterations for the one hole 3D cantilever.

\(J(\Omega_0)\)

\(J(\Omega_{n_{final}})\)

\(||\Omega_0|-|\Omega_{n_{final}}||\)

\(||\theta_{n_{final}}||_{H^1}\)

\(n_{final}\)

\(6.4813\)

\(3.7726\)

\(0.2582\)

\(4.7253e-2\)

\(250\)

-
-
-results hole 3D l0.5 t0.08 h0.8 a0.5b0.5 c2 n500 -
-
Figure 7. Evolution of the cost function, the volume of the domain \(\Omega_k\) and the norm of \(\theta_k\) of the 3D cantilever with holes with the following parameters: \(h=0.8\), \(t=0.08\), \(l=0.5\), \(a=b=0.5\) and \(c=2\). The iterations range from \(0\) to \(250\).
-
- - --- - - - - - -
Table 12. Visualisation of the results obtained during different iterations of the algorithm for the shape optimisation of the 3D cantilever with holes with the following parameters: \(h=0.8\), \(t=0.08\), \(l=0.5\), \(a=b=0.5\) and \(c=2\). The initial domain (\(k=0\)) is shown on the left. The intermediate domain (\(k=150\)) is displayed in the middle. The final domain (\(k=250\)) is displayed on the right.

3Dhole 0 -3Dhole 120 -3Dhole 250

-
-

It is important to note that the chosen optimal shape in the simulations effectively reduces the cost function while maintaining the volume. However, due to limitations, we had to prematurely halt the simulations. If allowed to continue, the cost function would have been further minimized, resulting in a more optimal shape. Nonetheless, the obtained shape presents several issues. Specifically, certain areas such as the feet and top part exhibit more prominent outgrowths, which can be attributed to volume conservation. The shape optimization process tends to hollow out the cantilever below the desired volume (between the four legs). As a result, to compensate for this volume loss, protrusions seem to emerge after a certain number of iterations. Another significant issue encountered is mesh collision and overlap. Additionally, the discontinuities observed in the curves of the different simulations correspond to the remeshing that occurs during the iterations.

-
-
-
-
-

3.2. A rigid body in a Stokes flow

-
-

We consider a rigid object \(S\) set in motion into an incompressible fluid with viscosity \(\mu\) at low Reynolds number. The fluid occupies a bounded domain \(\Omega\) in \(\mathbb{R}^N\). We suppose that the boundary is composed of two parts

-
-
-
-\[\begin{equation*} -\partial \Omega = \Gamma_S \cup \Gamma_B, -\end{equation*}\] -
-
-
-

where \(\Gamma_S\) is a variable part associated with the boundary of the rigid object and \(\Gamma_B\) is a fixed part corresponding to the boundary of the domain containing the fluid. Normally, we would consider an edge at an infinite distance from the object. However, in finite elements we are limited for this kind of hypothesis. It is impossible to consider an infinite wall, and the further away the edge of the domain containing the fluid is, the more elements are needed to mesh it, which has a huge impact on calculation time, particularly in 3D. We will therefore study shape optimisation in a context that is slightly different from [moreau_shapes_2022]. The Stokes problem is illustrated in the following figure.

-
-
-
-stokes -
-
Figure 8. Illustration of the stokes problem. The \(\Gamma_B\) boundaries is fixed. The \(\Gamma_S\) boundary is movable in order to optimize the shape.
-
-
-

Taking the same conditions as in [moreau_shapes_2022], we consider the folling PDE.

-
-
-
-
Rigid object in Stokes flow PDE
-
-

Let a linear background flow \(U^{\infty}\) and \(U\) the translational velocities of the object. The fluid velocity field \(u : \mathbb{R}^N \to \mathbb{R}^N\) and the pressure field \(p : \mathbb{R}^N \to \mathbb{R}\) is the solution of the Stokes equations

-
-
-
-\[\begin{equation} -\begin{cases} --\mu \Delta u + \nabla p = 0 &\qquad \text{in } \Omega,\\ -\nabla \cdot u = 0 &\qquad \text{in } \Omega,\\ -u = U &\qquad \text{on } \Gamma_S,\\ -u = U^{\infty} &\qquad \text{on } \Gamma_B. - \end{cases} -\end{equation}\] -
-
-
-
-
-
-
Rigid object in Stokes flow cost function
-
-

We wish to solve the following minimisation problem

-
-
-
-\[\begin{equation} -\inf_{\Omega \in \Omega_{ad}} \left\{ J_{\alpha}(\Omega) = \int_{\Gamma_S} \sigma(u,p)n\cdot \alpha \right\} -\end{equation}\] -
-
-
-

where \(n\) is the normal to \(\Gamma_S\) pointing outwards with respect to the fluid domain (therefore inward w.r.t the object) and \(\sigma\) is the stress tensor defined as

-
-
-
-\[\begin{equation*} - \sigma(u,p)=-pI + 2\mu e(u) -\end{equation*}\] -
-
-
-

where \(e\) is the deformation tensor, given by

-
-
-
-\[\begin{equation*} - e(u)=\frac{1}{2}\left(\nabla u + \nabla u^T\right). -\end{equation*}\] -
-
-
-

We define the set of admissible forms as

-
-
-
-\[\begin{equation} -\Omega_{ad}=\left\{\Omega \in C(\Omega_0) \mid \Gamma_B \subset \partial \Omega, ~ g(\Omega)=0\right\} -\end{equation}\] -
-
-
-

where \(g\) is defined as in Definition : Set of admissible domains.

-
-
-
-
-

By taking precise values for \(U\), \(\alpha\) and \(U^{\infty}\), one can express the cost function as an input to the large strength tensor explained in [moreau_shapes_2022]. This brings us back to a problem of resistance: what is the optimal shape to have the least possible resistance on the fluid. They correspond with the projection of the hydrodynamic drag force and torque - exerted by the moving particle to the fluid - along one of the basis vectors.

-
- - ------ - - - - - - - - - - - - - - - - - - - - - - - - - - -
Table 13. Cost function \(J_{\alpha}\) associated with the entries of the grand resistance tensor (see [moreau_shapes_2022]) with respect to the choice of \(U\), \(\alpha\) and \(U^{\infty}\).

\(J_{\alpha}\)

\(\alpha\)

\(U\)

\(U^{\infty}\)

\(K_{ij}\)

\(e_j\)

\(e_i\)

\(0\)

\(Q_{ij}\)

\(e_j \wedge (x,y,z)^T\)

\(e_i \wedge (x,y,z)^T\)

\(0\)

\(C_{ij}\)

\(e_j \wedge (x,y,z)^T\)

\(e_i\)

\(0\)

-
-

The problem setup is well presented in the following figure taken from the article [moreau_shapes_2022]. Several examples of resistance problem are illustrated on the right side of the figure for different entries of the grand resistance tensor.

-
-
-
-stokes flow shape privat -
-
Figure 9. Problem setup taken from [moreau_shapes_2022] : according the previous notation, we have \(S=\mathcal{S}\), \(B=\mathcal{B}\), \(\Omega=\mathcal{V}\) and the opposite direction of the normal, i.e \(-n\). The left figure illustrate the Stokes problem in three dimension. In the right side, several examples of resistance problem for different entries of the grand resistance tensor are presented.
-
-
-

3.2.1. Theoretical formulation of the problem

-
-

Let us start by defining the following property which is an important result by integration of parts in the field of fluid mechanics and which will be of great help to us in the future.

-
-
-
-
Integration by parts formula for Stokes
-
-

Let \(v\) and \(\phi\) in \(H^1(\mathbb{R}^N)\), we have

-
-
-
-\[\begin{equation*} - -\int_{\Omega}\left(\Delta v + \nabla(\nabla \cdot v)\right)\cdot \phi =2\int_{\Omega}e(v):e(\phi)-2\int_{\partial \Omega}e(v)n\cdot \phi. -\end{equation*}\] -
-
-
-
-
-Proof -
-
-

Let \(v\) and \(\phi\) in \(H^1(\mathbb{R}^N)\). -First of all, by integration by parts we obtain the following relation

-
-
-
-\[\begin{equation*} - -\int_{\Omega}\sum_{j=1}^N\frac{\partial}{\partial x_j}\left(\frac{\partial v_i}{\partial x_j}+ \frac{\partial v_j}{\partial x_i}\right)\phi_i = \int_{\Omega}\sum_{j=1}^{N}\left(\frac{\partial v_i}{\partial x_j}+\frac{\partial v_j}{\partial x_i}\right)\frac{\partial \phi_i}{\partial x_j}-\int_{\partial \Omega}\sum_{j=1}^N\left(\frac{\partial v_i}{\partial x_j}+\frac{\partial v_j}{\partial x_i}\right)\phi_i n_j -\end{equation*}\] -
-
-
-

By summing over the index \(i\) on the left and right, and using the relation obtain in the end of the proof of Cantilever PDE, we finally have that

-
-
-
-\[\begin{equation*} - -\int_{\Omega}\left(\Delta v + \nabla(\nabla \cdot v)\right)\cdot \phi = 2\int_{\Omega}e(v):e(\phi)-\int_{\partial \Omega} 2e(v)n\cdot \phi. -\end{equation*}\] -
-
-
-
-
-
-
Definition : Lagrangian of Rigid object in Stokes flow PDE and Rigid object in Stokes flow cost function
-
-

Since we are working with Dirichlet conditions, we must use the form of the equation Rigid object in Stokes flow PDE directly for the variational formulation and add to the Lagrangian, Lagrange multipliers associated with Dirichlet constraints. We denote the Lagrangian of this problem by

-
-
-
-\[\begin{equation*} -\begin{array}{rcl} -\mathcal{L}_{\alpha}:\Omega_{ad}\times \left(H^1(\mathbb{R}^N)\right)^4\times \left(L^2(\mathbb{R}^N)\right)^2&\to& \mathbb{R}\\ -(\Omega, v,h1,\lambda, \beta, q,h2) &\mapsto & \int_{\Gamma_S}\sigma(v,q)n\cdot \alpha + \int_{\Omega}\left[-\mu \Delta v + \nabla q\right]\cdot h_1 \\ -&& + \int_{\Omega}(\nabla \cdot v) h_2 + \int_{\Gamma_S}\lambda \cdot (v-U)\\&& + \int_{\Gamma_B}\beta \cdot (v-U^{\infty}) -\end{array} -\end{equation*}\] -
-
-
-
-
-
-
Dual PDE of the rigid object in Stokes flow problem
-
-

The dual PDE of the rigid object in Stokes flow problem is

-
-
-
-\[\begin{equation}\label{stockes:adjpb} -\begin{cases} --\mu \Delta u^{*} + \nabla p^{*} = 0 &\qquad \text{in } \Omega,\\ -\nabla \cdot u^{*} = 0 &\qquad \text{in } \Omega,\\ -u^{*} = \alpha &\qquad \text{on } \Gamma_S,\\ -u^{*} = 0 &\qquad \text{on } \Gamma_B, - \end{cases} -\end{equation}\] -
-
-
-
-
-Proof -
-
-

Let’s start by determining the dual problem. In all the following, for the sake of readability, we will omit the different variables of the Lagrangian. Let’s derive \(\mathcal{L}_{\alpha}\) w.r.t \(q\). For all \(\phi \in H^1(\mathbb{R}^N)\), we have

-
-
-
-\[\begin{equation*} - \left\langle \frac{\partial \mathcal{L}_{\alpha}}{\partial q}, \phi\right\rangle = \int_{\Gamma_S} \sigma(0,\phi)n\cdot \alpha + \int_{\Omega}\nabla \phi \cdot h_1. -\end{equation*}\] -
-
-
-

Using a integration by parts formula, then

-
-
-
-\[\begin{equation*} -\left\langle \frac{\partial \mathcal{L}_{\alpha}}{\partial q}, \phi\right\rangle = -\int_{\Gamma_S} \phi (n\cdot \alpha) - \int_{\Omega}\phi\nabla \cdot h_1 + \int_{\partial \Omega}\phi (n \cdot h_1). -\end{equation*}\] -
-
-
-

Taking \(\phi\) with compact support in \(\Omega\), the variables \(h_1=u^{*}\) and \(h_2=p^{*}\) and using Proposition : Dual problem thanks to the Lagrangian, we finally have

-
-
-
-\[\begin{equation}\label{stokes:adj2} - \nabla\cdot u^{*} = 0 \qquad \text{in } \Omega. -\end{equation}\] -
-
-
-

Let’s derive \(\mathcal{L}_{\alpha}\) w.r.t \(v\). For all \(\phi \in H^1(\mathbb{R}^N)\), then

-
-
-
-\[\begin{equation} - \left\langle \frac{\partial \mathcal{L}_{\alpha}}{\partial v},\phi \right\rangle = \int_{\Gamma_S}\sigma(\phi,0)n\cdot \alpha +\int_{\Omega}-\mu\Delta \phi \cdot h_1 + h_2\nabla \cdot \phi + \int_{\Gamma_S} \lambda \cdot v + \int_{\Gamma_B}\beta \cdot v. -\end{equation}\] -
-
-
-

Using twice integration by parts, we obtain the following formulation

-
-
-
-\[\begin{multline*} - \left\langle \frac{\partial \mathcal{L}_{\alpha}}{\partial v},\phi \right\rangle = -\int_{\Omega}\left[\mu \Delta h_1 + \nabla h_2\right]\cdot \phi + \int_{\partial \Omega}\mu \nabla h_1 n \cdot \phi - \mu \nabla \phi n \cdot h_1 + h_2(n\cdot \phi) \\+ \int_{\Gamma_S}\lambda \cdot v +2\mu e(\phi)n\cdot \alpha + \int_{\Gamma_B}\beta \cdot \phi. -\end{multline*}\] -
-
-
-

Taking \(\phi\) with compact support in \(\Omega\), the variables \(h_1=u^{*}\) and \(h_2=-p^{*}\) and using Proposition : Dual problem thanks to the Lagrangian, -even if it means changing the sign of \(p^{*}\), we finally have that

-
-
-
-\[\begin{equation}\label{stokes:adj1} - -\mu \Delta u^{*} + \nabla p^{*} = 0 \qquad \text{in } \Omega. -\end{equation}\] -
-
-
-

Let us start from the first equation for \(\frac{\partial \mathcal{L}_{\alpha}}{\partial v}\) and apply the Integration by parts formula for Stokes twice in order to obtain the following result

-
-
-
-\[\begin{multline*} -\left\langle \frac{\partial \mathcal{L}_{\alpha}}{\partial v},\phi \right\rangle = -\int_{\Omega}\mu \left[\Delta h_1 + \nabla(\nabla \cdot h_1) +\nabla h_2\right]\cdot \phi + \int_{\Omega}\mu \nabla(\nabla\cdot \phi)\cdot h_1 + \int_{\partial \Omega}h_2(n\cdot \phi)\\+2\mu\int_{\partial \Omega} e(h_1)n\cdot \phi-e(\phi)n\cdot h_1 + \int_{\Gamma_B}\beta \cdot \phi+\int_{\Gamma_S} \lambda\cdot \phi +2\mu e(\phi)n\cdot \alpha. -\end{multline*}\] -
-
-
-

Let us take the test functions \(\phi\) with zero divergence such that \(\phi\) vanishes on \(\partial \Omega\) and \(e(\phi)n\) describes \(L^2(\Gamma_S)\) (respectively \(L^2(\Gamma_B)\)). By taking the variables \(h_1=u^{*}\) and \(h_2=-p^{*}\) solution of [stokes:adj1] and [stokes:adj2], and using Proposition : Dual problem thanks to the Lagrangian, we obtain

-
-
-
-\[\begin{align} - u^{*}=\alpha \qquad \text{on } \Gamma_S, \label{stokes:adj41}\\ - u^{*}= 0 \qquad \text{on } \Gamma_B \label{stokes:adj42}. -\end{align}\] -
-
-
-

Using the previous equations, we find the result.

-
-
-
-
-
-
Lemma 2
-
-

The Lagrange multipliers associated with Dual PDE of the rigid object in Stokes flow problem are

-
-
-
-\[\begin{equation}\label{stockes:lagmult} -\begin{cases} -\lambda^{*} = -2\mu e(u^{*})n-p^{*}n &\qquad \text{on } \Gamma_S,\\ -\beta^{*} = -2\mu e(u^{*})n - p^{*}n &\qquad \text{on } \Gamma_B. - \end{cases} -\end{equation}\] -
-
-
-
-
-Proof -
-
-

Let’s derive \(\mathcal{L}_{\alpha}\) w.r.t \(v\). For all \(\phi \in H^1(\mathbb{R}^N)\), then

-
-
-
-\[\begin{equation} - \left\langle \frac{\partial \mathcal{L}_{\alpha}}{\partial v},\phi \right\rangle = \int_{\Gamma_S}\sigma(\phi,0)n\cdot \alpha +\int_{\Omega}-\mu\Delta \phi \cdot h_1 + h_2\nabla \cdot \phi + \int_{\Gamma_S} \lambda \cdot v + \int_{\Gamma_B}\beta \cdot v. -\end{equation}\] -
-
-
-

By applying the Integration by parts formula for Stokes twice, we have

-
-
-
-\[\begin{multline*} -\left\langle \frac{\partial \mathcal{L}_{\alpha}}{\partial v},\phi \right\rangle = -\int_{\Omega}\mu \left[\Delta h_1 + \nabla(\nabla \cdot h_1) +\nabla h_2\right]\cdot \phi + \int_{\Omega}\mu \nabla(\nabla\cdot \phi)\cdot h_1 + \int_{\partial \Omega}h_2(n\cdot \phi)\\+2\mu\int_{\partial \Omega} e(h_1)n\cdot \phi-e(\phi)n\cdot h_1 + \int_{\Gamma_B}\beta \cdot \phi+\int_{\Gamma_S} \lambda\cdot \phi +2\mu e(\phi)n\cdot \alpha. -\end{multline*}\] -
-
-
-

Let us take the test functions \(\phi\) with zero divergence such that \(e(\phi)n\) vanishes on \(\partial \Omega\) and \(\phi\) describes \(L^2(\Gamma_S)\) (respectively \(L^2(\Gamma_B)\)). By taking the variables \(h_1=u^{*}\) and \(h_2=-p^{*}\) solution of Dual PDE of the rigid object in Stokes flow problem, and using Proposition : Dual problem thanks to the Lagrangian, we obtain

-
-
-
-\[\begin{align} - 2\mu e(u^{*})n + p^{*}n + \lambda = 0 \qquad \text{on } \Gamma_S, \label{stokes:adj31}\\ - 2\mu e(u^{*})n + p^{*}n + \beta = 0 \qquad \text{on } \Gamma_B \label{stokes:adj32}. -\end{align}\] -
-
-
-

Hence the result.

-
-
-
-
-
-
Lemma 3
-
-

Let \(u\) be the solution to problem Rigid object in Stokes flow PDE, then

-
-
-
-\[\begin{equation*} - n\cdot \nabla (u-U)n= \nabla \cdot (u-U)=0. -\end{equation*}\] -
-
-
-

The relation is also valid for \(u=u^{*}\) solution of the problem Dual PDE of the rigid object in Stokes flow problem.

-
-
-
-
-Proof -
-
-

See [[moreau_shapes_2022], proof of Proposition 1] by taking the other convention for the normal.

-
-
-
-
-
-
Lemma 4
-
-

Let \(u\) be the solution to problem Rigid object in Stokes flow PDE and \(u^{*}\) solution to problem Dual PDE of the rigid object in Stokes flow problem, then

-
-
-
-\[\begin{equation*} - e(u^{*})n\cdot \nabla (u-U)n = e(u^{*}):e(u-U). -\end{equation*}\] -
-
-
-

The relation is also valid by exchanging the role of \(u\) and \(u^{*}\), and by taking \(U=\alpha\).

-
-
-
-
-Proof -
-
-

See [[moreau_shapes_2022], proof of Proposition 1] by taking the other convention for the normal.

-
-
-
-
-
-
Lemma 5
-
-

Let \((u^{*},p^{*})\) solution to problem Dual PDE of the rigid object in Stokes flow problem, \(u\) solution to problem Rigid object in Stokes flow PDE and \(f\) such as

-
-
-
-\[\begin{equation} - f(\Omega)=-\int_{\partial \Omega}\sigma(u^{*},\pm p^{*})n\cdot (u-U), -\end{equation}\] -
-
-
-

then we have

-
-
-
-\[\begin{equation*} - Df(\Omega)(\theta)=-2\mu\int_{\Gamma_S}\theta\cdot n \left[e(u^{*}):e(u-U)\right]. -\end{equation*}\] -
-
-
-

The relation is also valid by exchanging the role of \((u,p)\) and \((u^{*},p^{*})\) , and by taking \(U=\alpha\).

-
-
-
-
-Proof -
-
-

To simplify, we will note \(a_{\pm}=\sigma(u^{*},\pm p^{*})n\). Using the fact that \(\Gamma_S\) is fixed and Propostion 2, we obtain

-
-
-
-\[\begin{equation*} - Df(\Omega)(\theta)=-\int_{\Gamma_S}\theta\cdot n \left[\nabla\left(a_{\pm}\cdot (u-U)\right)\cdot n + (\nabla \cdot n)a_{\pm}\cdot (u-U)\right]. -\end{equation*}\] -
-
-
-

Thanks to the boundary conditions, the last term vanishes and we have

-
-
-
-\[\begin{equation*} - Df(\Omega)(\theta)=-\int_{\Gamma_S}\theta\cdot n \left[\nabla\left(a_{\pm}\cdot (u-U)\right)\cdot n\right]. -\end{equation*}\] -
-
-
-

However, we have that \(\nabla \left(a_{\pm}\cdot (u-U)\right)\cdot n = (\nabla a_{\pm})(u-U)\cdot n +\nabla (u-U)a_{\pm} \cdot n\). Thus, injecting this equation into the above integral while neglecting the terms due to the boundary conditions, we obtain

-
-
-
-\[\begin{align*} - Df(\Omega)(\theta) &=-\int_{\Gamma_S}\theta\cdot n\left[\nabla(u-U)a_{\pm}\cdot n\right]\\ - &=-\int_{\Gamma_S}\left[2\mu e(u^{*})\nabla (u-U)n\cdot n \pm p \nabla(u-U)n\cdot n\right]. -\end{align*}\] -
-
-
-

Lemma 3 allows us to vanish the last term when the first term can be rewritten using Lemma 4, which gives us

-
-
-
-\[\begin{equation*} - Df(\Omega)(\theta)=-2\mu\int_{\Gamma_S}\theta\cdot n\left[ e(u^{*}):e(u-U)\right]. -\end{equation*}\] -
-
-
-

The reasoning is exactly the same when we exchange the roles of \((u,p)\) and \((u^{*},p^{*})\), and by taking \(U=\alpha\).

-
-
-
-
-
-
Shape gradient for the rigid object in Stokes flow problem
-
-

We have the following shape gradient

-
-
-
-\[\begin{equation} - G(\Omega)=-2\mu\left[e(u^{*}):e(u)-e(u):e(\alpha)-e(u^{*}):e(U)\right]. -\end{equation}\] -
-
-
-

Of particular note, if we assume the that \(U\) and \(\alpha\) are taking in table of entries of grand resistance tensor, then

-
-
-
-\[\begin{equation} - G(\Omega)=-2\mu e(u^{*}):e(u). -\end{equation}\] -
-
-
-
-
-Proof -
-
-

We apply the Cea’s method assuming all the necessary regularities. We thus obtain

-
-
-
-\[\begin{multline} - \mathcal{L}_{\alpha}(\Omega,v,h_1,\lambda,\beta,q,h_2)=\int_{\Gamma_S}\sigma(v,q)n\cdot \alpha + \int_{\Omega}\left[-\mu \Delta v + \nabla q\right]\cdot h_1 + \int_{\Omega}(\nabla \cdot v)h_2 \\ - +\int_{\Gamma_S}\lambda \cdot (v-U) + \int_{\Gamma_B}\beta \cdot (v-U^{\infty}). -\end{multline}\] -
-
-
-

By taking as variables in the derivative of the Lagrangian the solutions of the primal problem \((u,p)\) and of the dual problem \((u^{*},p^{*})\), the terms \(-\Delta u + \nabla p\) as well as the divergence of \(u\) and \(u^{*}\) vanish. Moreover, since the boundary \(\Gamma_B\) is fixed, by deriving all the integrals \(\Gamma_B\) cancel out. Thus, we obtain the following relationship :

-
-
-
-\[\begin{equation*} - \frac{\partial \mathcal{L}_{\alpha}}{\partial \Omega}(\Omega,u,u^{*},\lambda^{*},\beta^{*},p,-p^{*})(\theta)= Df(\Omega)(\theta)+Dg(\Omega)(\theta), -\end{equation*}\] -
-
-
-

with

-
-
-
-\[\begin{equation*} -f(\Omega)=\int_{\Gamma_S}\sigma(v,q)n\cdot \alpha \quad \text{and}\quad - g(\Omega)=\int_{\Gamma_S}\lambda\cdot (v-U). -\end{equation*}\] -
-
-
-

To derive the function \(f\) and \(g\), we use Proposition Lemma 5. Thus, we have

-
-
-
-\[\begin{equation*} - \frac{\partial \mathcal{L}_{\alpha}}{\partial \Omega}(\Omega,u,u^{*},\lambda^{*},\beta^{*},p,-p^{*})(\theta)=2\mu \int_{\Gamma_S} \theta \cdot n \left[e(u):e(\alpha)-e(u^{*}):e(u-U)\right]. -\end{equation*}\] -
-
-
-

Finally, we obtain that

-
-
-
-\[\begin{equation}\label{stokes:shaped} - DJ(\Omega)(\theta)=-2\mu \int_{\Gamma_S}\theta \cdot n \left[e(u^{*}):e(u)-e(u):e(\alpha)-e(u^{*}):e(U)\right]. -\end{equation}\] -
-
-
-

Hence the results.

-
-
-
-
-
-

3.2.2. Experimental Evaluation

-
-

To align with the findings in [moreau_shapes_2022], we initially consider a sphere as the solid at the center of the fluid domain. The primal and dual problems are solved using \(\mathbb{P}_2\) continuous finite elements for velocity and \(\mathbb{P}_1\) continuous finite elements for pressure. The expansion problem solution is also expressed with a \(\mathbb{P}_1\) continuous finite elements. It is important to keep in mind that our assumptions differ from those in the paper [moreau_shapes_2022]. In particular, we assume that the \(\Gamma_B\) edge is infinitely far from the \(\Gamma_S\) edge, which is not achievable using finite elements. Consequently, the obtained results may differ. The results obtained for the GD method and the NSGF method will be presented.

-
-
-

2D simulation for \(K_{11}\) case :

-
-
-

We consider the 2D \(K_{11}\) resistance problem. The initial fluid domain is a square with a spherical hole at its center. The parameter values chosen are described in the following table and with the following initial domain.

-
-
-
-2D stokes init -
-
Figure 10. Initial geometry for the 2D case of Stokes.
-
- - ---- - - - - - - - - - - - - - - - - - - - - - - -
Table 14. Geometric and physics parameters of the initial domain for the 2D case of Stokes. The viscosity of the fluid is designed by \(\mu\). \(C\) represents the center of the sphere and the box, \(R\) the radius of the sphere and \(L\) the side length of the box.

Symbol

Value (dimensionless)

\(\mu\)

\(1\)

\(C\)

\((0,0)\)

\(R\)

\(1\)

\(L\)

\(10\)

-
-

For this specific test case, we employed an initial discretization of \(h=0.2\) around the spherical surface and \(h=1\) around the square. The optimization parameters we selected were \(l=20\), \(a=b=0.5\), and \(c=1e3\). In the NSGF method, remeshing plays a crucial role, particularly in the 2D case. To avoid convergence issues, we fixed the mesh at the edge of the solid, as the deformation was not significant. Consequently, the initial mesh was much more refined, with \(h=0.05\) at the solid boundary.

-
- - -------- - - - - - - - - - - - - - - - - - - - - - - - - - - -
Table 15. Comparison of the cost function between the initial domain and the final domain, the volume error between the two domains, the \(H^1\)-norm of the displacement field at the end and the number of iterations for the resistance problem \(K_{11}\) in 2D.

Method

\(J(\Omega_0)\)

\(J(\Omega_{n_{final}})\)

\(||\Omega_0|-|\Omega_{n_{final}}||\)

\(||\theta_{n_{final}}||_{H^1}\)

\(n_{final}\)

GD

\(16.7729\)

\(14.2940\)

\(5.2101e-05\)

\(1.0154e-4\)

\(500\)

NSGF

\(16.9087\)

\(14.5917\)

\(1.1554e-3\)

\(3.3883e-4\)

\(500\)

-
-
-results sphere2D 2D l20 t0.01 h0.2 a0.5b0.5 c1000 n500 -
-
-
-
-results nsgf t 2D h0.05 n020 n500 -
-
Figure 11. Evolution of the cost function, the volume of the domain \(\Omega_n\) and the \(H^1\)-norm of \(\theta_n\) for the 2D \(K_{11}\) resistance problem with the GD method (top) and the NSGF method (bottom).
-
-
-
-2D stokes initfinal -
-
Figure 12. Initial domain (left) and final domain (right) for the GD method.
-
-
-

The shape obtained in 2D closely resembles a rugby ball, as described in [moreau_shapes_2022]. It effectively preserves the initial volume and achieves a lower cost function compared to the initial domain. The results obtained with the two methods exhibit striking similarity.

-
-
-

3D simulations for various resistance problem :

-
-
-

The initial fluid representation is a cube with a spherical hole at its center.

-
- - ---- - - - - - - - - - - - - - - - - - - - - - - -
Table 16. Geometric and physics parameters of the initial domain for the 3D case of Stokes. The viscosity of the fluid is designed by \(\mu\). \(C\) represents the center of the sphere and the box, \(R\) the radius of the sphere and \(L\) the side length of the box.

Symbol

Value (dimensionless)

\(\mu\)

\(1\)

\(C\)

\((0,0,0)\)

\(R\)

\(1\)

\(L\)

\(10\)

-
-

The parameter values chosen are described in the table below.

-
-
-
-3D stokes init -
-
Figure 13. Illustration of the stokes problem. The \(\Gamma_B\) boundary is fixed. The \(\Gamma_S\) boundary is movable in order to optimize the shape.
-
-
-

In the following sections, we investigate the cases \(K_{11}\), \(K_{12}\), \(Q_{11}\), \(Q_{12}\), and \(C_{11}\) to compare our results with those presented in [moreau_shapes_2022]. For all simulations except the \(Q_{12}\) case, we use the optimization parameters \(t=0.03\), \(l=20\), \(a=b=0.5\), and \(c=1e3\). In the \(Q_{12}\) case, we employ \(t=5e-4\), \(l=20\), \(a=b=0.5\), and \(c=1e5\).

-
-
-

\(K_{11}\) resistance problem :

-
-
-

Comparison of the cost function between the initial domain and the final domain, the volume error between the two domains, the \(H^1\)-norm of the displacement field at the end and the number of iterations for the resistance problem \(K_{11}\) in 3D.

-
- -------- - - - - - - - - - - - - - - - - - - - - - - - - - - -

Method

\(J(\Omega_0)\)

\(J(\Omega_{n_{final}})\)

\(||\Omega_0|-|\Omega_{n_{final}}||\)

\(||\theta_{n_{final}}||_{H^1}\)

\(n_{final}\)

GD

\(29.6246\)

\(27.2228\)

\(2.9460e-4\)

\(5.0356e-4\)

\(125\)

NSGF

\(29.6246\)

\(27.6108\)

\(1.0395e-2\)

\(6.6623e-2\)

\(33\)

-
-
-results K11 3D l20 t0.03 h0.2 a0.5b0.5 c1000 n500 -
-
-
-
-results nsgf tK11 3D h0.2 n020 n502 -
-
Figure 14. Evolution of the cost function, the volume of the domain \(\Omega_n\) and the \(H^1\)-norm of \(\theta_n\) for the 3D \(K_{11}\) resistance problem with the GD method (top) and the NSGF method (bottom).
-
-
-
\(K_{11}\) simulation for the GD method.
-
- -
-
-
-

We successfully achieve the optimal rugby ball shape while preserving volume and reducing the cost function. The NSGF method seems to be highly sensitive to remeshing, particularly at the ends, resulting in a shape reminiscent of a lemon. As a result, the simulation needs to be terminated before the mesh collapses at the ends.

-
-
-

\(K_{12}\) resistance problem :

-
- - -------- - - - - - - - - - - - - - - - - - - - - - - - - - - -
Table 17. Comparison of the cost function between the initial domain and the final domain, the volume error between the two domains, the \(H^1\)-norm of the displacement field at the end and the number of iterations for the resistance problem \(K_{12}\) in 3D.

Method

\(J(\Omega_0)\)

\(J(\Omega_{n_{final}})\)

\(||\Omega_0|-|\Omega_{n_{final}}||\)

\(||\theta_{n_{final}}||_{H^1}\)

\(n_{final}\)

GD

\(-0.0188\)

\(-18.3755\)

\(0.2041\)

\(9.1405e-2\)

\(45\)

NSGF

\(-0.0188\)

\(-17.9411\)

\(8.5382e-3\)

\(2.5186e-3\)

\(850\)

-
-
-results K12 3D l20 t0.03 h0.2 a0.5b0.5 c1000 n500 -
-
-
-
-results nsgf tK12 Bfixed 3D h0.2 n020 n1000 -
-
Figure 15. Evolution of the cost function, the volume of the domain \(\Omega_n\) and the \(H^1\)-norm of \(\theta_n\) for the 3D \(K_{12}\) resistance problem with the GD method (top) and the NSGF method (bottom).
-
-
-
\(K_{12}\) simulation for the GD method.
-
- -
-
-
-

As depicted in the previous figures, the simulations are terminated before the cost function reaches a local minimum due to mesh collapse. There are minimal differences in the final shape obtained between the two methods. Nevertheless, we manage to achieve a shape that closely resembles the one presented in the reference article, although it should be noted that the problem we are studying differs slightly as the boundary is at a finite distance from the rigid object.

-
-
-

\(Q_{11}\) resistance problem :

-
- - -------- - - - - - - - - - - - - - - - - - - - - - - - - - - -
Table 18. Comparison of the cost function between the initial domain and the final domain, the volume error between the two domains, the \(H^1\)-norm of the displacement field at the end and the number of iterations for the resistance problem \(Q_{11}\) in 3D.

Method

\(J(\Omega_0)\)

\(J(\Omega_{n_{final}})\)

\(||\Omega_0|-|\Omega_{n_{final}}||\)

\(||\theta_{n_{final}}||_{H^1}\)

\(n_{final}\)

GD

\(24.1304\)

\(16.9092\)

\(5.4293e-3\)

\(5.0132e-3\)

\(500\)

NSGF

\(24.1304\)

\(16.9188\)

\(5.7716e-3\)

\(1.7381e-3\)

\(873\)

-
-
-results Q11 3D l20 t0.03 h0.2 a0.5b0.5 c1000 n500 -
-
-
-
-results nsgf tQ11 Bfixed 3D h0.2 n020 n1000 -
-
Figure 16. Evolution of the cost function, the volume of the domain \(\Omega_n\) and the \(H^1\)-norm of \(\theta_n\) for the 3D \(Q_{11}\) resistance problem with the GD method (top) and the NSGF method (bottom).
-
-
-
\(Q_{11}\) simulation for the GD method.
-
- -
-
-
-

In this case, premature termination of the simulation is less evident, and convergence is observed with good volume retention. The shapes obtained are virtually identical between the two resolution methods. However, a notable observation is that the ends of the object approach the edges of the cube until they collide with them. This phenomenon is not observed when the cube boundaries are assumed to be at infinity, as depicted in [moreau_shapes_2022].

-
-
-

\(Q_{12}\) resistance problem :

-
- - -------- - - - - - - - - - - - - - - - - - - - - - - - - - - -
Table 19. Comparison of the cost function between the initial domain and the final domain, the volume error between the two domains, the \(H^1\)-norm of the displacement field at the end and the number of iterations for the resistance problem \(Q_{12}\) in 3D.

Method

\(J(\Omega_0)\)

\(J(\Omega_{n_{final}})\)

\(||\Omega_0|-|\Omega_{n_{final}}||\)

\(||\theta_{n_{final}}||_{H^1}\)

\(n_{final}\)

GD

\(-9.0952e-3\)

\(-323.3069\)

\(3.9189e-2\)

\(3.7323e-2\)

\(980\)

NSGF

\(-9.0952e-3\)

\(-234.4491\)

\(1.5241e-2\)

\(6.3734e-3\)

\(1000\)

-
-
-results Q12 3D l20 t0.0005 h0.2 a0.5b0.5 c100000 n1000 -
-
-
-
-results nsgf tQ12 Bfixed 3D h0.2 n020 n1000 -
-
Figure 17. Evolution of the cost function, the volume of the domain \(\Omega_n\) and the \(H^1\)-norm of \(\theta_n\) for the 3D \(Q_{12}\) resistance problem with the GD method (top) and the NSGF method (bottom).
-
-
-
\(Q_{12}\) simulation for the GD method.
-
- -
-
-
-

The boundary of the cube plays a significant role in the obtained results. Contrary to the expected dumbbell shape, we observe the ends of the solid flattening, leading to mesh collapse and an abrupt termination of the simulation. Additionally, these ends are positioned too closely to the cube boundary, resembling the observations made in the \(Q_{11}\) case. The shapes obtained using the two different methods exhibit a notable similarity.

-
-
-

\(C_{11}\) resistance problem :

-
- - -------- - - - - - - - - - - - - - - - - - - - - - - - - - - -
Table 20. Comparison of the cost function between the initial domain and the final domain, the volume error between the two domains, the \(H^1\)-norm of the displacement field at the end and the number of iterations for the resistance problem \(C_{11}\) in 3D.

Method

\(J(\Omega_0)\)

\(J(\Omega_{n_{final}})\)

\(||\Omega_0|-|\Omega_{n_{final}}||\)

\(||\theta_{n_{final}}||_{H^1}\)

\(n_{final}\)

GD

\(2.7004e-3\)

\(-1.2593\)

\(4.2977e-2\)

\(4.4303e-2\)

\(74\)

NSGF

\(2.7004e-3\)

\(-1.7085\)

\(9.7921e-4\)

\(1.2600e-3\)

\(450\)

-
-
-results C11 3D l20 t0.03 h0.2 a0.5b0.5 c1000 n500 -
-
-
-
-results nsgf tC11 Bfixed 3D h0.2 n020 n1000 -
-
Figure 18. Evolution of the cost function, the volume of the domain \(\Omega_n\) and the \(H^1\)-norm of \(\theta_n\) for the 3D \(C_{11}\) resistance problem with the GD method (top) and the NSGF method (bottom).
-
-
-
\(C_{11}\) simulation for the GD method.
-
- -
-
-
-

Mesh collapse continues to be evident in this case, mainly due to sharp edges forming in specific regions. However, the dynamics align with the description provided in [moreau_shapes_2022], and volume conservation is satisfactory. Interestingly, we observe the emergence of four blades on the solid, each with varying prominence. It is worth noting that the difference between the two methods is apparent in figure below, where the blades are not necessarily in the same positions.

-
-
-
-C11 blades -
-
Figure 19. The black wireform corresponds to the result obtained using the GD method. The grey surface is the shape obtained using the NSGF method.
-
-
-
-
-
-
-

4. Implementation in Feel++

-
- -
-
-
-

References on geometric shape optimisation

-
-
-
    -
  • -

    [internship_Palazzolo] Lucas Palazzolo. Shape optimization for rigid objects in Stokes flow. Stage de M2, Univerité de Strasbourg/Inria Sophia-Antipolis, 2023. repo

    -
    • -
  • -

    [allaire_conception_2006] Grégoire Allaire. Conception optimale de structures. volume 58 of Mathématiques & applications. Springer Berlin Heidelberg, 2006

    -
    • -
  • -

    [cea_conception_1986] Jean Céa. Conception optimale ou identification de formes, calcul rapide de la dérivée directionnelle de la fonction coût. M2AN - Modélisation mathématique et analyse numérique, 20(3):371-402, 1986

    -
    • -
  • -

    [moreau_shapes_2022] Clément Moreau, Kenta Ishimoto and Yannick Privat. Shapes optimising grand resistance tensor entries for a rigid body in a Stokes flow. July 2022

    -
    • -
  • -

    [feppon_shape_2019] Florian Feppon. Shape and topology optimization of multiphysics systems. These de doctorat, Université Paris-Saclay (ComUE), Dec. 2019

    -
    • -
-
-
-
-
-
-
-
- - - - - - - - - - - - diff --git a/shapo/githubactions.html b/shapo/githubactions.html deleted file mode 100644 index 3eb5b57..0000000 --- a/shapo/githubactions.html +++ /dev/null @@ -1,365 +0,0 @@ - - - - - - Github Actions :: Feel++ ShapO // Docs - - - - - - - - - - - - - - - - - - - - - - -
- -
-
- -
- -
-
- - -
- -
-

Github Actions

-
-

Github actions are setup to

-
-
-
    -
  • -

    build the documentation using Antora and upload them to docs.feelpp.org[Feel++ docs] upon any change in the repo.

    -
    • -
  • -

    build a sample Feel++ code

    -
    • -
-
-
-
-
-
- - - - - - - - - - - - diff --git a/shapo/rename.html b/shapo/implementation.html similarity index 90% rename from shapo/rename.html rename to shapo/implementation.html index 6fd8fdf..d68e71b 100644 --- a/shapo/rename.html +++ b/shapo/implementation.html @@ -3,8 +3,8 @@ - Renaming the project :: Feel++ ShapO // Docs - + Implementation in Feel++ :: Feel++ ShapO // Docs + @@ -172,34 +172,10 @@

Feel++ ShapO

    • @@ -227,13 +203,12 @@

    Feel++ ShapO

    -
    Edit this Page
    +
    Edit this Page
    Download as @@ -253,12 +228,9 @@

    Feel++ ShapO

    -

    Renaming the project

    -
    -

    The script rename.sh renames the project.

    -
    +

    Implementation in Feel++

    -

    You have to answer a few questions to setup a new project out of this template.

    +

    The code is based on Lucas Palazzolo’s course. For more details, see [internship_Palazzolo].

    diff --git a/shapo/index.html b/shapo/index.html index 4a44ebf..8a6869d 100644 --- a/shapo/index.html +++ b/shapo/index.html @@ -1,10 +1,11 @@ - - + + - Feel++ ShapO Project :: Feel++ ShapO // Docs + Shape Optimisation :: Feel++ ShapO // Docs + @@ -92,8 +93,8 @@ --> - - + +
    +
    - +
    -
    - + +
    + + +