This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License
本作品采用 Creative Commons Attribution-ShareAlike 4.0 International License 进行许可

The corresponding files can be obtained from:
相应的文件可以从以下位置获得：

Jupyter Notebook: simp_topology_optimization.ipynb
Jupyter 笔记本： simp_topology_optimization.ipynb
Python script: simp_topology_optimization.py
Python 脚本：simp_topology_optimization.py

Topology optimization using the SIMP method¶
使用 SIMP 方法进行拓扑优化 ¶

This numerical tour investigates the optimal design of an elastic structure through a topology optimization approach. It consists in finding the optimal distribution of a material in a computational domain which minimizes the compliance (or, equivalently, maximizes the stiffness) of the resulting structure under a fixed volume fraction constraint.
此数值导览通过拓扑优化方法研究弹性结构的优化设计。它包括在计算域中找到材料的最佳分布，该分布使在固定体积分数约束下最终结构的柔度最小（或等效地，使刚度最大化）。

For instance, in the bottom animation, the optimal material density (a continuous field $θ (x) \in [0; 1]$ ) of a cantilever beam evolves with the topoology optimization algorithm iterations. The first 20 iterations correspond to a “non-penalized” distribution where intermediate densities (i.e. $0 < θ (x) < 1$ ) are allowed. However, such a gray-level distribution is hard to realize in practice and one would like to obtain a black-and-white design of material/void distribution. This is what is achieved in the last iterations where some penalty is applied to intermediate densities, exhibiting finally a truss-like design. Note that the final design is not a global optimum but only a local one. Hence, different designs can well be obtained if changing the initial conditions or the mesh density for instance.
例如，在底部动画中，悬臂梁的最佳材料密度（连续场 $θ (x) \in [0; 1]$ ）随着拓扑优化算法的迭代而变化。前 20 次迭代对应于允许中等密度（即 $0 < θ (x) < 1$ ）的 “非惩罚” 分布。然而，这样的灰度分布在实践中很难实现，人们希望获得材料/空隙分布的黑白设计。这就是在最后一次迭代中实现的，其中对中等密度施加了一些惩罚，最终展示了类似桁架的设计。请注意，最终设计不是全局最优设计，而只是局部设计。因此，例如，如果更改初始条件或网格密度，则可以很好地获得不同的设计。

The present tour relies on one of the most widely implemented approach in topology optimization, namely the Solid Isotropic Material with Penalisation (SIMP) method developed by Bendsoe and Sigmund [BEN03]. See also the monograph [ALL07] and Gregoire Allaire’s class on *Optimal Design of Structures* at Ecole Polytechnique. In particular, this tour refers to lesson 8 and is the direct adaptation of the SIMP Topology Optimization example of the CMAP group toolbox written in Freefem++.
本导览依赖于拓扑优化中应用最广泛的方法之一，即由 Bendsoe 和 Sigmund [BEN03] 开发的具有惩罚作用的固体各向同性材料（SIMP）方法。另请参见专著 [ALL07] 和 Gregoire Allaire 在巴黎综合理工学院的 *最佳结构设计* 课程。具体而言，本教程参考了第 8 课，并且直接改编了用 Freefem++ 编写的 CMAP 组工具箱的 SIMP 拓扑优化示例。

This tour should also work in parallel.
这次旅行也应该同时进行。

Theory¶ 理论¶

Compliance optimization on binary shapes¶
二进制 shape 的合规性优化¶

Let us consider a computational domain $D$ in which we seek an elastic material domain $Ω$ ( $D ∖ Ω$ representing void) so that the compliance under prescribed loading is minimal when subjected to a fixed volume constraint $Vol (Ω) = η Vol (D)$ . Replacing the shape $Ω$ by its characteristic function $χ (x) \in {0, 1}$ , the previous problem can be formulated as the following optimization problem:
让我们考虑一个计算域，在该域中 $D$ ，我们寻求一个弹性材料域 $Ω$ （ $D ∖ Ω$ 表示空隙），以便在受到固定体积约束时，在规定载荷下的柔度最小。 $Vol (Ω) = η Vol (D)$ 用 shape 的特征函数 $χ (x) \in {0, 1}$ 替换 shape $Ω$ ，前面的问题可以表述为以下优化问题：

