Zenn
Closed5

TopOpt.jlを使って何か遊びつつドキュメントを残す

Julia
topopt
トポロジー最適化
MarieΔMarieΔ

TopOpt.jlについて

TopOpt.jlはJuliaのトポロジー最適化用ライブラリです。

以下はgithubに書いてあった情報をまとめています。

下記のように様々なペナルティの取り方に対応しています。

  • Solid isotropic material with penalization (SIMP)
  • Rational approximation of material properties (RAMP)
  • Continuation SIMP/RAMP
  • Bi-directional evolutionary structural optimization (BESO) with soft-kill
  • Topology optimization of binary structures (TOBS)

また最適化の手法に関しても以下のようなものが利用可能です。

  • Method of moving asymptotes
  • All the algorithms in NLopt
  • Ipopt
  • First order augmented Lagrangian algorithm
  • Nonlinear semidefinite programming for buckling constrained optimization
  • Basic surrogate assisted optimization and Bayesian optimization
  • Integer nonlinear programming (design variables guaranteed to be integer)
  • Sequential integer linear programming in the topology optimization for binary structures (TOBS) algorithm
    (Nonconvex.jlというものを利用しているので、それに準拠するものなら全て使えそうです。)
MarieΔMarieΔ

TopOpt.jl のinstallの仕方について

using Pkg
pkg"add TopOpt"

描画に関しては、Makie.jlを使っています。このライブラリは別途入れてください。

環境について

OS X(intel), julia 1.6.4(LTS)では動作を確認しています。
逆にDockerでこのバージョンで環境構築しようとした場合うまくいかなかったです。(1.6-1.8まで確認済み)

MarieΔMarieΔ

点荷重片持ち梁問題

ここからはTutorialに移っていきます。

ライブラリのインポート

using TopOpt
using Makie
using GLMakie
using TopOpt.TopOptProblems.Visualization: visualize

(結果を視覚化して出力するので、TopOpt以外のライブラリも入れています。単純に計算のみではTopOptだけ入れればいいです)

パラメータの定義

E = 1.0 # ヤング率
v = 0.3 # ポアソン比
f = 1.0; # 点荷重

nels = (30, 20, 10) # 設計空間
problem = PointLoadCantilever(Val{:Linear}, nels, (1.0, 1.0, 1.0), E, v, f); # 点荷重片持ち梁問題呼び出し

ちなみに、PointLoadCantilierの詳細は以下のようになっています。

///**********************************
///*                                *
///*                                * |
///*                                * |
///********************************** v


@params struct PointLoadCantilever{dim, T, N, M} <: StiffnessTopOptProblem{dim, T}
    rect_grid::RectilinearGrid{dim, T, N, M}
    E::T
    ν::T
    ch::ConstraintHandler{<:DofHandler{dim, <:Cell{dim,N,M}, T}, T}
    force::T
    force_dof::Integer
    black::AbstractVector
    white::AbstractVector
    varind::AbstractVector{Int}
    metadata::Metadata
end

束縛パラメータの設定

V = 0.2 # 体積の束縛
xmin = 1e-6 # 最小密度
rmin = 2.0; # 密度に関するフィルタ(おそらく平滑化)の半径

ペナルティとソルバーの定義

penalty = TopOpt.PowerPenalty(3.0) ### pρKのpの設定
solver = FEASolver(Direct, problem; xmin=xmin, penalty=penalty)

コンプライアンスマトリクスの設定

comp = TopOpt.Compliance(solver) # コンプライアンスマトリクスの計算
filter = DensityFilter(solver; rmin=rmin) # フィルタの設定
obj = x -> comp(filter(PseudoDensities(x)))

体積制約の設定

volfrac = TopOpt.Volume(solver) # 体積の計算
constr = x -> volfrac(filter(PseudoDensities(x))) - V # 体積制約との誤差計算

最適化処理の設定

x0 = fill(V, length(solver.vars)) # 配列生成
model = Model(obj) 
addvar!(model, zeros(length(x0)), ones(length(x0)))
add_ineq_constraint!(model, constr)
alg = MMA87() #最適化アルゴリズムの定義
convcriteria = Nonconvex.KKTCriteria() # Karush-Kuhn-Tucker 残差
options = MMAOptions(;
    maxiter=3000, tol=Nonconvex.Tolerance(; x=1e-3, f=1e-3, kkt=0.001), convcriteria
) # イテレーション, 終了条件, 誤差函数の設定
r = optimize(model, alg, x0; options) # 最適化器の設定

@show obj(r.minimizer)

実際にプログラムが走るとこのようなログが出力されます。

