first commit
|
@ -0,0 +1,21 @@
|
|||
MIT License
|
||||
|
||||
Copyright (c) 2022 LIYu
|
||||
|
||||
Permission is hereby granted, free of charge, to any person obtaining a copy
|
||||
of this software and associated documentation files (the "Software"), to deal
|
||||
in the Software without restriction, including without limitation the rights
|
||||
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
||||
copies of the Software, and to permit persons to whom the Software is
|
||||
furnished to do so, subject to the following conditions:
|
||||
|
||||
The above copyright notice and this permission notice shall be included in all
|
||||
copies or substantial portions of the Software.
|
||||
|
||||
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
||||
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
||||
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
||||
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
||||
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
||||
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
|
||||
SOFTWARE.
|
|
@ -0,0 +1,122 @@
|
|||
# Topology-Optimization-in-Julia
|
||||
|
||||
Julia Codes for Structural Topology Optimization Design
|
||||
|
||||
## Codes
|
||||
|
||||
A personal [Julia](https://epubs.siam.org/doi/10.1137/141000671) code is given mainly based on a compact and efficient Matlab implementation **top99neo** of compliance topology optimization (TO) for 2D continua[^1], which is a v3.0 version of the celebrated **top99** Matlab code developed by Sigmund[^2] and **top88** by its heir[^3].
|
||||
|
||||
Assemble just one half of the sysmetric stiffness matrix, thus substantial speedups are acheived.
|
||||
|
||||
`top_oc/` and `top_mma/` contain corresponding files related to the TO with OC and MMA algorithms, respectively.
|
||||
|
||||
Running codes could be tried out as:
|
||||
|
||||
```julia
|
||||
include("./top99neo_mma.jl")
|
||||
setup = SetUp() # problem setup
|
||||
mat = Mat() # material property
|
||||
disfeature = DiscretizationFeature(setup, mat) # model discretization
|
||||
load = LoadsSupportsBCs(setup, disfeature) # boudary conditions
|
||||
ini = Initialization(setup, disfeature, mat) # initial conditions
|
||||
filter = Filter(setup) # filtering
|
||||
xPhys, opt_hist, vf_hist, anim = Optimization(setup, mat, load, filter, ini, disfeature) # optimization process
|
||||
gif(anim, "./res/des_hist.gif", fps=20) # design result visulization
|
||||
```
|
||||
|
||||
## Results
|
||||
|
||||
A benchmark MBB example is presented. TO design results are saved in `./res/` folder and a evolution history is shown as below.
|
||||
|
||||
<!-- ![TO design evolution with MMA](./top_mma/res/des_hist.gif) -->
|
||||
|
||||
<div align=center>
|
||||
<img src=./top_mma/res/des_hist.gif width="500">
|
||||
👍 💯
|
||||
</div>
|
||||
|
||||
## Packages
|
||||
|
||||
Run [the Julia REPL](https://docs.julialang.org/en/v1/stdlib/REPL/), enter `]` to bring up Julia's [package manager](https://docs.julialang.org/en/v1/stdlib/Pkg/),
|
||||
and add the listed packages:
|
||||
|
||||
> julia> ]
|
||||
>
|
||||
> (@v1.7) pkg> add #pkg_name#
|
||||
|
||||
- Scientific computing
|
||||
- LinearAlgebra
|
||||
- SparseArrays
|
||||
- Statistics : mean
|
||||
- Image process
|
||||
- ImageFiltering: imfilter
|
||||
- Modelling
|
||||
- [Gmsh](https://onlinelibrary.wiley.com/doi/10.1002/nme.2579)
|
||||
- FEM
|
||||
- [Gridap](https://joss.theoj.org/papers/10.21105/joss.02520)
|
||||
- AD
|
||||
- [ForwardDiff](https://arxiv.org/abs/1607.07892)
|
||||
- [Zygote](https://arxiv.org/abs/1810.07951)
|
||||
- Optimization
|
||||
- [NLopt](https://github.com/stevengj/nlopt)
|
||||
- Visualization
|
||||
- Plots
|
||||
|
||||
## TODOs
|
||||
|
||||
- [x] [`top99neo.jl`](./top_oc/top99neo.jl)
|
||||
top99neo.m rewritten in Julia
|
||||
- [x] [`MMA.jl`](./top_mma/MMA.jl)
|
||||
MMA algorithm (mmasub.m + subsolve.m) rewritten in Julia
|
||||
- [x] [`top99neo_mma.jl`](./top_mma/top99neo_mma.jl)
|
||||
2D code (top99neo + MMA) written in Julia
|
||||
- [ ] `top99neo_AD.jl`
|
||||
Sensitivity Analysis using Automatic Differentiation
|
||||
- [ ] `top99neo_NLopt.jl`
|
||||
Optimization solved with NLopt
|
||||
- [ ] `top3D.jl`
|
||||
3D code (top3D125 + MMA) written in Julia
|
||||
- [ ] `top_flux.jl`
|
||||
Combine TO with machine learning through [Flux](https://arxiv.org/abs/1811.01457)
|
||||
|
||||
## Acknowledgements
|
||||
|
||||
- [TopOpt Group](https://www.topopt.mek.dtu.dk/) 🇩🇰
|
||||
Matlab codes for topology optimization
|
||||
|
||||
- v1.0 [**top99.m**](https://www.topopt.mek.dtu.dk/Apps-and-software/A-99-line-topology-optimization-code-written-in-MATLAB)
|
||||
Educatianal TO enlightenment for every beginners
|
||||
- v2.0 [**top88.m**](https://www.topopt.mek.dtu.dk/Apps-and-software/Efficient-topology-optimization-in-MATLAB)
|
||||
Loop vectorization and memory preallocation
|
||||
- v3.0 [**top99neo.m**](https://www.topopt.mek.dtu.dk/Apps-and-software/New-99-line-topology-optimization-code-written-in-MATLAB)
|
||||
Half matrix assembly operation, filter implementation and volume-preserving density projection
|
||||
|
||||
- Prof. [Krister Svanberg](https://people.kth.se/~krille/) 🇸🇪
|
||||
- Freely available [Matlab code](http://www.smoptit.se/) for CCSA/MMA[^4] and GCMMA
|
||||
|
||||
## Author ©️
|
||||
|
||||
📧 Please contact to liyu_npu@outlook.com
|
||||
|
||||
⚠️ Disclaimer: The author reserves all rights but does not guarantee that the code is free from errors. Furthermore, we shall not be liable in any event.
|
||||
|
||||
| Name | Info. | Hobby | Food |
|
||||
| ----- | :--------: | :-----------: | :-------: |
|
||||
| Yu Li | 🇨🇳 🎓 1️⃣9️⃣9️⃣0️⃣ ♑ | 🎧 🃏 🎮 🏀 🏊 🏃 🚴♂️ | 🍦 🦞 🍣 🌽 🍌 |
|
||||
|
||||
```bibtex
|
||||
@misc{Yu2022,
|
||||
author = {Yu Li},
|
||||
title = {Topology Optimization in Julia},
|
||||
year = {2022},
|
||||
publisher = {GitHub},
|
||||
journal = {GitHub repository},
|
||||
howpublished = {\url{https://github.com/yuloveyet/Topology-Optimization-in-Julia}},
|
||||
}
|
||||
```
|
||||
|
||||
**References**
|
||||
[^1]: Ferrari, F., & Sigmund, O. (2020). A new generation 99 line Matlab code for compliance topology optimization and its extension to 3D. Structural and Multidisciplinary Optimization, 62(4), 2211-2228.
|
||||
[^2]:Sigmund, O. (2001). A 99 line topology optimization code written in Matlab. Structural and multidisciplinary optimization, 21(2), 120-127.
|
||||
[^3]:Andreassen, E., Clausen, A., Schevenels, M., Lazarov, B. S., & Sigmund, O. (2011). Efficient topology optimization in MATLAB using 88 lines of code. Structural and Multidisciplinary Optimization, 43(1), 1-16.
|
||||
[^4]: Svanberg, K. (2002). A class of globally convergent optimization methods based on conservative convex separable approximations. SIAM journal on optimization, 12(2), 555-573.
|
|
@ -0,0 +1,407 @@
|
|||
########################################################################################################
|
||||
### MMA-Julia (1.7.1)
|
||||
########################################################################################################
|
||||
### MODIFIED BY YULI 2022/5/4
|
||||
########################################################################################################
|
||||
export mmasub
|
||||
# Loading modules
|
||||
using LinearAlgebra
|
||||
using SparseArrays
|
||||
|
||||
########################################################################################################
|
||||
### MMA FUNCTIONS ###
|
||||
########################################################################################################
|
||||
|
||||
# Function for the MMA sub problem
|
||||
function mmasub(m::Int, n::Int, iter::Int, xval::Array{Float64}, xmin::Array{Float64},
|
||||
xmax::Array{Float64}, xold1::Array{Float64}, xold2::Array{Float64}, f0val,
|
||||
df0dx::Array{Float64}, df0dx2:: Array{Float64},
|
||||
fval::Float64, dfdx::Array{Float64}, dfdx2::Array{Float64},
|
||||
low::Array{Float64}, upp::Array{Float64}, a0::Float64,
|
||||
a::Array{Float64}, c::Array{Float64}, d::Array{Float64})
|
||||
|
||||
# """
|
||||
# This function mmasub performs one MMA-iteration, aimed at solving the nonlinear programming problem:
|
||||
#
|
||||
# Minimize f_0(x) + a_0*z + sum( c_i*y_i + 0.5*d_i*(y_i)^2 )
|
||||
# subject to f_i(x) - a_i*z - y_i < = 0, i = 1, ..., m
|
||||
# xmin_j < = x_j < = xmax_j, j = 1, ..., n
|
||||
# z > = 0, y_i > = 0, i = 1, ..., m
|
||||
# INPUT:
|
||||
#
|
||||
# m = The number of general constraints.
|
||||
# n = The number of variables x_j.
|
||||
# iter = Current iteration number ( =1 the first time mmasub is called).
|
||||
# xval = Column vector with the current values of the variables x_j.
|
||||
# xmin = Column vector with the lower bounds for the variables x_j.
|
||||
# xmax = Column vector with the upper bounds for the variables x_j.
|
||||
# xold1 = xval, one iteration ago (provided that iter>1).
|
||||
# xold2 = xval, two iterations ago (provided that iter>2).
|
||||
# f0val = The value of the objective function f_0 at xval.
|
||||
# df0dx = Column vector with the derivatives of the objective function
|
||||
# f_0 with respect to the variables x_j, calculated at xval.
|
||||
# fval = Column vector with the values of the constraint functions f_i, calculated at xval.
|
||||
# dfdx = (m x n)-matrix with the derivatives of the constraint functions
|
||||
# f_i with respect to the variables x_j, calculated at xval.
|
||||
# dfdx(i,j) = the derivative of f_i with respect to x_j.
|
||||
# low = Column vector with the lower asymptotes from the previous
|
||||
# iteration (provided that iter>1).
|
||||
# upp = Column vector with the upper asymptotes from the previous
|
||||
# iteration (provided that iter>1).
|
||||
# a0 = The constants a_0 in the term a_0*z.
|
||||
# a = Column vector with the constants a_i in the terms a_i*z.
|
||||
# c = Column vector with the constants c_i in the terms c_i*y_i.
|
||||
# d = Column vector with the constants d_i in the terms 0.5*d_i*(y_i)^2.
|
||||
#
|
||||
# OUTPUT:
|
||||
#
|
||||
# xmma = Column vector with the optimal values of the variables x_j
|
||||
# in the current MMA subproblem.
|
||||
# ymma = Column vector with the optimal values of the variables y_i
|
||||
# in the current MMA subproblem.
|
||||
# zmma = Scalar with the optimal value of the variable z
|
||||
# in the current MMA subproblem.
|
||||
# lam = Lagrange multipliers for the m general MMA constraints.
|
||||
# xsi = Lagrange multipliers for the n constraints alfa_j - x_j < = 0.
|
||||
# eta = Lagrange multipliers for the n constraints x_j - beta_j < = 0.
|
||||
# mu = Lagrange multipliers for the m constraints -y_i < = 0.
|
||||
# zet = Lagrange multiplier for the single constraint -z < = 0.
|
||||
# s = Slack variables for the m general MMA constraints.
|
||||
# low = Column vector with the lower asymptotes, calculated and used
|
||||
# in the current MMA subproblem.
|
||||
# upp = Column vector with the upper asymptotes, calculated and used
|
||||
# in the current MMA subproblem.
|
||||
# """
|
||||
|
||||
epsimin = sqrt(m + n) * 10^(-9)
|
||||
feps = 0.000001
|
||||
asyinit = 0.5
|
||||
asyincr = 1.05
|
||||
asydecr = 0.65
|
||||
albefa = 0.1
|
||||
een = ones(n, 1)
|
||||
zeron = zeros(n, 1)
|
||||
|
||||
|
||||
if iter <= 2
|
||||
low = xval - asyinit .* (xmax - xmin)
|
||||
upp = xval + asyinit .* (xmax - xmin)
|
||||
else
|
||||
zzz = (xval - xold1) .* (xold1 - xold2)
|
||||
factor = copy(een)
|
||||
factor[findall(zzz .> 0.0)] .= asyincr
|
||||
factor[findall(zzz .< 0.0)] .= asydecr
|
||||
low = xval - factor .* (xold1 - low)
|
||||
upp = xval + factor .* (upp - xold1)
|
||||
lowmin = xval - 10.0 .* (xmax - xmin)
|
||||
lowmax = xval - 0.01 .* (xmax - xmin)
|
||||
uppmin = xval + 0.01 .* (xmax - xmin)
|
||||
uppmax = xval + 10.0 .* (xmax - xmin)
|
||||
low = max.(low, lowmin)
|
||||
low = min.(low, lowmax)
|
||||
upp = min.(upp, uppmax)
|
||||
upp = max.(upp, uppmin)
|
||||
end
|
||||
|
||||
zzz = low + albefa .* (xval - low)
|
||||
alfa = max.(zzz, xmin)
|
||||
|
||||
zzz = upp - albefa .* (upp - xval)
|
||||
beta = min.(zzz, xmax)
|
||||
|
||||
|
||||
ux1 = upp - xval
|
||||
ux2 = ux1 .* ux1
|
||||
ux3 = ux2 .* ux1
|
||||
xl1 = xval - low
|
||||
xl2 = xl1 .* xl1
|
||||
xl3 = xl2 .* xl1
|
||||
ul1 = upp - low
|
||||
ulinv1 = een ./ ul1
|
||||
|
||||
uxinv1 = een ./ ux1
|
||||
xlinv1 = een ./ xl1
|
||||
uxinv3 = een ./ ux3
|
||||
xlinv3 = een ./ xl3
|
||||
diap = (ux3 .* xl1) ./ (2 * ul1)
|
||||
diaq = (ux1 .* xl3) ./ (2 * ul1)
|
||||
|
||||
p0 = copy(zeron)
|
||||
q0 = copy(zeron)
|
||||
p0[findall(df0dx .> 0.0)] = df0dx[findall(df0dx .> 0.0)]
|
||||
p0 = p0 + 0.001 .* abs.(df0dx) + feps * ulinv1
|
||||
p0 = p0 .* ux2
|
||||
# q0 = zeron
|
||||
q0[findall(df0dx .< 0.0)] = -df0dx[findall(df0dx .< 0.0)]
|
||||
q0 = q0 + 0.001 * abs.(df0dx) + feps * ulinv1
|
||||
q0 = q0 .* xl2
|
||||
dg0dx2 = 2 * (p0 ./ ux3 + q0 ./ xl3)
|
||||
del0 = df0dx2 - dg0dx2
|
||||
delpos0 = copy(zeron)
|
||||
delpos0[findall(del0 .> 0.0)] = del0[findall(del0 .> 0.0)]
|
||||
p0 = p0 + delpos0 .* diap
|
||||
q0 = q0 + delpos0 .* diaq
|
||||
P = spzeros(m, n)
|
||||
P[findall(dfdx .> 0.0)] = dfdx[findall(dfdx .> 0.0)]
|
||||
# P = P * diag(ux2);
|
||||
# P = P * spdiagm(n, n, 0 => ux2)
|
||||
P = P * sparse(Diagonal(ux2[:]))
|
||||
Q = spzeros(m, n)
|
||||
Q[findall(dfdx .< 0.0)] = -dfdx[findall(dfdx .< 0.0)]
|
||||
# Q = Q * diag(xl2);
|
||||
# Q = Q * spdiagm(n, n, 0 => xl2)
|
||||
Q = Q * sparse(Diagonal(xl2[:]))
|
||||
# dgdx2 = 2.0*(P*diag(uxinv3) + Q*diag(xlinv3));
|
||||
# dgdx2 = P * spdiagm(n, n, 0 => uxinv3) + Q * spdiagm(n, n, 0 => xlinv3)
|
||||
dgdx2 = P * sparse(Diagonal(uxinv3[:])) + Q * sparse(Diagonal(xlinv3[:]))
|
||||
dgdx2 = 2.0 * dgdx2
|
||||
del = dfdx2 - dgdx2
|
||||
delpos = zeros(m, n)
|
||||
delpos[findall(del .> 0.0)] = del[findall(del .> 0.0)]
|
||||
# P = P + delpos*diag(diap);
|
||||
# P = P + delpos * spdiagm(n, n, 0 => diap)
|
||||
P = P + delpos * sparse(Diagonal(diap[:]))
|
||||
# Q = Q + delpos*diag(diaq);
|
||||
# Q = Q + delpos * spdiagm(n, n, 0 => diap)
|
||||
Q = Q + delpos * sparse(Diagonal(diap[:]))
|
||||
b = P * uxinv1 + Q * xlinv1 .- fval
|
||||
|
||||
|
||||
# b = P * uxinv + Q * xlinv - fval
|
||||
# Solving the subproblem by a primal-dual Newton method
|
||||
xmma, ymma, zmma, lam, xsi, eta, mu, zet, s =
|
||||
subsolv(m, n, epsimin, low, upp, alfa, beta, p0, q0, P, Q, a0, a, b, c, d)
|
||||
# Return values
|
||||
return xmma, ymma, zmma, lam, xsi, eta, mu, zet, s, low, upp
|
||||
end
|
||||
|
||||
|
||||
# Function for solving the subproblem (can be used for MMA and GCMMA)
|
||||
function subsolv(m::Int, n::Int, epsimin::Float64, low::Array{Float64}, upp::Array{Float64},
|
||||
alfa::Array{Float64}, beta::Array{Float64}, p0::Array{Float64}, q0::Array{Float64},
|
||||
P::Array{Float64}, Q::Array{Float64}, a0::Float64, a::Array{Float64}, b::Array{Float64},
|
||||
c::Array{Float64}, d::Array{Float64})
|
||||
|
||||
# """
|
||||
# This function subsolv solves the MMA subproblem:
|
||||
#
|
||||
# minimize SUM[p0j/(uppj-xj) + q0j/(xj-lowj)] + a0*z + SUM[ci*yi + 0.5*di*(yi)^2],
|
||||
#
|
||||
# subject to SUM[pij/(uppj-xj) + qij/(xj-lowj)] - ai*z - yi < = bi,
|
||||
# alfaj < = xj < = betaj, yi > = 0, z > = 0.
|
||||
#
|
||||
# Input : m, n, low, upp, alfa, beta, p0, q0, P, Q, a0, a, b, c, d.
|
||||
# Output: xmma, ymma, zmma, slack variables and Lagrange multiplers.
|
||||
# """
|
||||
|
||||
een = ones(Float64, n)
|
||||
eem = ones(Float64, m)
|
||||
epsi = 1.0
|
||||
epsvecn = epsi .* een
|
||||
epsvecm = epsi .* eem
|
||||
x = 0.5 .* (alfa + beta)
|
||||
y = copy(eem)
|
||||
z = 1.0
|
||||
lam = copy(eem)
|
||||
xsi = een ./ (x - alfa)
|
||||
xsi = max.(xsi, een)
|
||||
eta = een ./ (beta - x)
|
||||
eta = max.(eta, een)
|
||||
mu = max.(eem, 0.5 .* c)
|
||||
zet = 1.0
|
||||
s = copy(eem)
|
||||
itera = 0
|
||||
|
||||
while epsi > epsimin # Start while epsi>epsimin
|
||||
epsvecn = epsi .* een
|
||||
epsvecm = epsi .* eem
|
||||
ux1 = upp - x
|
||||
xl1 = x - low
|
||||
ux2 = ux1 .* ux1
|
||||
xl2 = xl1 .* xl1
|
||||
uxinv1 = een ./ ux1
|
||||
xlinv1 = een ./ xl1
|
||||
|
||||
plam = p0 + (P)' * lam
|
||||
qlam = q0 + (Q)' * lam
|
||||
gvec = P * uxinv1 + Q * xlinv1
|
||||
dpsidx = (plam ./ ux2) - (qlam ./ xl2)
|
||||
|
||||
rex = dpsidx - xsi + eta
|
||||
rey = c + d .* y - mu - lam
|
||||
rez = a0 .- zet .- (a)' * lam
|
||||
relam = gvec - a .* z - y + s - b
|
||||
rexsi = xsi .* (x - alfa) - epsvecn
|
||||
reeta = eta .* (beta - x) - epsvecn
|
||||
remu = mu .* y - epsvecm
|
||||
rezet = zet * z - epsi
|
||||
res = lam .* s - epsvecm
|
||||
|
||||
residu1 = [rex' rey' rez]'
|
||||
residu2 = [relam' rexsi' reeta' remu' rezet res']'
|
||||
residu = [residu1' residu2']'
|
||||
residunorm = norm(residu, 2)
|
||||
residumax = maximum(abs.(residu))
|
||||
|
||||
ittt = 0
|
||||
while (residumax > 0.9 * epsi) && (ittt < 100) # Start while (residumax>0.9*epsi) and (ittt<100)
|
||||
ittt = ittt + 1
|
||||
itera = itera + 1
|
||||
|
||||
if ittt == 100
|
||||
println("max inner iter reached")
|
||||
end
|
||||
ux1 = upp - x
|
||||
xl1 = x - low
|
||||
ux2 = ux1 .* ux1
|
||||
xl2 = xl1 .* xl1
|
||||
ux3 = ux1 .* ux2
|
||||
xl3 = xl1 .* xl2
|
||||
uxinv1 = een ./ ux1
|
||||
xlinv1 = een ./ xl1
|
||||
uxinv2 = een ./ ux2
|
||||
xlinv2 = een ./ xl2
|
||||
plam = p0 + (P)' * lam
|
||||
qlam = q0 + (Q)' * lam
|
||||
gvec = P * uxinv1 + Q * xlinv1
|
||||
# GG = P .* transpose(uxinv2) - Q .* transpose(xlinv2)
|
||||
# GG = P .* spdiagm(n,n, 0 => uxinv2) - Q .* spdiagm(n, n, 0 => xlinv2)
|
||||
GG = P * sparse(Diagonal(uxinv2[:])) - Q * sparse(Diagonal(xlinv2[:]))
|
||||
dpsidx = (plam ./ ux2) - (qlam ./ xl2)
|
||||
delx = dpsidx - epsvecn ./ (x - alfa) + epsvecn ./ (beta - x)
|
||||
dely = c + d .* y - lam - epsvecm ./ y
|
||||
delz = a0 .- transpose(a) * lam .- epsi / z
|
||||
dellam = gvec - a .* z - y - b + epsvecm ./ lam
|
||||
diagx = plam ./ ux3 + qlam ./ xl3
|
||||
diagx = 2.0 .* diagx + xsi ./ (x - alfa) + eta ./ (beta - x)
|
||||
diagxinv = een ./ diagx
|
||||
diagy = d + mu ./ y
|
||||
diagyinv = eem ./ diagy
|
||||
diaglam = s ./ lam
|
||||
diaglamyi = diaglam + diagyinv
|
||||
|
||||
if m < n # Start if m < n
|
||||
blam = dellam .+ dely ./ diagy .- GG * (delx ./ diagx)
|
||||
bb = [blam' delz]'
|
||||
# Alam = spdiagm(0 => diaglamyi) + (GG .* transpose(diagxinv)) * transpose(GG)
|
||||
# Alam = spdiagm(m, m, 0 => diaglamyi) + (GG .* spdiagm(n, n, 0 => diagxinv)) * transpose(GG)
|
||||
Alam = sparse(Diagonal(diaglamyi[:])) + (GG * sparse(Diagonal(diagxinv[:]))) * (GG)'
|
||||
AA = [Alam a
|
||||
a' -zet/z]
|
||||
solut = AA \ bb
|
||||
dlam = solut[1:m]
|
||||
dz = solut[m+1]
|
||||
dx = -delx ./ diagx - ((GG)' * dlam) ./ diagx
|
||||
else
|
||||
diaglamyiinv = eem ./ diaglamyi
|
||||
dellamyi = dellam + dely ./ diagy
|
||||
# Axx = spdiagm(0 => diagx) + (transpose(GG) .* transpose(diaglamyiinv)) * GG
|
||||
# Axx = spdiagm(n,n, 0 => diagx) + (transpose(GG) .* spdiagm(m, m, 0 => diaglamyiinv)) * GG
|
||||
Axx = sparse(Diagonal(diagx[:])) + ((GG)' .* sparse(Diagonal(diaglamyiinv[:]))) * GG
|
||||
azz = zet / z + (a)' * (a ./ diaglamyi)
|
||||
axz = -(GG)' * (a ./ diaglamyi)
|
||||
bx = delx + (GG)' * (dellamyi ./ diaglamyi)
|
||||
bz = delz - (a)' * (dellamyi ./ diaglamyi)
|
||||
# AAr1 = [Axx axz]
|
||||
# AAr2 = [transpose(axz) azz]
|
||||
# AA = [AAr1; AAr2]
|
||||
AA = [Axx axz
|
||||
axz' azz]
|
||||
bb = [-bx' -bz]'
|
||||
solut = AA \ bb
|
||||
dx = solut[1:n]
|
||||
dz = solut[n+1]
|
||||
dlam = (GG * dx) ./ diaglamyi - dz .* (a ./ diaglamyi) + dellamyi ./ diaglamyi
|
||||
end # End if m<n
|
||||
|
||||
dy = -dely ./ diagy + dlam ./ diagy
|
||||
dxsi = -xsi .+ epsvecn ./ (x - alfa) - (xsi .* dx) ./ (x - alfa)
|
||||
deta = -eta .+ epsvecn ./ (beta - x) + (eta .* dx) ./ (beta - x)
|
||||
dmu = -mu .+ epsvecm ./ y - (mu .* dy) ./ y
|
||||
dzet = -zet .+ epsi / z - zet * dz / z
|
||||
ds = -s .+ epsvecm ./ lam - (s .* dlam) ./ lam
|
||||
xx = [y' z lam' xsi' eta' mu' zet s']'
|
||||
dxx = [dy' dz dlam' dxsi' deta' dmu' dzet ds']'
|
||||
#
|
||||
stepxx = -1.01 .* dxx ./ xx
|
||||
stmxx = maximum(stepxx)
|
||||
stepalfa = -1.01 .* dx ./ (x - alfa)
|
||||
stmalfa = maximum(stepalfa)
|
||||
stepbeta = 1.01 .* dx ./ (beta - x)
|
||||
stmbeta = maximum(stepbeta)
|
||||
stmalbe = max.(stmalfa, stmbeta)
|
||||
stmalbexx = max.(stmalbe, stmxx)
|
||||
stminv = max.(stmalbexx, 1.0)
|
||||
steg = 1.0 / stminv
|
||||
#
|
||||
xold = copy(x)
|
||||
yold = copy(y)
|
||||
zold = copy(z)
|
||||
lamold = copy(lam)
|
||||
xsiold = copy(xsi)
|
||||
etaold = copy(eta)
|
||||
muold = copy(mu)
|
||||
zetold = copy(zet)
|
||||
sold = copy(s)
|
||||
#
|
||||
itto = 0
|
||||
resinew = 2.0 * residunorm
|
||||
# Start: while (resinew>residunorm) and (itto<50)
|
||||
while (resinew > residunorm) && (itto < 50)
|
||||
itto = itto + 1
|
||||
|
||||
x = xold + steg .* dx
|
||||
y = yold + steg .* dy
|
||||
z = zold + steg * dz
|
||||
lam = lamold + steg .* dlam
|
||||
xsi = xsiold + steg .* dxsi
|
||||
eta = etaold + steg .* deta
|
||||
mu = muold + steg .* dmu
|
||||
zet = zetold + steg * dzet
|
||||
s = sold + steg .* ds
|
||||
ux1 = upp - x
|
||||
xl1 = x - low
|
||||
ux2 = ux1 .* ux1
|
||||
xl2 = xl1 .* xl1
|
||||
uxinv1 = een ./ ux1
|
||||
xlinv1 = een ./ xl1
|
||||
plam = p0 + (P)' * lam
|
||||
qlam = q0 + (Q)' * lam
|
||||
gvec = P * uxinv1 + Q * xlinv1
|
||||
dpsidx = plam ./ ux2 - qlam ./ xl2
|
||||
|
||||
rex = dpsidx - xsi + eta
|
||||
rey = c .+ d .* y - mu - lam
|
||||
rez = a0 - zet .- (a)' * lam
|
||||
relam = gvec - a .* z - y + s - b
|
||||
rexsi = xsi .* (x - alfa) - epsvecn
|
||||
reeta = eta .* (beta - x) - epsvecn
|
||||
remu = mu .* y - epsvecm
|
||||
rezet = zet * z - epsi
|
||||
res = lam .* s - epsvecm
|
||||
|
||||
residu1 = [rex' rey' rez]'
|
||||
residu2 = [relam' rexsi' reeta' remu' rezet res']'
|
||||
residu = [residu1' residu2']'
|
||||
resinew = norm(residu, 2)
|
||||
steg = steg / 2.0
|
||||
end # End: while (resinew>residunorm) and (itto<50)
|
||||
residunorm = copy(resinew)
|
||||
residumax = maximum(abs.(residu))
|
||||
steg = 2.0 * steg
|
||||
end # End: while (residumax>0.9*epsi) and (ittt<200)
|
||||
epsi = 0.1 * epsi
|
||||
end # End: while epsi>epsimin
|
||||
xmma = copy(x)
|
||||
ymma = copy(y)
|
||||
zmma = copy(z)
|
||||
lamma = lam
|
||||
xsimma = xsi
|
||||
etamma = eta
|
||||
mumma = mu
|
||||
zetmma = zet
|
||||
smma = s
|
||||
# Return values
|
||||
return xmma, ymma, zmma, lamma, xsimma, etamma, mumma, zetmma, smma
|
||||
end
|
||||
|
|
@ -0,0 +1,770 @@
|
|||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# module TopOpt99neo\n",
|
||||
"\n",
|
||||
"# include(\"FastConv.jl/src/FastConv.jl\")\n",
|
||||
"# import .FastConv.fastconv\n",
|
||||
"\n",
|
||||
"using LinearAlgebra, SparseArrays\n",
|
||||
"using Plots\n",
|
||||
"# : heatmap, savefig, @animate\n",
|
||||
"using ImageFiltering: imfilter\n",
|
||||
"using Statistics: mean\n",
|
||||
"\n",
|
||||
"using BenchmarkTools\n",
|
||||
"using Zygote\n",
|
||||
"\n",
|
||||
"BLAS.set_num_threads(1)\n",
|
||||
"\n",
|
||||
"# abstract type top99neo end\n",
|
||||
"\n",
|
||||
"include(\"utils.jl\")\n",
|
||||
"include(\"MMA.jl\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"mutable struct SetUp\n",
|
||||
" nelx::Int\n",
|
||||
" nely::Int\n",
|
||||
"\n",
|
||||
" eta::Float64\n",
|
||||
" beta::Int\n",
|
||||
" betaCnt::NTuple{4,Int}\n",
|
||||
"\n",
|
||||
" move::Float64\n",
|
||||
" maxit::Int\n",
|
||||
" maxchang::Float64\n",
|
||||
"\n",
|
||||
" pasS::Array{Int}\n",
|
||||
" pasV::Array{Int}\n",
|
||||
"\n",
|
||||
" function SetUp()\n",
|
||||
" nelx = 30\n",
|
||||
" nely = 10\n",
|
||||
" # volfrac = 0.3\n",
|
||||
" maxit = 300\n",
|
||||
" move = 0.1\n",
|
||||
" beta = 2\n",
|
||||
" eta = 0.5\n",
|
||||
" maxchang = 1e-6\n",
|
||||
" # penalCnt = {maxit+1, 3, 25, 0.25}\n",
|
||||
" betaCnt = (1, 32, 10, 2) # continuation scheme on beta parCont = { istart,maxPar,steps,deltaPar }\n",
|
||||
" \n",
|
||||
" # elNrs = reshape(1:nely*nelx, nely, nelx)\n",
|
||||
" # a1 = elNrs[Int.(nely/4:nely/2), Int.(nelx/4:nelx/2)]\n",
|
||||
" # pasS, pasV = Array([]), a1[:]\n",
|
||||
" pasS, pasV = Array([]), Array([])\n",
|
||||
" \n",
|
||||
" new(nelx, nely, eta, beta, betaCnt, move, maxit, maxchang, pasS, pasV)\n",
|
||||
" end\n",
|
||||
"\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"setup = SetUp()\n",
|
||||
"# setup.nelx = 20\n",
|
||||
"# setup.nely = 10\n",
|
||||
"# setup.maxit = 100\n",
|
||||
"# typeof(setup.nodeNrs)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"\n",
|
||||
"mutable struct Mat\n",
|
||||
" E0::Float64\n",
|
||||
" Emin::Float64\n",
|
||||
" ν::Float64\n",
|
||||
" penal::Float64\n",
|
||||
" volfrac::Float64\n",
|
||||
"\n",
|
||||
" function Mat()\n",
|
||||
" E0 = 1\n",
|
||||
" Emin = 1e-9\n",
|
||||
" ν = 0.3\n",
|
||||
" penal = 3.0\n",
|
||||
" volfrac = 0.3\n",
|
||||
" new(E0, Emin, ν, penal, volfrac)\n",
|
||||
" end\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"mat = Mat()\n",
|
||||
"# mat.volfrac"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"mutable struct DiscretizationFeature\n",
|
||||
" nEl::Int\n",
|
||||
" nodeNrs::Array{Int,2}\n",
|
||||
" nDof::Int\n",
|
||||
" act\n",
|
||||
"\n",
|
||||
" cMat::Array{Int}\n",
|
||||
" Iar::Array{Int}\n",
|
||||
" Ke::Array{Float64,1}\n",
|
||||
" Ke0::Array{Float64,2}\n",
|
||||
" \n",
|
||||
" function DiscretizationFeature(setup::SetUp, mat::Mat)\n",
|
||||
" nelx = setup.nelx\n",
|
||||
" nely = setup.nely\n",
|
||||
" pasS, pasV = setup.pasS, setup.pasV\n",
|
||||
" ν = mat.ν\n",
|
||||
" \n",
|
||||
" nEl = nely * nelx\n",
|
||||
" nodeNrs = reshape(1:(1+nelx)*(1+nely), nely + 1, nelx + 1)\n",
|
||||
" nDof = (nely + 1) * (nelx + 1) * 2\n",
|
||||
" act = setdiff(collect(1:nEl), union(pasS, pasV))\n",
|
||||
" \n",
|
||||
" cVec = reshape(2 * nodeNrs[1:end-1, 1:end-1] .+ 1, nEl, 1)\n",
|
||||
" # cMat = cVec + Int.([0 1 2 * nely .+ [2 3 0 1] -2 -1])\n",
|
||||
" cMat = Int.(repeat(cVec, 1, 8) + repeat([0 1 2 * nely .+ [2 3 0 1] -2 -1], nelx * nely, 1))\n",
|
||||
" # nDof = (nely + 1) * (nelx + 1) * 2\n",
|
||||
" FuckRow = [1 2 3 4 5 6 7 8]\n",
|
||||
" sI::Array{Int}, sII::Array{Int} = copy(FuckRow), fill(1, 1, 8)\n",
|
||||
" for j in 2:8\n",
|
||||
" sI = cat(sI, FuckRow[j:8]'; dims=2)\n",
|
||||
" # sI = append!(sI, convert(Array{Int},[j:8]))\n",
|
||||
" # sII = cat(2, sII, repmat(j, 1, 8 - j + 1))\n",
|
||||
" sII = cat(sII, fill(j, 1, 8 - j + 1); dims=2)\n",
|
||||
" end\n",
|
||||
" iK::Array{Int,2}, jK::Array{Int,2} = cMat[:, sI][:, 1, :]', cMat[:, sII][:, 1, :]'\n",
|
||||
" Iar = sort([iK[:] jK[:]]; dims=2, rev=true) # comma is a newline\n",
|
||||
" # iK[:], jK[:] .= 0.0, 0.0\n",
|
||||
" c1 = [12, 3, -6, -3, -6, -3, 0, 3, 12, 3, 0, -3, -6, -3, -6, 12, -3, 0, -3, -6, 3, 12, 3, -6, 3, -6, 12, 3, -6, -3, 12, 3, 0, 12, -3, 12]\n",
|
||||
" c2 = [-4, 3, -2, 9, 2, -3, 4, -9, -4, -9, 4, -3, 2, 9, -2, -4, -3, 4, 9, 2, 3, -4, -9, -2, 3, 2, -4, 3, -2, 9, -4, -9, 4, -4, -3, -4]\n",
|
||||
" Ke = 1 / (1 - ν^2) / 24 .* (c1 .+ ν .* c2) # half-KE vector\n",
|
||||
" # full KE\n",
|
||||
" # Ke0::Array{Float64} = zeros(8, 8)\n",
|
||||
" # start_id, end_id = 1, 8\n",
|
||||
" # for i in 1:8\n",
|
||||
" # Ke0[i:8, i] = Ke[start_id:end_id]\n",
|
||||
" # start_id, end_id = end_id + 1, 2 * end_id - start_id\n",
|
||||
" # end\n",
|
||||
" \n",
|
||||
" Ke0::Array{Float64} = zeros(8, 8)\n",
|
||||
" # Index::Array{Int} = [sI' sII']\n",
|
||||
" # Ke0[sI, sII] = Ke\n",
|
||||
" Index = findall(isequal(1), tril(ones(8, 8)))\n",
|
||||
" Ke0[Index] = Ke'\n",
|
||||
" # Ke0 = reshape(Ke0, 8, 8)\n",
|
||||
" Ke0 = Ke0 + Ke0' - diagm(diag(Ke0))\n",
|
||||
" new(nEl, nodeNrs, nDof, act, cMat, Iar, Ke, Ke0)\n",
|
||||
" end\n",
|
||||
"end\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"\n",
|
||||
"disfeature = DiscretizationFeature(setup, mat)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"mutable struct LoadsSupportsBCs\n",
|
||||
" # setup::SetUp\n",
|
||||
" lcDof::Array{Int}\n",
|
||||
" F\n",
|
||||
" free::Array{Int}\n",
|
||||
" fixed::Array{Int}\n",
|
||||
"\n",
|
||||
" function LoadsSupportsBCs(setup::SetUp, disfeature::DiscretizationFeature)\n",
|
||||
" nelx = setup.nelx\n",
|
||||
" nely = setup.nely\n",
|
||||
" nodeNrs = disfeature.nodeNrs\n",
|
||||
" nDof = disfeature.nDof\n",
|
||||
" \n",
|
||||
" load_position::Symbol = :half_MBB\n",
|
||||
" if load_position == :half_MBB\n",
|
||||
" load_nodey, load_nodex = 1, 1\n",
|
||||
" # fixed = union([1:2:2*(nely+1)], 2 * nodeNrs[nely+1, nelx+1])\n",
|
||||
" fixed = union(collect(1:2:2*(nely+1)), 2 * nodeNrs[end, end])\n",
|
||||
" elseif load_position == :cantilever\n",
|
||||
" load_nodey = nely + 1\n",
|
||||
" load_nodex = nelx / 2 + 1\n",
|
||||
" fixed = 1:2*(nely+1)\n",
|
||||
" end\n",
|
||||
" \n",
|
||||
" F = spzeros(nDof)\n",
|
||||
" \n",
|
||||
" load_type::Symbol = :pin\n",
|
||||
" if load_type == :pin # 1 point\n",
|
||||
" lcDof = collect(2 * nodeNrs[load_nodey, load_nodex])\n",
|
||||
" F[2, 1] = -1.0\n",
|
||||
" elseif load_type == :points # 5 points\n",
|
||||
" lcDof = collect(2 * nodeNrs[load_nodey, load_nodex], nodeNrs[load_nodey, load_nodex-1], nodeNrs[load_nodey, load_nodex-2], nodeNrs[load_nodey, load_nodex+1], nodeNrs[load_nodey, load_nodex+2])\n",
|
||||
" F[lcDof', ones(length(lcDof'))] .= -1.0\n",
|
||||
" elseif load_type == :line\n",
|
||||
" lcDof = [2:2*(nely+1):nDof]\n",
|
||||
" F = spzeros(nDof, 1)\n",
|
||||
" # F = sparse(lcDof', ones(length(lcDof')), -1.0)\n",
|
||||
" end\n",
|
||||
" all = collect(1:nDof)\n",
|
||||
" free = setdiff(all, fixed)\n",
|
||||
" \n",
|
||||
" new(lcDof, F, free, fixed)\n",
|
||||
" end\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"load = LoadsSupportsBCs(setup, disfeature)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"mutable struct Initialization \n",
|
||||
" # setup::SetUp\n",
|
||||
" x::Array{Float64}\n",
|
||||
" xPhys::Array{Float64}\n",
|
||||
" xOld::Array{Float64}\n",
|
||||
" ch::Float64\n",
|
||||
" loop::Int\n",
|
||||
" # U::Array{Float64}\n",
|
||||
" # dsK::Array{Float64}\n",
|
||||
" # dV::Array{Float64}\n",
|
||||
"\n",
|
||||
" function Initialization(setup::SetUp, disfeature::DiscretizationFeature, mat::Mat)\n",
|
||||
" pasV = setup.pasV\n",
|
||||
" pasS = setup.pasS\n",
|
||||
" act = disfeature.act\n",
|
||||
" volfrac = mat.volfrac\n",
|
||||
" nEl = disfeature.nEl\n",
|
||||
" nDof = disfeature.nDof\n",
|
||||
" # column vectors\n",
|
||||
" x = zeros(nEl, 1)\n",
|
||||
" # dV[act, 1] .= 1.0 / nEl / volfrac\n",
|
||||
" x[act] .= volfrac\n",
|
||||
" x[pasS] .= 1.0\n",
|
||||
" xPhys, xOld, ch, loop = copy(x), ones(nEl, 1), 1.0, 0\n",
|
||||
" # x̅ x̃\n",
|
||||
" \n",
|
||||
" new(x, xPhys, xOld, ch, loop)\n",
|
||||
" end\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"ini = Initialization(setup, disfeature, mat)\n",
|
||||
"# typeof(ini.x)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"mutable struct Filter \n",
|
||||
" rmin::Float64\n",
|
||||
" ft::Int\n",
|
||||
" h::Array{Float64}\n",
|
||||
" Hs::Array{Float64}\n",
|
||||
" dHs::Array{Float64}\n",
|
||||
"\n",
|
||||
" function Filter(setup::SetUp)\n",
|
||||
" nelx = setup.nelx\n",
|
||||
" nely = setup.nely\n",
|
||||
" # bcF = setup.bcF\n",
|
||||
" rmin = 6.5\n",
|
||||
" ft = 3\n",
|
||||
" dy, dx = meshgrid(-ceil(rmin)+1:ceil(rmin)-1, -ceil(rmin)+1:ceil(rmin)-1)\n",
|
||||
" h = max.(0, rmin .- sqrt.(dx .^ 2 + dy .^ 2))\n",
|
||||
" Hs = imfilter(ones(nely, nelx), h, \"symmetric\")\n",
|
||||
" dHs = Hs\n",
|
||||
" new(rmin, ft, h, Hs, dHs)\n",
|
||||
" end\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"\n",
|
||||
"\n",
|
||||
"filter = Filter(setup)\n",
|
||||
"# filter.Hs"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 15,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"FiniteElementAnalasys (generic function with 1 method)"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"function FiniteElementAnalasys(mat::Mat, disfeature::DiscretizationFeature, load::LoadsSupportsBCs, xPhys::Array{Float64})\n",
|
||||
" nEl, nDof, Iar, Ke = disfeature.nEl, disfeature.nDof, disfeature.Iar, disfeature.Ke\n",
|
||||
" act = disfeature.act\n",
|
||||
" Emin, penal, E0 = mat.Emin, mat.penal, mat.E0\n",
|
||||
" F, free = load.F, load.free\n",
|
||||
"\n",
|
||||
" sK = Emin .+ xPhys .^ penal .* (E0 - Emin)\n",
|
||||
" sK = reshape(Ke[:] * sK', length(Ke) * nEl, 1)\n",
|
||||
" K = sparse(Iar[:, 1], Iar[:, 2], vec(sK), nDof, nDof)\n",
|
||||
" U = zeros(nDof)\n",
|
||||
" #~\n",
|
||||
" # U[free] = cholesky(Symmetric(K[free, free], :L), check=false) \\ F[free]\n",
|
||||
" K0 = K + K' - diagm(diag(K))\n",
|
||||
" U[free] = K0[free, free] \\ F[free]\n",
|
||||
" #~\n",
|
||||
" Obj = F' * U\n",
|
||||
" Vf = mean(xPhys[act])\n",
|
||||
" return U, Obj, Vf\n",
|
||||
"end\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"1-element Vector{Float64}:\n",
|
||||
" 4558.123953305818"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"U, C, vf = FiniteElementAnalasys(mat, disfeature, load, xval)\n",
|
||||
"C"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"function SensitivityAnalasys(setup::SetUp, filter::Filter, mat::Mat, disfeature::DiscretizationFeature, \n",
|
||||
" U::Array{Float64}, xPhys::Array{Float64})\n",
|
||||
" nelx, nely, act = setup.nelx, setup.nely, disfeature.act\n",
|
||||
" dHs, h = filter.dHs, filter.h\n",
|
||||
" E0, Emin, penal = mat.E0, mat.Emin, mat.penal\n",
|
||||
" nEl, cMat, Ke0 = disfeature.nEl, disfeature.cMat, disfeature.Ke0\n",
|
||||
" act = disfeature.act\n",
|
||||
"\n",
|
||||
" dsK, dV = zeros(nEl, 1), zeros(nEl, 1)\n",
|
||||
" dV[act, 1] .= 1.0 / length(act)\n",
|
||||
" dsK[act] = -penal * (E0 - Emin) .* xPhys[act] .^ (penal - 1)\n",
|
||||
"\n",
|
||||
" dc = dsK .* sum((U[cMat] * Ke0) .* U[cMat], dims=2)\n",
|
||||
" dc = imfilter(reshape(dc, nely, nelx) ./ dHs, h, \"symmetric\")\n",
|
||||
" dV0 = imfilter(reshape(dV, nely, nelx) ./ dHs, h, \"symmetric\")\n",
|
||||
"\n",
|
||||
" return reshape(dc, nEl, 1), reshape(dV0, nEl, 1)\n",
|
||||
"end\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"function AutomaticDifferentiation(mat::Mat, disfeature::DiscretizationFeature, load::LoadsSupportsBCs, xPhys::Array{Float64})\n",
|
||||
" # fea = FiniteElementAnalasys()\n",
|
||||
" dc_AD = gradient(x -> FiniteElementAnalasys(mat, disfeature, load, x), xPhys)\n",
|
||||
" # dv_AD = gradient(x)\n",
|
||||
"\n",
|
||||
" return dc_AD\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"xval = ini.x\n",
|
||||
"dc_AD = AutomaticDifferentiation(mat, disfeature, load, xval)\n",
|
||||
"typeof(dc_AD)\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"typeof(dc_AD)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"function Visualization(setup::SetUp, x::Array{Float64}, loop::Int)\n",
|
||||
" nelx, nely = setup.nelx, setup.nely\n",
|
||||
" # cmap = cgrad(:Blues_9, rev=false)\n",
|
||||
" plot = heatmap(reshape(x, nely, nelx), c=:Blues_9, aspect_ratio=:equal, yflip=true, grid=false, axis=:off, tick=false, colorbar=false, border=nothing, dpi=300, size=(400,nely/nelx*400) , legend=:none, display_type=:gui)\n",
|
||||
" display(plot)\n",
|
||||
" savefig(plot, \"./top/res_$loop.pdf\")\n",
|
||||
" # PLOT FINAL DESIGN\n",
|
||||
" # heatmap(1.0 .- x[end:-1:1, :], yaxis=false, xaxis=false, legend=:none,color=:greys, grid=false, border=nothing, aspect_ratio=:equal)\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"function Optimization(setup::SetUp, mat::Mat, load::LoadsSupportsBCs,\n",
|
||||
" filter::Filter, ini::Initialization, disfeature::DiscretizationFeature)\n",
|
||||
" ch, loop = ini.ch, ini.loop\n",
|
||||
" x, xPhys, xOld = ini.x, ini.xPhys, ini.xOld\n",
|
||||
"\n",
|
||||
" maxit, maxchang = setup.maxit, setup.maxchang\n",
|
||||
" nely, nelx = setup.nely, setup.nelx\n",
|
||||
" eta, beta, betaCnt = setup.eta, setup.beta, setup.betaCnt\n",
|
||||
" volfrac, penal, nEl = mat.volfrac, mat.penal, disfeature.nEl\n",
|
||||
" act = disfeature.act\n",
|
||||
"\n",
|
||||
" Hs, h, ft = filter.Hs, filter.h, filter.ft\n",
|
||||
"\n",
|
||||
" opt_hist = []\n",
|
||||
" vf_hist = []\n",
|
||||
"\n",
|
||||
" ### Initiation of MMA ###\n",
|
||||
" xval = copy(x[act])\n",
|
||||
" xold1 = copy(xval)\n",
|
||||
" xold2 = copy(xold1)\n",
|
||||
"\n",
|
||||
" low = copy(xval)\n",
|
||||
" upp = copy(low)\n",
|
||||
" #########################\n",
|
||||
"\n",
|
||||
" anim = @animate while ch > maxchang && loop < maxit || beta < betaCnt[2]\n",
|
||||
" @time begin\n",
|
||||
" loop = loop + 1\n",
|
||||
" # COMPUTE PHYSICAL DENSITY FIELD \n",
|
||||
" xTilde = imfilter(reshape(x, nely, nelx), h, \"symmetric\") ./ Hs\n",
|
||||
" xPhys[act] = copy(xTilde[act])\n",
|
||||
" if ft > 1\n",
|
||||
" f = (mean(prj(xPhys[act], eta, beta)) .- volfrac) * (ft == 3)\n",
|
||||
" while abs(f) > maxchang\n",
|
||||
" eta = eta - f / mean(deta(xPhys[:], eta, beta))\n",
|
||||
" f = mean(prj(xPhys[act], eta, beta)) - volfrac\n",
|
||||
" end\n",
|
||||
" filter.dHs = Hs ./ reshape(dprj(xTilde, eta, beta), nely, nelx)\n",
|
||||
" xPhys = prj(xPhys, eta, beta)\n",
|
||||
" end\n",
|
||||
" ch = norm(xPhys - xOld) ./ sqrt(nEl)\n",
|
||||
" xOld = copy(xPhys)\n",
|
||||
" #~ SETUP AND SOLVE EQUILIBRIUM EQUATIONS\n",
|
||||
" U, C, Vf = FiniteElementAnalasys(mat, disfeature, load, xPhys)\n",
|
||||
" push!(opt_hist, C)\n",
|
||||
" push!(vf_hist, Vf)\n",
|
||||
" #~ COMPUTE SENSITIVITIES\n",
|
||||
" dc, dV0 = SensitivityAnalasys(setup, filter, mat, disfeature, U, xPhys)\n",
|
||||
" #~ MMA iteration\n",
|
||||
" xmma, low, upp = MMAupdate(xval, low, upp, xold1, xold2, C, Vf, dc, dV0, loop, Mat().volfrac)\n",
|
||||
" # Some vectors are updated:\n",
|
||||
" xold2 = copy(xold1)\n",
|
||||
" xold1 = copy(xval)\n",
|
||||
" xval = copy(xmma)\n",
|
||||
"\n",
|
||||
" x[act] = xval\n",
|
||||
"\n",
|
||||
" #~ CONTINUATION\n",
|
||||
" beta = cnt(beta, betaCnt, loop, ch, maxchang)\n",
|
||||
"\n",
|
||||
" heatmap(reshape(xPhys, nely, nelx), c=:Blues_9, aspect_ratio=:equal, yflip=true, grid=false, axis=:off, tick=false, colorbar=false, border=nothing, dpi=300, size=(400, nely / nelx * 400), legend=:none)\n",
|
||||
"\n",
|
||||
" end\n",
|
||||
" if mod(loop, 10) == 0\n",
|
||||
" println(\"It.: $loop C.: $C Vf.: $Vf ch.: $ch, p.: $penal beta.:$beta eta.: $eta \")\n",
|
||||
" Visualization(setup, xPhys, loop)\n",
|
||||
" end\n",
|
||||
" end\n",
|
||||
" gif(anim, \"./top/top_hist.gif\", fps=8)\n",
|
||||
" return xPhys, opt_hist, vf_hist, anim\n",
|
||||
"end\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# @time xPhys, opt_hist, vf_hist, loop = Optimization(setup, mat, load, filter, ini, disfeature)\n",
|
||||
"Optimization(setup, mat, load, filter, ini, disfeature)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"function MMAupdate(xval::Array{Float64}, low::Array{Float64}, upp::Array{Float64}, \n",
|
||||
" xold1::Array{Float64}, xold2::Array{Float64},\n",
|
||||
" Obj::Float64, Vf::Float64, dc::Array{Float64}, dV0::Array{Float64}, loop::Int, volfrac::Float64)\n",
|
||||
" move = 0.1\n",
|
||||
" ### Initiation of MMA ###\n",
|
||||
" m = 1\n",
|
||||
" # active design variable\n",
|
||||
" n = length(act)\n",
|
||||
" onen = ones(n, 1)\n",
|
||||
" onem = ones(m, 1)\n",
|
||||
" zeron = zeros(n, 1)\n",
|
||||
" zerom = zeros(m, 1)\n",
|
||||
" a_mma = zerom\n",
|
||||
" c_mma = 1.0e3 * onem\n",
|
||||
" d_mma = zerom\n",
|
||||
" a0 = 1.0\n",
|
||||
" # column vector\n",
|
||||
" # xval = xval\n",
|
||||
" xmin = max.(xval .- move, zeron)\n",
|
||||
" xmax = min.(xval .+ move, onen)\n",
|
||||
"\n",
|
||||
" # low = low\n",
|
||||
" # upp = upp\n",
|
||||
" # objective function \n",
|
||||
" f0val = Obj\n",
|
||||
" df0dx = dc[act]\n",
|
||||
" df0dx2 = 0.0 * df0dx\n",
|
||||
" # constraint function\n",
|
||||
" fval = Vf / volfrac - 1.0 # column vector\n",
|
||||
" dfdx = reshape(dV0[act], 1, length(act)) ./ volfrac # (m * n)\n",
|
||||
" dfdx2 = 0.0 * dfdx\n",
|
||||
"\n",
|
||||
" # The MMA subproblem is solved at the point xval:\n",
|
||||
" xmma, ymma, zmma, lam, xsi, eta, mu, zet, s, low, upp =\n",
|
||||
" mmasub(m, n, loop, xval, xmin, xmax, xold1, xold2,\n",
|
||||
" f0val, df0dx, df0dx2, fval, dfdx, dfdx2, low, upp, a0, a_mma, c_mma, d_mma)\n",
|
||||
" \n",
|
||||
" return xmma, low, upp\n",
|
||||
" # xmma =\n",
|
||||
" # mmasub(m, n, loop, xval, xmin, xmax, xold1, xold2,\n",
|
||||
" # f0val, df0dx, df0dx2, fval, dfdx, dfdx2, low, upp, a0, a_mma, c_mma, d_mma)\n",
|
||||
" # return xmma\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"include(\"MMA.jl\")\n",
|
||||
"### Initiation of MMA ###\n",
|
||||
"x = ones(disfeature.nEl, 1) * 0.5\n",
|
||||
"act = disfeature.act\n",
|
||||
"xval = copy(x[act])\n",
|
||||
"xold1 = copy(xval)\n",
|
||||
"xold2 = copy(xold1)\n",
|
||||
"\n",
|
||||
"low = copy(xval)\n",
|
||||
"upp = copy(low)\n",
|
||||
"\n",
|
||||
"xval, low, upp = MMAupdate(xval, low, upp, xold1, xold2, 1.0, 0.3, x[act],\n",
|
||||
" reshape(x[act], 1, length(act)) , 1, Mat().volfrac)\n",
|
||||
"# ux2 = MMAupdate(xval, low, upp, xold1, xold2, 1.0, 0.3, x[act],\n",
|
||||
"# reshape(x[act], 1, length(act)) , 1, Mat().volfrac)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"size(xval)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"using NLopt\n",
|
||||
"\n",
|
||||
"function myfunc(x::Vector, grad::Vector)\n",
|
||||
" if length(grad) > 0\n",
|
||||
" grad[1] = 0\n",
|
||||
" grad[2] = 0.5/sqrt(x[2])\n",
|
||||
" end\n",
|
||||
" return sqrt(x[2])\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"function myconstraint(x::Vector, grad::Vector, a, b)\n",
|
||||
" if length(grad) > 0\n",
|
||||
" grad[1] = 3a * (a*x[1] + b)^2\n",
|
||||
" grad[2] = -1\n",
|
||||
" end\n",
|
||||
" (a*x[1] + b)^3 - x[2]\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"opt = Opt(:LD_MMA, 2)\n",
|
||||
"opt.lower_bounds = [-Inf, 0.]\n",
|
||||
"opt.xtol_rel = 1e-4\n",
|
||||
"\n",
|
||||
"opt.min_objective = myfunc\n",
|
||||
"inequality_constraint!(opt, (x,g) -> myconstraint(x,g,2,0), 1e-8)\n",
|
||||
"inequality_constraint!(opt, (x,g) -> myconstraint(x,g,-1,1), 1e-8)\n",
|
||||
"\n",
|
||||
"(minf,minx,ret) = optimize(opt, [1.234, 5.678])\n",
|
||||
"numevals = opt.numevals # the number of function evaluations\n",
|
||||
"println(\"got $minf at $minx after $numevals iterations (returned $ret)\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"using NLopt\n",
|
||||
"function update_NLopt(act, xval, low, upp, C, Vf, dc, dV0)\n",
|
||||
" alg = \":LD_MMA\"\n",
|
||||
" n = length(act)\n",
|
||||
" opt = Opt(alg, n)\n",
|
||||
" opt.lower_bounds = low\n",
|
||||
" opt.upper_bounds = upp\n",
|
||||
"\n",
|
||||
" opt.min_objective = ObjectiveFunction\n",
|
||||
"\n",
|
||||
" inequality_constraint!(opt, (x, g) -> VolumeConstraint(x, g;), 1e-8)\n",
|
||||
"\n",
|
||||
" return\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"function ObjectiveFunction(x, grad; act, obj, dobj)\n",
|
||||
" if length(grad) > 0\n",
|
||||
" grad[:] = dobj[act]\n",
|
||||
" end\n",
|
||||
" return obj\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"function VolumeConstraint(x, g; )\n",
|
||||
" if length(grad) > 0\n",
|
||||
" grad[:] = dV0\n",
|
||||
" end\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"f(x::Matrix, y) = sum(sin, x) + prod(tan, x) * sum(sqrt, x)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"using Zygote\n",
|
||||
"grad = gradient(y -> f(y), rand(3))\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"typeof(grad)"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Julia 1.7.2",
|
||||
"language": "julia",
|
||||
"name": "julia-1.7"
|
||||
},
|
||||
"language_info": {
|
||||
"file_extension": ".jl",
|
||||
"mimetype": "application/julia",
|
||||
"name": "julia",
|
||||
"version": "1.7.2"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 4
|
||||
}
|
After Width: | Height: | Size: 8.2 MiB |
After Width: | Height: | Size: 12 KiB |
After Width: | Height: | Size: 12 KiB |
After Width: | Height: | Size: 13 KiB |
After Width: | Height: | Size: 12 KiB |
After Width: | Height: | Size: 12 KiB |
After Width: | Height: | Size: 9.7 KiB |
After Width: | Height: | Size: 12 KiB |
After Width: | Height: | Size: 13 KiB |
After Width: | Height: | Size: 12 KiB |
After Width: | Height: | Size: 11 KiB |
After Width: | Height: | Size: 10 KiB |
After Width: | Height: | Size: 9.9 KiB |
After Width: | Height: | Size: 9.4 KiB |
After Width: | Height: | Size: 9.0 KiB |
After Width: | Height: | Size: 8.7 KiB |
After Width: | Height: | Size: 8.5 KiB |
After Width: | Height: | Size: 10 KiB |
After Width: | Height: | Size: 8.2 KiB |
After Width: | Height: | Size: 7.9 KiB |
After Width: | Height: | Size: 7.7 KiB |
After Width: | Height: | Size: 7.5 KiB |
After Width: | Height: | Size: 7.3 KiB |
After Width: | Height: | Size: 7.2 KiB |
After Width: | Height: | Size: 10 KiB |
After Width: | Height: | Size: 11 KiB |
|
@ -0,0 +1,224 @@
|
|||
%PDF-1.2
|
||||
%äãÏÒ
|
||||
1 0 obj
|
||||
<<
|
||||
/Creator (GKS)
|
||||
/CreationDate (D:20220504164635)
|
||||
/Producer (GKS 5 PDF driver)
|
||||
>>
|
||||
endobj
|
||||
2 0 obj
|
||||
<<
|
||||
/Type /Catalog
|
||||
/Pages 4 0 R
|
||||
/Outlines 3 0 R
|
||||
>>
|
||||
endobj
|
||||
3 0 obj
|
||||
<<
|
||||
/Type /Outlines
|
||||
/Count 0
|
||||
>>
|
||||
endobj
|
||||
4 0 obj
|
||||
<<
|
||||
/Type /Pages
|
||||
/Count 1
|
||||
/Kids [5 0 R]
|
||||
>>
|
||||
endobj
|
||||
7 0 obj
|
||||
<<
|
||||
/Length 12
|
||||
/Filter/ASCIIHexDecode
|
||||
>>
|
||||
stream
|
||||
000000ffffff
|
||||
endstream
|
||||
5 0 obj
|
||||
<<
|
||||
/Type /Page
|
||||
/Parent 4 0 R
|
||||
/Resources << /Font << >>
|
||||
/ExtGState <<
|
||||
/GS255 << /CA 1 /ca 1 >>
|
||||
>>
|
||||
/Pattern <<
|
||||
>>
|
||||
/XObject <<
|
||||
>>
|
||||
>>
|
||||
/MediaBox [0 0 400 200]
|
||||
/Contents 6 0 R
|
||||
>>
|
||||
endobj
|
||||
6 0 obj
|
||||
<<
|
||||
/Length 39798
|
||||
/Filter [/FlateDecode]
|
||||
>>
|
||||
stream
|
||||
xœì½9Ó+Û¶\çç¯8¶^ô<>ÉG…Hz’hÐFK¤"(9úùBó%PíjòkPµ*ãE$pîÝw<C39D>³fFæš5ûw³Ù?³{ü¬f³×»ÿÎw‹Ù×ûú»ÿÿõŸÿÿú<C3BF>ÿe±^ÿóßþ_Ìÿ¹ÿßÿóßð?1‹ý‘«Ú¹ªÿ‘ÑÿvÏ_èú¿àÿø<C3BF>¿Ï|þïv·?ê‘ÿËÝîßÍ÷_oÿûëíóå¿ó?Ì·<C38C>¿ÓÿüÖöüë,׋nÿæÛÍ?‹åí¿ûÏéàßþ3þõ_ÿYÏfø×z¾üÛÿþþÙá_ÿá¿üó¯ÿó?þþõ¿ýó¯ÿŸþ?üçÿõŸÙn9ÛõÄwþË÷Ä÷þ÷7Gü,Â)šóf¢ýÑ+Ï)‰„ÿPßÿÀ#ñõ×+fŒ`#a\<08>o
|
||||
´?JýF¸åj;Çó%<25>ûí>”¸¿¬V»HÎëËaÓ•h~´:ï§C8O‹óþr~&ªoªy¹\××S8¯Ûë×ëýEO|ç¿|O|ïßFø4˜Oª¿ó½ h~]¶ÍgBÕ|Z~1«`#ÁãBèùÖ€ªùL¨šOˇªùL¨šÏ„ªùŸw¤4ÏÓFóIõW¾znÂ$i>3Qó™P5¿”U0‚Œ`„‘!Dn¥ühÛ©ùL¨šÏ„ªùL¨šÏ„ªùÖ<#ADóIõW¾znÂdi>3ts>»™'KóǾ
|
||||
F0‚Œ02„”_‚#?Ú5ŸùnæÉÔ|&TͧÍCÕüW3<57>ªùÖ<#!aÌ'Õ_ù^è¹ #i~èæ|v3<76>¤ùc]#ÁFB¢æ§}/„›ò¿šyDÍgBÕ|&TͧÍCÕ|kžŒ<>ˆ0æ“ê¯|/ÎÑlæÉÖüœ<C3BC>¶ÁÖü®æÿ$ÍÛ*ÁF0ÂÈ25?ü½<C3BC> ùµ<C3B9>¶ŠæÓò[;mS5Ÿ Uó™P5ßšg#DÆ|Rý•ï…„3uX¶và*Í<½ÍÿYš?–U0‚Œ`„‘!ˆš?ëü^ÈÐ|&Úå
Ú<>ªùµfEó™P5ßšg#ô Œù¤ú+ß)gê lgÚÉiæémþ—4è«`#ÁãBèiòÌž§ Uó™P5ÿµÓVÕ|Z>Tͯ5óÑ
|
||||
Ð<#LaÄ'Õßù^ÈÐüY§l{l¾ÉF0Âäzš<=6ßcó<63>`<04>ÍÎ÷‚ ù¡Îùì<C3B9>¶’æÓò‹Y#ÁFBÏ·†Çæ{l¾Œð9„1ŸTå{¡ç&ŒÇæüíl#ÁßFˆÜ"òØ|<7C>Í7‚þaÌ'Õ_ù^è¹ ã±ù;ÁF0·wàzl¾ÇæÁ‡0æ“ê¯|/ôÜ„ñØü‚¿<E2809A><C2BF>`#áÛ‰šŸö½à±ù›o#üBâÙ)iVØ,$›ãø^ˆÞœ¿¿Yl.‰Ù{¦šëí"š½?ÚV3áo„¯7iW
—<?¾
|
||||
Ñ‹µT%ýtþ?"çXè<:Æq,XóŒ0
„;Oø^H¹k³Ø\ۉΰÞ.£™26yÞ^b‰Õfwû—Ïùr}=lÛ‰ê›õÝ÷gἩþöº~'ªo˜·ÿ;\/IiÍ3‚ú¢'8Á‘Fü½<C3BC>¢¿¨š?{Ý„5?|½<>~ù‘4ù?Q¶xþH¶?w½0šca:šg„I#$ÎI?ݦühÛ©ùL¨š›Ÿ¢ùôx¨šÏ„ªùL¨šoÍ3‚z§6Ù‘Fü½<C3BC>¤ùôeÌDͯߜ4Öy½<79>¡ùÄ MþO”í^…?’íÏ]/ŒæX˜ŒæaÚ ‡sF‡fx<rPó™P5Ÿ Uó™P5Ÿ–Uó™P5ßšg#4ºÎs9ÙåH#þ^ÈÒ|&Tͧåc&j~ýzAÐ|fWó–µþØ*ü‘lîza4ÇÂd4ÏÓFH8œ³4ŸÙµÓ6IóiùÏïAóëÍ<‚æ3[Í<©šO›‡ªù´|¨šoÍ3‚žÉóØ<C3B3>:Òˆ¿$Íg¶›y5Ÿ Uóiùñß"ùÕféæôwWat7糯Fs,LFóŒ0m„„ÃYÒü®<C3BC>¶YšÏ„ªù´|ÌDÍo5óäj>ªæÓã¡j¾5ÏSGhžÔ~Ä‘Æû½Ð{g;IóÛÍ<™šÏ„ªùÌþßwà¦Ëö¯Â_Éöç®Ær,LGóŒ0m„øáœô½<C3B4>4YÒ|&TÍgBÕüW3<57>ªùL¨šÏ„ªùSÖ<#L¡ïÔöGñ÷Bï<42>í,ÍŸ½d[Ô|æ»™'Sóû_È´—í^…¿’íÏ]/ŒåX˜ŽæaÚ‰šû^H´Uó™ïfžLÍgBÕüÚN[Eó_Í<ªæ3›Í<Ö<#!Só™Š#<23>ø{¡÷ζ¤ùic0và¦]/$ý¾ ÎÓì—í^…¿’íÏ]/ŒåX˜ŽæaÚ™šß÷½<C3B7>=Oªæ3“Æ懺u^Í<ê<M¨š_Ûi«h~³™Çšg#ˆšÏôØ|<7C>Í_vÊö¯Â_ÉöÇ®Fs,LGóŒ0mQóg<C3B3>ï<05>Í÷Øü<C398>jž¦<>=Á Ž4âï…¤¸›ï±ù©ÇBðèö±0Í3¤礸›ï±ùÕ<#L!pj“iÄßYƒv<6ßcóó¯Fs,LFóŒ0m„„Ã9kÐŽÇæ{lþÀ4ÏÓFè:Ïå¤Çæ{lþ<6C>®ÂÉöç®Fs,LFóŒ0m„„ÃYš§é±ù›?Í3¤Bg¤®ì=ͽo‹÷4ý{!i†är;‹æy{íN4>:®×ûHί›ëqÅDõÍ;‡íáxˆäl{ûñ| ç|~;Ñö$øævRßïg±\^wóz¢ýÑò¸=%æiûX…ßtfä[ë9D¢ý·HÎ.„”¼<þ–¨¾iäi{ŽæÛãã_I(o«ÞŸÿãŽLúõ»ÝD³w$l;ƒ‡¾Þ¤]µöL/L<y5ò3¿ò|îÂs8<08><>0Q<30>î´é„›0Ï\>¤:<3A>©ß‘›0¼©¿EÿÌÍvÅDõM-·Ûu47Û5ž/<2F>ŒüA¨¼I8‘¯Îûu=ÑühµÝ¯O»p^f—Óõ+Q}ÓÎÛåá&5<>œ²æaÚÍsd¬½0ú™ú½ð÷÷‘R¾àk‚‘ùLÌDÍgBÕ|&TÍgBÕ|&TÍgBÕ|Z>B½”„ªùL¨šO›ÇLÔ|&TÍŸuþúŸ¡ùí‘°™šOËOÙT>¤Ÿ¼ºoTLN¶?‡<>t)Só™}…ÕüYâ÷Bô§ÛøÍùÍgBÕ|&TÍg>ÏH‚æ3¡j>ªæOYóŒ0m„¾3e_{aòwdì{áï´MÖ|&TÍÿÊ#TÍgBÕ|&TͧÇCÕ|fâ÷B[ó™P5Ÿ©Þœ¯!HšÏ„ªùL¨šÏDô?”ôë¿ ùÌ®ç;'i~ÊF˜¤_"Sv<53><>œŒl!é>’¨ù}…”¬ù±ï…´Íà3m#šÏÄLÔ|&TÍŸ}<7D>©¡j>ªæ3¡jþ5ÏÓFˆ<46>/›'ÕìïȾoòÌÖ|fóGÛdͧÍCÕ|Z>TÍgBÕ|&TÍOü^è×üv3O¦æ«7ç³›yz5¿ÚÌ#i>ªæ3û›yRæoß<3Qó»šy²4ŸÙÕƒ‘ÕpÚužÔvX¼l!£>ü»sÆlálÍïû^HÖüþgÚ&j>ªæ·›y25Ÿ UóióP5¿ÑÌ“¯ùSÒ<#L!p¦ìlæÉÖü¾ï…_ÜQÕÓ³“r7/©[8[ó™P5Ÿ6UókÍ<Šæ3¡j~ä{!®ùÌÞǵÇ4?÷æ|v3OTóg_Uó«Í<’æ÷7ó$j>ªæW›y$ÍïêÁ<C3AA>öuí:ÏÙ]T®lAØéÚý»sú4‚%TÍŸ5¾²5¿ýLÛLÍïêœÏÒ|fÚõB 5¿ñ£mºæ3¡j>-ªæOAóŒ0m„Îsdà›òÇn…ýÚNÛèÜBÒÜæNÛdÍgBÕ|&TͧåCÕüF3O¾æ÷|/¤k~ïãÚS5ŸùãÍ<ɚτªù3y§mo3O¦æ3¡jþ¬áHÙš_íÁøÖø€öåÿÔeûs‚æ3ë¿;gh~m§¢ù´üð³vàBÕüj缤ùi× ;p¡j>ªæWwÚJš_²æaÚý§Íð7åwo…ýýĶ¤»yYƒv05¿ÒÌ£i>ªæ¿vÚªšOËOé„Iù^È×üÖãÚs5ÿÇ›y²5ŸÙlæIÖüü<C3BC>¶½Í<¢æ3ßÍ<™škxN´ïÙŽdreûcI÷‘’~w4Ÿ UóÃÏ.Ì´Uó™P5?|½<>1hªæךyÍgBÕü5ÏÓFˆŸ6Ãß”ê°ßÛi›¬ù]Ò<M¨šOˇªùÕfIó™P5?¥&in{§mæ ¨šÿSÍ<µ<>¶Ò<M¨šŸ¾Ó6ÒšUó™];m³í4›y²çiö÷lG4Ÿ9yÙþBÒ}¤Ä¸ñgFZó¡j>ªæÓòß;mÅyší<C5A1>¶™ƒvºšy²í@Õ|Z>TÍgBÕü’4Ï“FàyÑcó=6ßcó¿×Ìã±ù›_„l!é>Ræ <1D>Í÷Ø|<7C>Í7”šçH<C3A7>Í÷Ø|<7C>Íךy<6ßcó‹<C3B3>íÏ!$ÝGçizl¾ÇæOUóŒ0m„¾3¥Çæ{l¾Çæ¿$Í÷Ø|<7C>ͯæXdûsI÷‘DÍ÷Øüwzlþ´4ÏÓFˆ<46>/=6ßcó=6?ízÁcó=6¿ÙþBÚ/. ¿;{l~Ò\<5C>ÍŸ†æaÚ<08>3eg3O¶æ÷}/üÚNÛÞÖ|<7C>Í÷Øü_mæñØ|<7C>Í/B¶?‡<>¡ùáß<C3A1>=6?Ió™›_¶æaÚ³ûWø.1#³Âžù×שß)ã‘W‹Ý:œËÓî„ÓòþËÓê|¸\ºï7—ÕåzÝGò|]]OÕDû£ËýT–œˆþ‡Î×u(ñx9\×c$oÿ^¯‡v¢ã¬o¹<6F>äövet¸'ž/<2F>¹ü/‡ó€„¿|ð<>@í£Ý-‘Ü>0_‰öG·L©ˆÅårû‚[àùÌÙmÉzÏ—ó})7‘\^V×e3Qÿèx9ï.§”ÜO³óywÄó¥;WÇýî´'‹í⸠ç|v»v]u'êÍn‡úþËÕéá+§çîÂÊ›Zö‹h¶¯<C2B6>{s³;u'*o–»]4ç<34>ÓY%Ñþè–«h.v›g¢ú¦3·»s(qIø7Zo³®DåÍq½NÌÝzþø£o‰ê›z&Îñaá+øzo^Šr¯WͯDû£øÌ##û’“?‘¤\òœ7d‰ê›Ån»;îSò²¼<C2B2>ì‰kåM=§«yF˜6B®æ3»f…%i>-ÿG:ÏÒ|fÇxä4ͧÇCÕ|&TÍgBÕ|&TͧåCÕ|&TÍgBÕ|f!åÚCÕ|&TÍgBÕ|Z>TÍgBÕ|&TͧÍCÕ|&TÍgBÕüÄ<C3BC>0ýšÏ„ªùL¨šÏ„ªù´|¨šÏ„ªùL¨šÏ/P5¿Þ›'h>ªæ§l„Iɤ'…uj>ªæ3¡j~‰šg„i#¨šß5+,Kó™?òÓ¹¤ùã‘ó4Ÿ UóiùP5Ÿ¿9I¨šO!†ªùL¨šÏ„ªù!„$ͧåCÕ|&TÍgBÕ|Ú<TͧåCÕ|&TÍgBÕ|Z>TÍgBÕüÈF˜¸æÓã¡j>3`ªaÍgBÕ|fß*D5ŸUó™P5Ÿ Uóiù¡gT5Ÿ Uóg<C3B3>›óÙšŸò¤° æÓò¡j>ªæ—¤yF˜4
ªæ·ç?~ÌVÂõBDó+Í<šæ3¡j>ªæÇoÎg7ódj>³Ý “¨ù´|¨š_mæ‘4¿!KóióP5Ÿ–Uó»šy²4Ÿ UóéñP5¿ÖÌ£h>ªæÓ)¡jþìÕÌ#j~ÏF˜tÍgBÕü€©¦i~»™'SóûV!Yó™P5¿ÒÌ£i>ªæ3û;ç#šOËô©†5_mæIzRX’æÓã¡j~¥™g²šg„i#4mªæ¿ç?~ÌVÂõB¢æÓæ¡j~£™'_ó™MÍKÖüþ›óÙÍ<¢æ·;a25Ÿ UóiùP5¿Š i>ªæsáûå‰h~µ™GÒ|&TÍgBÕ|Z>Tͯ5ó(šÏìl#IÑ|Ú<TÍgBÕ|f³™'Yó;L5Oó™P5¿¹
|
||||
ÙšÏìnæIÐüæ*dk~¥™GÓüþ<C3BC>¶‰šÏìèSMÓüÜfž¤'…ei>ªæÓò§¨yF˜6BŸÓw7ó¤ì«ùëÇl%\/dj>ªæÓò¡j~Só²5?}§m¤5ªæ¿;aDÍgBÕüzç¼ ùÝ;m34Ÿ UóûåIÔüÆNÛ|ͯ6óHšÏl5ó¤jþk§ªù´|¨šßÙF’£ùL¨šOˇªù³F3O¶æ÷µ<C3B7>'k~m§¢ù´|¨šßÝÌ“¡ù±–ª„&§þ§`$i~õU’æwô©æi>3ÖÌ“ô¤0Ió™P5Ššg„i#ôiþ¬³™'Cóÿì1[ ×¢æ3kÍ<9š_Ùi«i>-ªæëÍ<
§‡ªù];m³4Ÿ Uó™Í›óÉš»^H´UóÛ¿òdj>ªæÓò¡j~«™'Wó™P5¿±Ó6_ókm$Šæ3¡j~u§¤ù´|¨šÏ„ªùL¨š_Ýi+i~½™GÐü¾–ªdÍïz
|
||||
F–æÏ^Í<¢æ·‡Î‰ûmûšy’í@ÕüF3O¾æOIóŒ0m„˜æ×›yÍgþÚc¶âЗå'ŒG·æCÕ|&Tͯ4?¿™§ç>Tͯ4Ÿ–Uó›7ç³5¿ïz!{ž&TÍÿÊ#j>ªæ3¡jþ«™GÕ|&TÍgBÕ|Zþ{§8Oªæ3¡j~}§ ù´üöNÛÌA;P5Ÿ UóiùP5Ÿ Uó«OÁ<4F>4¿¾ÓVÐ|¦ºÓ¶·™'{ž&TͧåCÕü)hž¦<><C2A6>ªù›?óØ|<7C>Í÷Ø|<7C>Í÷Ø|<7C>Í÷Ø|<7C>Í“æaÚ¹šÏôØ|<7C>ͯï´õØ|<7C>Í÷Ø|<7C>Í÷Ø|<7C>Í÷Øü<C398>iž¦<> j~×Æy<C386>Í÷Ø|<7C>Í÷Ø|<7C>Í÷Ø|<7C>Í÷Øüدÿ›oS5Â_ Ì:wÚzl~–æ{lþ;=6ßcó=6ßcó=6ÿ<36>›ï±ù6U#|a¾ßì—w¢ú¦ž<C2A6>ï¥y,±Üm··ƒ4œI{ÿӮΩg¤ívÎùis<Îæ'<_ºó¼¸,oRL<_nÿ·»^"¹¹®×f¢ñÑ
|
||||
©€cÁƒCØ<6|Ó“‡Ð;ñxYݾä/IÙ:FÑqØ.nyŽäüñHTß4r–’÷‰ç¤+õËãò»•¨}”r×àTO´?J¼<15>4˜÷™íÿ<C3AD>JâëÍîq9ÌóùvYºm&Ú]·§Íùx¹†s»?nÎûíÏ—î\,·ËÛ5e0;‡õeÚØü”õ<>9y31’Zm#Y¿›ÿjßlWõDû£{Þeû0•@ö7yV²ü“ª&‡ÐçôP5Ÿ6Uó™©<E284A2>Ùji~_{a²æ3¡j>-ªæ3¡j¾<6A>#”€P±y¨šOˇªùL¨šÏ„ªù×/Ë<>÷æõh>ªæ7e;[óÃy4Ÿ–Uó™P5Ÿ UóiðIO
|
||||
K›Ÿ©ùÍ9ٚόýtžx½ h>ªæw5yfiþØOªF˜BŸæ3¡j>ªæ‡ÎHIšÏ„ªùL¨šO-€ªù´üRÇ‚‡ÐáôèøIšO<C5A1>‡ªùÕÃYÒ|&®¢æÇ{ó"šOË
æ
j>ªæ3»ž²<C5BE>¤ù´y¨šÏ„ªùL¨šŸô¤°´±ù¢æw=ÔænÎçfä÷… æ3¡j~ÏĶòOªF˜BLó+Í<šæ3¡j~ì1[QͧåCÕüZ3<5A>¢ù¡_ÿ“4…TÀ±`„Á!ºuªÍ<YšÏ„ªù´|¨šÏ„ªùýç5ŸÙ5˜7Ióiù¡›óÙÍ<YšÏ„ªù´|¨š_kæQ4?ø¤°¤A;½Ï´MÔüªlÿÀÍyAócwó¢šOˇªù<C2AA>‰måŸT<C5B8>09„TͧåCÕüz3Ï7nBt]þ'i>mªæÓò¡j~ׯÿYš?àB*àX0ÂàbMùÇwc[¶æ3›Í<ɚτªù4x¨šßÞ8Ÿ©ù]ƒy³4?ts>»™GÒ|&TͧÍCÕ|Z>TÍïlæÉÑüða2fb|ïæ|v3Oä÷…,ͧÍCÕ|Z~ñ'U#L!Wó_;mUͧåÇ÷þ'>f+[ó™P5¿±Ó6_ó«¿þKš?ÀB*àX0Âà4¿¾ÓVÐüf3O¶æ3¡j~m§¢ù×W3<57>¨ùÌw3O¦æçì´
¶æCÕüJ3<4A>¦ùL¨šÿÚi«j~™GÑ|&TÍÿ¥<C3BF>¶šÏäݼlÍgBÕüâOªF˜‚ªùL¨šß}FÊ?;<3B>SäÑÛUó™P5Ÿ–»!9†B*àX0Âà24¿oJU²æÓòS6λu j>³³™'gÐTÍgÆd;ºWiæéhþ4Ÿ–Uó™P5Ÿ Uóiùñ±ù‘_Å¡j>ó—›yvàBÕ|&TÍ/ö¤j„©!tì´ÕæiBÕüþöÂÌA;±yDwàBÕ|&TÍ<54>Ý<EFBFBD>A!p,ax‚æ3¡j~ÊÆù¤¸Õ<C2B8>¶Ò<M¨šÏ„ªù1ÙN´“ÓÌÓù?ÓlæÉž§ Uó_Í<ªæ3¡j>³l~bókÊX›_ÜiÛ{½<> ù<C2A0><C3B9>¶ùš_iæÑ4¿¸“ª¦†à±ù›ïcÁEðØ|<7C>Í÷Ø|<7C>Í÷Ø|/AGðØ|<7C>Í÷±`„<>"xl¾Çæ{l¾Çæ{l¾¿Œ #xl¾ÇæûX0Â@<6ßcó=6ßcó=6ßßFÐ<6ßcó},a ›ï±ù›ï±ù›ïï#è›ï±ù>Œ0P<04>Í÷Ø|<7C>Í÷Ø|<7C>Í÷÷‚t„Óù2»ÎSr{>îÎ&ªoÞ9_®¯‡m8oÇåqû¸ü¿½rµ<72>‡•7Ëí,š§ÇZI4?Z^v—Ó*œ—ÿ½'Œ¨~”THF0‚Œ0„´ùL÷Lù•ç²¿]·ÂyØÞ.¢VÏDõM5S¾ï‰C|ˆÅ3#}ªÏLQîY=Ñþ(/±9>Þ„3ðGàëMúUÃz»¨&Ú-VÛÕqήkçVŽéX(àp6ÂwTÍgBÕ|&TÍgBÕ|&TÍ<54>íåI¾Q_B!ÁF0¢[¶ß™ò+O®æ3¡jþ!¸¯-AóëCçÍoÊv~~Ù<Tͧå§ÿ¾Ðpz&TÍ¿öo„ѱPÀál„ï"Ðã¡j>ªæ3¡j>=ªæ3¡j~ß^žì~œ’Œ`#!ºeû<65>)GV5Ÿ Uóg<C3B3>ûÚ24Ÿ–ïœ<C3AF>Þœ5Ÿ UóÓ›ÿ{4Ÿ Uó;6ÂŒèX(àp6Âwš6UóiùP5ŸùnæÉÔ|&TͧÍCÕüæ^©í~¤…d#Áƒ@ˆnÙ~g‚æÓæ¡jþ«™GÕüþ!Õ‰šÏìïœO4øôfžžøP5?ŽÑ|Z>Tͯ4óŒèX(àp6Âwúœ¾£™'Mó+í…šæÏ^Í<¢æ3›Í<ɚτªù´ü¾gÚ&?øjD…d#Áƒ@ˆnÙ~”¡ùL¨šÏ/Á´[a<>n<EFBFBD>w3O¦æ÷ï´Íìʉ7óDZóS¯2&y&j>ªæÓòÇq,p8á»±¦|¨šÏ츃‘¦ù];m³4ÖhæÉÖ|&TÍï{¦möómGPHF0‚Œ0„¤¸OAó™P5?íVXÂ\¨šßî„‘wÚŠšŸz½<7A>1É3Só™P5ÇB‡³¾‹<C2BE>ºªæwÜÁÈÓüêN[IóiùP5¿ÖÌ£h>³kô_)…d#Áƒ@È´}8r¬5ªæ‡o…eÚ<>ªùL}§íW¶oÎgÚé»^Șä)j>ªæûX(àp6Âwrí@Õüú¬0Aó™]Í<Yƒv j>-ªæ3¡jþ€ÉF0‚<06> ÍÓì}8rê\¨šß}+L˜§ Uóõ<C3B3>¶<EFBFBD>„ªù}×Ùó4¡j>-¿ýȹÌA;Ã:
|
||||
8œ<EFBFBD>ð]<04>Í÷Ø|F0‚ÊAðØ|<7C>Í÷Øübg#|Ácó=6ßÇ‚Œ`„b<6ßcó=6¿˜ÃÙßEðØ|<7C>Í÷±`#¡<04>Í÷Ø|<7C>Í/æp6Âw<6ßcó},ÁF(Ácó=6ßcó‹9œ<39>ðm„õãtÖJ´?ºÜ_Oá<.΋Ëîžx¾tå|µ™"y?œÛÃíðLÉ΃ºÙ›wÜù¸ÿ+9½Õ7ÌËòFºOI<4F>Í7‚Œ`„±#ô5óT󼼬®Ëpί›Ëñ–x¾„òö<C3B2>´ßîú¯7išn&Ú%÷Æ|³™'ûz!ù÷…ÎVVþE¢ûßíf}8œ—ἩÎáéH‡‡ug‘Ç‚F‰Ð£ùL¨šÏ„ªùôx¨š?kôæek>ªæ{l¾Œ`#Œ¡¯™GÕ|&TÍ›ª¨ùL¨š¯7ó4ò¯25Ÿ Uó™P5ŒÇ‚F‰Ñ|&TÍgBÕ|&TͧåCÕ|&¶¢æ3=6ßF0‚Æ…iì<69>ªùL¨šŸ:6?ª¿ý<C2BF>0‰í=¿ps>ûzAÔüÆO$ùšÏ„ªùc:Œ0J„DͧåCÕüJ3<4A>¦ùÌ®fž$ͧÍCÕüí—åCÕ|<7C>Í7‚Œ`„q!$ôïß=ªæw4óäi~ll~\¿<>e¬Í_ÝœÏÍ@?R’æÓò¡j>ªæ<C2AA>áX0Â(25Ÿ UóiùP5¿«™'Kó™P5ŸUó=6ßF0‚Æ<><C386>¨ù´y¨š_iæÑ4ŸÙ›Ÿ¬ùÌÐX›¿º9/j~W?R–æÓæ¡j>-¿ÈcÁ£D5Ÿ Uó™P5¿ÚÌ#i>ªæ3¡j~Q…d#Á"dj>ªæÓòõÒıùÙšÏìjæùôNÛDÍgBÕ|&TÍ/òX0ÂÚ;mó4Ÿ Uó™P5ÖØi›ù¯fUó™P5¿„B2‚Œ`„DÍgBÕ|fÏL˜¸æÓò»ïlgÚ<>ªùÌ?hæ‰ |