\begin{aligned} \begin{aligned} min_{χ} & \int_{D} f \cdot u (χ) dx + \int_{Γ_{N}} T \cdot u (χ) dS \\ s.t. & \int_{D} χ (x) dx = η Vol (D) \\ χ (x) \in {0, 1} \end{aligned} \end{aligned}

where $u (χ)$ is the solution of the following $χ (x) -$ dependent elasticity problem:
其中 $u (χ)$ 是以下 $χ (x) -$ 依赖弹性问题的解：

\begin{aligned} \begin{aligned} div σ + f = 0 & in D \\ σ = χ (x) C : \nabla^{s} u & in D \\ σ \cdot n = T & on Γ_{N} \\ u = 0 & on Γ_{D} \end{aligned} \end{aligned}

where $f$ is a body force, $T$ corresponds to surface tractions on $Γ_{N}$ and $C$ is the elastic material stiffness tensor. Owing to the elastic nature of the problem, the compliance can also be expressed using the elastic stress energy density:
其中 $f$ 是体力， $T$ 对应于 ON $Γ_{N}$ 的表面牵引力， $C$ 是弹性材料刚度张量。由于问题的弹性性质，柔度也可以用弹性应力能量密度来表示：

\int_{D} f \cdot u (χ) dx + \int_{Γ_{N}} T \cdot u (χ) dS = \int_{D} σ (χ) : (χ (x) C)^{- 1} : σ (χ) dx = \int_{D} \nabla^{s} u (χ) : (χ (x) C) : \nabla^{s} u (χ) dx

so that the above optimization problem can be reformulated as:
这样上面的优化问题就可以改述为：

\begin{aligned} \begin{aligned} min_{χ, u} & \int_{D} σ (χ) : (χ (x) C)^{- 1} : σ (χ) dx \\ s.t. & div σ + f = 0 in D \\ σ \cdot n = T on Γ_{N} \\ \int_{D} χ (x) dx = η Vol (D) \\ χ (x) \in {0, 1} \end{aligned} \end{aligned}

Continuous relaxation and SIMP penalization¶
Continuous relaxing 和 SIMP 惩罚¶

The binary constraint $χ (x) \in {0, 1}$ makes the above problem extremely difficult to solve so that a classical remedy consists in relaxing the constraint by a continuous constraint $θ (x) \in [0, 1]$ where $θ (x)$ will approximate the characteristic function $χ (x)$ by a gray-level continuous function taking values between 0 and 1.
二进制约束 $χ (x) \in {0, 1}$ 使上述问题极难解决，因此经典的补救措施包括通过连续约束来放松约束，其中将通过取值在 0 和 1 之间的灰度级连续函数来近似特征函数。 $θ (x) \in [0, 1]$ $χ (x)$ $θ (x)$

However, in order to still obtain a final binary density distribution, the binary modulus $χ (x) C$ will be replaced by $θ (x)^{p} C$ with $p > 1$ in order to penalize intermediate densities, yielding the so-called SIMP formulation:
然而，为了仍然获得最终的二进制密度分布，二进制模 $χ (x) C$ 量将被替换为 $θ (x)^{p} C$ with $p > 1$ 以惩罚中间密度，从而产生所谓的 SIMP 公式 ：

\begin{aligned} \begin{aligned} min_{θ, σ} J_{p} (θ, σ) = min_{θ, σ} & \int_{D} σ : (θ (x)^{- p} C^{- 1}) : σ dx \\ s.t. & div σ + f = 0 in D \\ σ \cdot n = T on Γ_{N} \\ \int_{D} θ (x) dx = η Vol (D) \\ 0 \leq θ (x) \leq 1 \end{aligned} \end{aligned}

Optimization of the SIMP formulation¶
SIMP 公式的优化¶

Unfortunately jointly minimizing $J_{p}$ over $(θ, σ)$ is hard to do in practice since this functional is non-convex (except if $p = 1$ ). However, it is convex over each variable $θ$ and $σ$ when fixing the other variable. This makes alternate minimization an attractive method for finding a local optimum. Besides, minimizing directly for a fixed value of $p > 1$ does not work well in practice. A better solution consists in performing alternate minimization steps and progressively increase $p$ from 1 to a maximum value (typically 3 or 4) using some heuristic (see later).
不幸的是，在实践中很难共同最小化 $J_{p}$ over $(θ, σ)$ ，因为这个函数是非凸的（除非 $p = 1$ ）。但是，它在每个变量 $θ$ 上和 $σ$ 固定另一个变量时都是凸的。这使得交替最小化成为查找局部最优值的有吸引力的方法。此外，直接对固定值 of 进行最小化在实践中效果不佳。 $p > 1$ 更好的解决方案包括执行交替的最小化步骤，并使用一些启发式方法（见下文）逐渐从 1 增加到 $p$ 最大值（通常为 3 或 4）。

