Julia MeshCat SPH 流体シミュレーション (Mesh)

update 2022/03/03
Julia MeshCat SPH 流体シミュレーション (Mesh) の記事のupdate です。
GeometryTypes ライブラリーが GeometryBasics に変更されました。
この変更により、MeshCat でメッシュの表示ができなくなっていました。
この変更に合わせて、Meshing ライブラリーも update され、新たな書式で
表示できるようになっています。
Meshデータの生成は、直接 pov ファイルを生成するように変更しました。
（変更後のプログラムを掲載しています。）

Julia MeshCat ライブラリを用いた、水柱崩壊のシミュレーションです。
流体のシミュレーションには、SPH (Smoothed Particle Hydrodynamics) 法を使用しています。
前回は、PointCloudを用いて、粒子による描画を行いました。
今回は、Marching cubes法を用いて、meshによる描画を行います。
水柱崩壊のシミュレーション（meshによる描画）

水柱崩壊のシミュレーション（povrayによる描画）

1. meshの描画
　粒子による描画をmeshによる描画にするには、粒子の集まりの表面にmeshを
張ります。具体的には、粒子の密度と静止状態の密度の差(density - rest_density)
がゼロとなる境界にmeshを張っていきます。
これは、「修士論文 水の実時間アニメーション 天田崇」
　(https://library.naist.jp/mylimedio/dllimedio/showpdf2.cgi/DLPDFR003344_P1-56)の方法を参考にしました。
・meshの生成
　①distance field 値を計算する。
　　陰関数には、粒子の密度を用いる。
　②distance field 値から、Marching cubes法を用いて、meshのvertexとfaceを生成する。

## MeshCat-sph-mc.jl

using MeshCat
vis = Visualizer()
#open(vis)

using Blink
AtomShell.isinstalled() || AtomShell.install()
open(vis, Window())

using Meshing
using ColorTypes
using CoordinateTransformations
using Rotations
using LinearAlgebra
using LinearAlgebra: norm
using GeometryBasics
using GeometryBasics: Mesh

mutable struct Particle
    position::Point3
    velocity::Point3
    density
    pressure
    force::Point3
end

# weight function
w_rho(r, h) = r < h ? (h^2 - r^2)^3 : 0.0
w_pressure(r,h) = r < h ? (h - r)^2 : 0.0
w_viscosity(r,h) = r < h ? (h- r) : 0.0

# Gravity
accel_g = [0.0, 0.0, -9.8]

# constants
mass = 0.00020543
pressure_stiffness = 1.0
rest_density = 600.0
h = 0.01
l0 = (mass / rest_density)^(1.0 / 3.0)
simscale = 0.004
d = l0 / simscale * 0.87
visc = 0.2
limit = 200.0
radius = 0.004
epsiron = 0.00001
extstiff = 10000.0
extdamp = 256.0

w_rho_coef = 315.0 / (64.0 * pi * h^9)
w_pressure_coef = 45.0 / (pi * h^6)
w_visc_coef = w_pressure_coef

init_min = [ 0.0,  0.0, 0.0]
init_max = [10.0, 20.0, 10.0] 
pos_min = [0.0, 0.0, 0.0]
pos_max = [10.0, 20.0, 20.0]

nx = 5; ny = 5; nz = 8

particle_list = []

for iz = 0:nz+1, iy = 0:ny+1, ix = 0:nx+1
    x = init_min[1] + d*ix; y = init_min[2] + d*iy; z = init_min[3] + d*iz

    position = Point3(x, y, z)
    velocity = Point3(0, 0, 0)
    density = 0.0
    pressure = 0.0
    force = Point3(0.0, 0.0, 0.0)
    push!(particle_list,
       Particle(position, velocity, density, pressure, force))
end
particles = length(particle_list)
   println("number of particles: ", particles)
## number of particles: (nx+2)*(ny+2)*(nz+2)

function Simulation(DT)
   # neghiber_list
    neighbor_list = []
    for i = 1:particles
        position_i = particle_list[i].position
        neighbors = []
        for j = 1:particles
            j == i && continue
            position_j = particle_list[j].position
            distance = norm(position_j - position_i) * simscale
            if distance < h
                push!(neighbors, j)
            end
        end
        push!(neighbor_list, neighbors)
    end

    # density(rho)
    for i = 1:particles
        position_i = particle_list[i].position

        sum = 0.0
        map(neighbor_list[i]) do j
            position_j = particle_list[j].position
            distance = norm(position_j - position_i) * simscale
            sum += w_rho(distance, h)
        end
        # density
        particle_list[i].density = sum * mass * w_rho_coef

        # pressure
        particle_list[i].pressure =
            pressure_stiffness * (particle_list[i].density - rest_density)
    end
    
    # force: pressure gradient, viscosity and adj(force at wall)
    for i = 1:particles
        pressure_gradient = [0.0, 0.0, 0.0]
        visc_term = [0.0, 0.0, 0.0]
        position_i = particle_list[i].position
        
        map(neighbor_list[i]) do j
            position_j = particle_list[j].position
            distance = norm(position_j - position_i) * simscale
            
            diffji =  (position_j - position_i) * simscale
            pweight = w_pressure_coef * w_pressure(distance, h) * diffji / distance
            
            pji = particle_list[j].pressure + particle_list[i].pressure

            rho_i = particle_list[i].density
            rho_j = particle_list[j].density
            pgradient = -0.5 * pji / (rho_i * rho_j) * pweight
            pressure_gradient .+= pgradient

            vweight = w_visc_coef * w_viscosity(distance, h)
            vji = particle_list[j].velocity - particle_list[i].velocity
            vterm = visc  * vji / (rho_i * rho_j) * vweight 
            
            visc_term .+= vterm
        end
        particle_list[i].force = pressure_gradient + visc_term
    end
    
    # adj force
    accel_adj = Point3(0.0, 0.0, 0.0)
    diff_min = Point3(0.0, 0.0, 0.0)
    diff_max = Point3(0.0, 0.0, 0.0)
    map(particle_list) do particle
        accel = particle.force * mass

        speed = norm( accel )^2;
        if speed > limit^2
            accel *= limit / sqrt(speed)
        end
        
        diff_min = 2.0 * radius .- ( particle.position - pos_min ) * simscale
        diff_max = 2.0 * radius .- ( pos_max - particle.position ) * simscale
        
        # Z
        if diff_min[3] > epsiron
            normal = Point3( 0.0, 0.0, 1.0 )
            adj = extstiff * diff_min[3] - extdamp * dot( normal, particle.velocity )
            accel += adj * normal
        end
        if diff_max[3] > epsiron
            normal = Point3( 0.0, 0.0, -1.0 )
            adj = extstiff * diff_max[3] - extdamp * dot( normal, particle.velocity )
            accel += adj * normal
        end

        # X
        if diff_min[1] > epsiron
            normal = Point3( 1.0, 0.0, 0.0 )
            adj = extstiff * diff_min[1] - extdamp * dot( normal, particle.velocity )
            accel += adj * normal
        end
        if diff_max[1] > epsiron
            normal = Point3( -1.0, 0.0, 0.0 )
            adj = extstiff * diff_max[1] - extdamp * dot( normal, particle.velocity )
            accel += adj * normal
        end
        
        # Y
        if diff_min[2] > epsiron
            normal = Point3( 0.0, 1.0, 0.0 )
            adj = extstiff * diff_min[2] - extdamp * dot( normal, particle.velocity )
            accel += adj * normal
        end
        if diff_max[2] > epsiron
            normal = Point3( 0.0, -1.0, 0.0 )
            adj = extstiff * diff_max[2] - extdamp * dot( normal, particle.velocity )
            accel += adj * normal
        end
        
        accel += accel_g

        veloc = particle.velocity
        veloc = veloc .+ accel * DT

        posit = particle.position
        posit = posit .+ veloc * DT / simscale
        particle.position = posit
        particle.velocity = veloc
    end
    positions = map(particle_list) do particle
        particle.position
    end
    positions
end

function create_mesh(nloop)
    println("nloop ", nloop)

    DT = 0.004
    positions = Simulation(DT)

    function sdfunc(v)
        sum = 0.0
        map(positions) do position
            distance = norm(v - position) * simscale
            sum += w_rho(distance, h)
        end
        # density
        density = sum * mass * w_rho_coef
        density - rest_density
    end

    m1 = Mesh(sdfunc, Rect(Vec(-2.,-2,-2),Vec(14.,24,20)),
            MarchingCubes(), samples=(20,20,20))
end

### Mesh settings
color = RGBA{Float32}(0.0, 0.84, 1.0, 0.7)
material = MeshPhongMaterial(color = color)

N = 200

for i = 1:N    
    m = create_mesh(i)

    m_vertices = decompose(Point{3, Float32}, m)
    m_vertices .*= 0.1

    ### display Mesh
    setobject!(vis, m, material)

    trans = Translation(0., -1, 0)
    rot = LinearMap(RotYZ(0.2, -0.2))
    transform = compose(rot, trans)
    settransform!(vis, transform)
    
    #sleep(0.01)
    #Base.prompt("continue?") # for capture
end

2. povrayによる描画
　povファイルは、「粒子法入門　越塚他著」にあるpovファイルを参考にして
います。

meshデータの生成
meshデータ（vertex_vectorsとface_indices）は以下のようにして生成します。

### create pov file
function writeData(i, mesh)
    string1 = ["#include \"colors.inc\"", "#include \"glass.inc\"", "",
    "global_settings{max_trace_level 20}",
    "camera{ location <-2.5, 0.7, -2.5> look_at <0.5, 0.5, -1.5>}",
    "light_source{ <0.5, 1, 0.25> color rgb<1, 1, 1>}",
    "sky_sphere{",
    "   pigment{", "      gradient y",
    "      color_map{",
    "         [0.0 color rgb<1.0,1.0,1.0>*0.5]",
    "         [0.2 color rgb<1.0,1.0,1.0>*0.2]", "      }", "   }", "}", "",
    "plane { y, 0.0",
    "   texture {",
    "      pigment{checker color White color Gray70 scale 0.2}", "   }",
    "   finish { ambient 0.5}", "}", "", "mesh2{"]

    string2 = ["   rotate x*(-90)",
        "   texture { T_Glass3 } interior{ I_Glass ior 1.33 }", "}"]

    filename = string("out", lpad(i,4,"0"), ".pov")
    out = open(filename, "w")
    # strin1 (pre)
    foreach(line -> println(out, line), string1)

    # vertices
    vertices = decompose(Point{3, Float32}, mesh)
    println(out, "vertex_vectors { ", length(vertices), ",")
    foreach(vertex -> println(out, "<", vertex[1], ",", vertex[2], ",", vertex[3], ">"),
        vertices)
    println(out, "}")
    
    # face indices
    ### notice
    ### for i = 1:length(faces)
    ### => Parse Error: Mesh face index out of range.
    ### => "for i = 1:(length(faces) - 1)"
    faces = decompose(TriangleFace{Int}, mesh)
    println(out, "face_indices { ", length(faces) - 1, ",")
    for i = 1:(length(faces) - 1)
        println(out, "<", faces[i][1], ",", faces[i][2], ",", faces[i][3], ">")
    end
    println(out, "}")

    # strin2 (post)
    foreach(line -> println(out, line), string2)

    close(out)
end

このwrireData関数を上のプログラムに追加し、
m_vertices .*= 0.1 の下で、writeData(i, m) を実行します。