[ Info:   iter       obj      Δobj  violation  kkt_residual  
[ Info:      0   5.0e+03       Inf   1.1e-14   8.0e+03
[ Info:      1   2.3e+03   2.7e+03   0.0e+00   2.1e+03
[ Info:      2   1.1e+03   1.2e+03   0.0e+00   1.9e+03
[ Info:      3   3.8e+02   7.6e+02   0.0e+00   5.6e+02
[ Info:      4   1.3e+02   2.5e+02   0.0e+00   8.2e+01
[ Info:      5   7.6e+01   5.4e+01   0.0e+00   9.7e+00
[ Info:      6   6.0e+01   1.6e+01   0.0e+00   3.1e+00
.
.
.
obj(r.minimizer) = 36.48359966639993

値はプログラムを走らせる事で当然変わるので、表記のみ参考にしてください。

グラフの保存

fig = visualize(
    problem; topology = r.minimizer,
    default_exagg_scale = 0.07, scale_range = 10.0,
    vector_linewidth = 1, vector_arrowsize = 0.3,
)

save("opt.png", fig)

実際に実行するとこのようなグラフが得られます。無事構造を求めることができました。
(みづらいですが、x, yが底面になっており、z方向に荷重がかかっています。)

MarieΔMarieΔ

系の切り替えについて

デフォルトで定義されている系は以下です。

TopOpt.TopOptProblems.PointLoadCantilever
TopOpt.TopOptProblems.HalfMBB
TopOpt.TopOptProblems.LBeam
TopOpt.TopOptProblems.TieBeam

これ以上の場合、別途系を定義する必要があります。

系の定義には、Typeとmethodどちらの定義も必要となります。
具体的には以下のように作成します。

using Ferrite
using TopOpt
using TopOpt.TopOptProblems: RectilinearGrid, Metadata
using TopOpt.TopOptProblems:
    left,
    right,
    bottom,
    middley,
    middlez,
    nnodespercell,
    nfacespercell,
    find_black_and_white,
    find_varind
using TopOpt.Utilities: @params

@params struct Hoge{dim,T,N,M} <: StiffnessTopOptProblem{dim,T}
    rect_grid::RectilinearGrid{dim,T,N,M}
    E::T
    ν::T
    ch::ConstraintHandler{<:DofHandler{dim,<:Cell{dim,N,M},T},T}
    load_dict::Dict{Int,Vector{T}}
    black::AbstractVector
    white::AbstractVector
    varind::AbstractVector{Int}
    metadata::Metadata
end

function NewPointLoadCantilever(
    ::Type{Val{CellType}},
    nels::NTuple{dim,Int},
    sizes::NTuple{dim},
    E=1.0,
    ν=0.3,
    force=1.0,
) where {dim,CellType}
.
.
.
    return Hoge(
        rect_grid, E, ν, ch, load_dict, black, white, varind, metadata
    )
end
TopOptProblems.nnodespercell(p::Hoge) = nnodespercell(p.rect_grid)

function TopOptProblems.getcloaddict(p::Hoge{dim,T}) where {dim,T}
    return p.load_dict
end

また、Methodの書き方ですが、PointLoadCantileverの例を以下に載せておきます。

function PointLoadCantilever(
    ::Type{Val{CellType}},
    nels::NTuple{dim,Int},
    sizes::NTuple{dim},
    E=1.0,
    ν=0.3,
    force=1.0,
) where {dim,CellType}
    iseven(nels[2]) && (length(nels) < 3 || iseven(nels[3])) ||
        throw("Grid does not have an even number of elements along the y and/or z axes.")

    _T = promote_type(eltype(sizes), typeof(E), typeof(ν), typeof(force))
    if _T <: Integer
        T = Float64
    else
        T = _T
    end
    if CellType === :Linear || dim === 3
        rect_grid = RectilinearGrid(Val{:Linear}, nels, T.(sizes))
    else
        rect_grid = RectilinearGrid(Val{:Quadratic}, nels, T.(sizes))
    end

    if haskey(rect_grid.grid.facesets, "fixed_all")
        pop!(rect_grid.grid.facesets, "fixed_all")
    end
    #addfaceset!(rect_grid.grid, "fixed_all", x -> left(rect_grid, x));
    addnodeset!(rect_grid.grid, "fixed_all", x -> left(rect_grid, x))

    if haskey(rect_grid.grid.nodesets, "down_force")
        pop!(rect_grid.grid.nodesets, "down_force")
    end
    addnodeset!(
        rect_grid.grid, "down_force", x -> right(rect_grid, x) && middley(rect_grid, x)
    )

    # Create displacement field u
    dh = DofHandler(rect_grid.grid)
    if CellType === :Linear || dim === 3
        push!(dh, :u, dim) # Add a displacement field
    else
        ip = Lagrange{2,RefCube,2}()
        push!(dh, :u, dim, ip) # Add a displacement field        
    end
    close!(dh)

    ch = ConstraintHandler(dh)

    #dbc = Dirichlet(:u, getfaceset(rect_grid.grid, "fixed_all"), (x,t) -> zeros(T, dim), collect(1:dim))
    dbc = Dirichlet(
        :u, getnodeset(rect_grid.grid, "fixed_all"), (x, t) -> zeros(T, dim), collect(1:dim)
    )
    add!(ch, dbc)
    close!(ch)
    t = T(0)
    update!(ch, t)

    metadata = Metadata(dh)

    fnode = Tuple(getnodeset(rect_grid.grid, "down_force"))[1]
    node_dofs = metadata.node_dofs
    force_dof = node_dofs[2, fnode]

    N = nnodespercell(rect_grid)
    M = nfacespercell(rect_grid)

    black, white = find_black_and_white(dh)
    varind = find_varind(black, white)

    return PointLoadCantilever(
        rect_grid, E, ν, ch, force, force_dof, black, white, varind, metadata
    )
end
このスクラップは2023/10/18にクローズされました