Iteration $n + 1$ of the algorithm, knowing a previous pair $(θ_{n}, σ_{n})$ and a given value of the penalty exponent $p_{n}$ , therefore consists in:
算法的迭代 $n + 1$ ，知道前一对 $(θ_{n}, σ_{n})$ 和给定的惩罚指数 $p_{n}$ 值，因此包括：

minimizing $J_{p_{n}} (θ_{n}, σ)$ yielding $σ_{n + 1} = (θ_{n})^{p_{n}} C : \nabla^{s} u_{n + 1}$
最小化 $J_{p_{n}} (θ_{n}, σ)$ 产量 $σ_{n + 1} = (θ_{n})^{p_{n}} C : \nabla^{s} u_{n + 1}$
minimizing $J_{p_{n}} (θ, σ_{n + 1})$ yielding $θ_{n + 1}$
最小化 $J_{p_{n}} (θ, σ_{n + 1})$ 产量 $θ_{n + 1}$
potentially update the value of the exponent using some heuristic $p_{n} \to p_{n + 1}$
可能使用一些启发式 $p_{n} \to p_{n + 1}$ 方法更新指数的值

Both minimization problems are quite easy to solve since the first one is nothing else than a heterogeneous elastic problem with some known modulus $(θ_{n} (x))^{p_{n}} C$ .
这两个最小化问题都很容易解决，因为第一个问题只不过是具有一些已知模量的异质弹性问题。 $(θ_{n} (x))^{p_{n}} C$

The second problem can be rewritten using the Lagrange multiplier $λ$ corresponding to the toal volume constraint, as follows:
第二个问题可以使用对应于 toal 体积约束的拉格朗日乘子 $λ$ 重写，如下所示：

min_{θ \in [0; 1]} \int_{D} θ (x)^{- p_{n}} e (σ_{n + 1}) dx + λ (\int_{D} θ (x) dx - η Vol (D))

where $e (σ_{n + 1}) = σ_{n + 1} : C^{- 1} : σ_{n + 1} = (θ_{n})^{p_{n}} σ_{n + 1} : \nabla^{s} u_{n + 1}$ . The optimality conditions for this problem yields the following explicit and local condition:
其中 $e (σ_{n + 1}) = σ_{n + 1} : C^{- 1} : σ_{n + 1} = (θ_{n})^{p_{n}} σ_{n + 1} : \nabla^{s} u_{n + 1}$ .此问题的最优性条件产生以下显式局部条件：

- p_{n} θ_{n + 1}^{- p_{n} - 1} e (σ_{n + 1}) + λ = 0

which along with the $[θ_{m i n}; 1]$ constraint gives:
它与约束一起 $[θ_{m i n}; 1]$ 得到：

θ_{n + 1} = min {1; max {θ_{m i n}; {(\frac{p_{n} e (σ_{n + 1})}{λ})}^{1 / (p_{n} + 1)}}}

where we replaced the $0$ constraint by a minimum density value $θ_{m i n} > 0$ to avoid degeneracy issue with void material.
其中， $0$ 我们将约束替换为 Minimum Density 值 $θ_{m i n} > 0$ ，以避免 void 材质的简并问题。

Note that the above expression assumes that $λ$ is known. Its value is found by satsifying the volume constraint $\int_{D} θ (x) dx = η Vol (D)$ .
请注意，上面的表达式假定 $λ$ this is known 。它的值是通过满足 volume constraint $\int_{D} θ (x) dx = η Vol (D)$ 来找到的。

FEniCS implementation¶ FEniCS 实现¶

We first define some parameters of the algorithm. The most important ones concern the number of total alternate minimizations (AM) niter and the parameters controlling the exponent update strategy:
我们首先定义算法的一些参数。最重要的参数涉及总交替最小化（AM） 的数量和控制指数更新策略的参数：

niternp corresponds to the number of first (AM) for which $p = 1$ . These are non-penalized iterations yielding a diffuse gray-level field $θ (x)$ at convergence (note that we solve this convex problem with AM iterations although one could use a dedicated convex optimization solver).
niternp 对应于 . $p = 1$ 这些是未受惩罚的迭代，在收敛时产生漫反射灰度场 $θ (x)$ （请注意，我们用 AM 迭代解决了这个凸问题，尽管可以使用专用的凸优化求解器）。
pmax is the maximum value taken by the penalty exponent $p$
pmax 是罚指数 $p$ 取的最大值
exponent_update_frequency corresponds to the minimum number of AM iterations between two consecutive updates of the exponent
exponent_update_frequency 对应于指数的两次连续更新之间的最小 AM 迭代次数

[8]:

%matplotlib notebook
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np

# Algorithmic parameters
niternp = 20 # number of non-penalized iterations
niter = 80 # total number of iterations
pmax = 4        # maximum SIMP exponent
exponent_update_frequency = 4 # minimum number of steps between exponent update
tol_mass = 1e-4 # tolerance on mass when finding Lagrange multiplier
thetamin = 0.001 # minimum density modeling void

We now define the problem parameters, in particular the target material density $(η = 40 %$ here). The mesh consists of a rectangle of dimension $4 \times 1$ , clamped on its left side and loaded by a uniformly distributed vertical force on a line of length 0.1 centered around the center of the right side.
现在，我们在此处定义问题参数，特别是目标材料密度 $(η = 40 %$ ）。网格由一个维度 $4 \times 1$ 的矩形组成，该矩形夹在其左侧，并由均匀分布的垂直力加载在以右侧中心为中心的长度为 0.1 的线上。

Finally, we initialized the SIMP penalty exponent to $p = 1$ and initialized also the density field and the Lagrange multiplier $λ$ .
最后，我们将 SIMP 惩罚指数初始化为密度场和拉格朗日乘子，并初始化了密度场和拉格朗日乘子 $λ$ 。 $p = 1$

[9]:

# Problem parameters
thetamoy = 0.4 # target average material density
E = Constant(1)
nu = Constant(0.3)
lamda = E*nu/(1+nu)/(1-2*nu)
mu = E/(2*(1+nu))
f = Constant((0, -1)) # vertical downwards force

# Mesh
mesh = RectangleMesh(Point(-2, 0), Point(2, 1), 50, 30, "crossed")
# Boundaries
def left(x, on_boundary):
    return near(x[0], -2) and on_boundary
def load(x, on_boundary):
    return near(x[0], 2) and near(x[1], 0.5, 0.05)
facets = MeshFunction("size_t", mesh, 1)
AutoSubDomain(load).mark(facets, 1)
ds = Measure("ds", subdomain_data=facets)

# Function space for density field
V0 = FunctionSpace(mesh, "DG", 0)
# Function space for displacement
V2 = VectorFunctionSpace(mesh, "CG", 2)
# Fixed boundary condtions
bc = DirichletBC(V2, Constant((0, 0)), left)

p = Constant(1) # SIMP penalty exponent
exponent_counter = 0 # exponent update counter
lagrange = Constant(1) # Lagrange multiplier for volume constraint

thetaold = Function(V0, name="Density")
thetaold.interpolate(Constant(thetamoy))
coeff = thetaold**p
theta = Function(V0)

volume = assemble(Constant(1.)*dx(domain=mesh))
avg_density_0 = assemble(thetaold*dx)/volume # initial average density
avg_density = 0.

We now define some useful functions for formulating the linear elastic variational problem.
现在，我们定义了一些有用的函数来公式化线性弹性变分问题。

[10]:

def eps(v):
    return sym(grad(v))
def sigma(v):
    return coeff*(lamda*div(v)*Identity(2)+2*mu*eps(v))
def energy_density(u, v):
    return inner(sigma(u), eps(v))

# Inhomogeneous elastic variational problem
u_ = TestFunction(V2)
du = TrialFunction(V2)
a = inner(sigma(u_), eps(du))*dx
L = dot(f, u_)*ds(1)

We now define the function for updating the value of $θ$ . Note that this function will be called many times in each step since we are finding $λ$ through a dichotomy procedure. For this reason, we make use of the local_project function (see Elasto-plastic analysis of a 2D von Mises material and Efficient projection on DG or Quadrature spaces) for fast projection on local spaces such as DG0 for the present case.
现在，我们定义用于更新的值 $θ$ 的函数。请注意，此函数在每个步骤中将被多次调用，因为我们是通过二分法过程查找 $λ$ 的。出于这个原因，我们利用 local_project 函数（参见 2D von Mises 材料的弹塑性分析和 DG 或正交空间上的高效投影）在局部空间上进行快速投影，例如本例中的 DG0。

[11]:

def local_project(v, V):
    dv = TrialFunction(V)
    v_ = TestFunction(V)
    a_proj = inner(dv, v_)*dx
    b_proj = inner(v, v_)*dx
    solver = LocalSolver(a_proj, b_proj)
    solver.factorize()
    u = Function(V)
    solver.solve_local_rhs(u)
    return u
def update_theta():
    theta.assign(local_project((p*coeff*energy_density(u, u)/lagrange)**(1/(p+1)), V0))
    thetav = theta.vector().get_local()
    theta.vector().set_local(np.maximum(np.minimum(1, thetav), thetamin))
    theta.vector().apply("insert")
    avg_density = assemble(theta*dx)/volume
    return avg_density

We now define a function for finding the correct value of the Lagrange multiplier $λ$ . First, a rough bracketing of $λ$ is obtained, then a dichotomy is performed in the interval [lagmin,lagmax] until the correct average density is obtained to a certain tolerance.
现在，我们定义一个函数来查找拉格朗日乘数 $λ$ 的正确值。首先，获得粗略的括号 $λ$ ，然后在区间 [lagmin，lagmax] 中执行二分法，直到获得正确的平均密度到一定的容差。

[12]:

def update_lagrange_multiplier(avg_density):
    avg_density1 = avg_density
    # Initial bracketing of Lagrange multiplier
    if (avg_density1 < avg_density_0):
        lagmin = float(lagrange)
        while (avg_density < avg_density_0):
            lagrange.assign(Constant(lagrange/2))
            avg_density = update_theta()
        lagmax = float(lagrange)
    elif (avg_density1 > avg_density_0):
        lagmax = float(lagrange)
        while (avg_density > avg_density_0):
            lagrange.assign(Constant(lagrange*2))
            avg_density = update_theta()
        lagmin = float(lagrange)
    else:
        lagmin = float(lagrange)
        lagmax = float(lagrange)

    # Dichotomy on Lagrange multiplier
    inddico=0
    while ((abs(1.-avg_density/avg_density_0)) > tol_mass):
        lagrange.assign(Constant((lagmax+lagmin)/2))
        avg_density = update_theta()
        inddico += 1;
        if (avg_density < avg_density_0):
            lagmin = float(lagrange)
        else:
            lagmax = float(lagrange)
    print("   Dichotomy iterations:", inddico)

Finally, the exponent update strategy is implemented:
最后，实现指数更新策略：

first, $p = 1$ for the first niternp iterations
首先， $p = 1$ 对于第一个 Niternp 迭代
then, $p$ is increased by some amount which depends on the average gray level of the density field computed as $g = \frac{1}{Vol(D)} \int_{D} 4 (θ - θ_{m i n}) (1 - θ) dx$ , that is $g = 0$ is $θ (x) = θ_{m i n}$ or $1$ everywhere and $g = 1$ is $θ = (θ_{m i n} + 1) / 2$ everywhere.
然后， $p$ 增加一定量，这取决于密度场的平均灰度级别，计算公式为 $g = \frac{1}{Vol(D)} \int_{D} 4 (θ - θ_{m i n}) (1 - θ) dx$ ，即 $g = 0$ is $θ (x) = θ_{m i n}$ 或 $1$ everywhere 和 $g = 1$ is $θ = (θ_{m i n} + 1) / 2$ everywhere。
Note that $p$ can increase only if at least exponent_update_frequency AM iterations have been performed since the last update and only if the compliance evolution falls below a certain threshold.
请注意，仅当自上次更新以来至少执行了 exponent_update_frequency 次 AM 迭代，并且合规性演变低于特定阈值时，才能增加。 $p$

[13]:

def update_exponent(exponent_counter):
    exponent_counter += 1
    if (i < niternp):
        p.assign(Constant(1))
    elif (i >= niternp):
        if i == niternp:
            print("\n Starting penalized iterations\n")
        if ((abs(compliance-old_compliance) < 0.01*compliance_history[0]) and
            (exponent_counter > exponent_update_frequency) ):
            # average gray level
            gray_level = assemble((theta-thetamin)*(1.-theta)*dx)*4/volume
            p.assign(Constant(min(float(p)*(1+0.3**(1.+gray_level/2)), pmax)))
            exponent_counter = 0
            print("   Updated SIMP exponent p = ", float(p))
    return exponent_counter

Finally, the global loop for the algorithm is implemented consisting, at each iteration, of the elasticity problem resolution, the corresponding compliance computation, the update for $θ$ and its associated Lagrange multiplier $λ$ and the exponent update procedure.
最后，实现算法的全局循环，包括在每次迭代中，弹性问题解决、相应的柔度计算、更新 $θ$ 及其相关的拉格朗日乘子 $λ$ 和指数更新过程。

[ ]:

u = Function(V2, name="Displacement")
old_compliance = 1e30
ffile = XDMFFile("topology_optimization.xdmf")
ffile.parameters["flush_output"]=True
ffile.parameters["functions_share_mesh"]=True
compliance_history = []
for i in range(niter):
    solve(a == L, u, bc, solver_parameters={"linear_solver": "cg", "preconditioner": "hypre_amg"})
    ffile.write(thetaold, i)
    ffile.write(u, i)

    compliance = assemble(action(L, u))
    compliance_history.append(compliance)
    print("Iteration {}: compliance =".format(i), compliance)

    avg_density = update_theta()

    update_lagrange_multiplier(avg_density)

    exponent_counter = update_exponent(exponent_counter)

    # Update theta field and compliance
    thetaold.assign(theta)
    old_compliance = compliance

The final density is represented as well as the convergence history of the compliance. One can note that the final compliance obtained after the first non-penalized iterations is smaller than the final one. This initially obtained topology is therefore more optimal than the final one, although it is made of large diffuse gray regions (see XMDF outputs or the animation at the beginning of the tour) contrary to the final one which is close to being binary.
表示最终密度以及柔度的收敛历史。可以注意到，在第一次未受惩罚的迭代后获得的最终合规性小于最终的合规性。因此，这种最初获得的拓扑比最终的拓扑更优化，尽管它由大型漫反射灰色区域组成（参见 XMDF 输出或游览开始时的动画），而最终的拓扑则接近二进制。

[15]:

plot(theta, cmap="bone_r")
plt.title("Final density")
plt.show();

plt.figure()
plt.plot(np.arange(1, niter+1), compliance_history)
ax = plt.gca()
ymax = ax.get_ylim()[1]
plt.plot([niternp, niternp], [0, ymax], "--k")
plt.annotate(r"$\leftarrow$ Penalized iterations $\rightarrow$", xy=[niternp+1, ymax*0.02], fontsize=14)
plt.xlabel("Number of iterations")
plt.ylabel("Compliance")
plt.title("Convergence history", fontsize=16)
plt.show();

References¶ 参考资料¶

[ ]:

Topology optimization using the SIMP method¶
使用 SIMP 方法进行拓扑优化 ¶

Theory¶ 理论¶

Compliance optimization on binary shapes¶
二进制 shape 的合规性优化¶

Continuous relaxation and SIMP penalization¶
Continuous relaxing 和 SIMP 惩罚¶

Optimization of the SIMP formulation¶
SIMP 公式的优化¶

FEniCS implementation¶ FEniCS 实现¶

References¶ 参考资料¶

Table Of Contents 目录

Search 搜索

Topology optimization using the SIMP method¶使用 SIMP 方法进行拓扑优化 ¶

Theory¶ 理论¶

Compliance optimization on binary shapes¶二进制 shape 的合规性优化¶

Continuous relaxation and SIMP penalization¶Continuous relaxing 和 SIMP 惩罚¶

Optimization of the SIMP formulation¶SIMP 公式的优化¶

FEniCS implementation¶ FEniCS 实现¶

References¶ 参考资料¶

Table Of Contents 目录

Search 搜索

Topology optimization using the SIMP method¶
使用 SIMP 方法进行拓扑优化 ¶

Compliance optimization on binary shapes¶
二进制 shape 的合规性优化¶

Continuous relaxation and SIMP penalization¶
Continuous relaxing 和 SIMP 惩罚¶

Optimization of the SIMP formulation¶
SIMP 公式的优化¶