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

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

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

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

1. meshの描画
　粒子による描画をmeshによる描画にするには、粒子の集まりの表面（のようなもの）に
　meshを張ります。具体的には、粒子の密度がゼロとなる境界にmeshを張っていきます。
　（「修士論文 水の実時間アニメーション 天田崇」
　https://library.naist.jp/mylimedio/dllimedio/showpdf2.cgi/DLPDFR003344_P1-56
　の方法を参考にしました。）

・meshの生成
　①distance field 値を計算する。
　　陰関数には、粒子の密度を用いる。
　②distance field 値から、Marching cubes法を用いて、meshのvertexとfaceを
　　生成する。
　GeometryTypesとMeshingライブラリを使用します。

using Meshing
using GeometryTypes

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

    DT = 0.004
    positions = Simulation(DT)
    resolution  = 0.75
    sdf = SignedDistanceField(HyperRectangle(Vec(-2.,-2,-2),Vec(14.,30,20)), resolution) do 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
    
    m = SimpleMesh(sdf, MarchingCubes())
end

# 描画
N = 200

for i = 1:N
    m = create_mesh(i)
    m.vertices .*= 0.1
    setobject!(vis, m, material)
    # transform 関連
end

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

povファイル

#include "colors.inc"
#include "glass.inc"

global_settings{max_trace_level 20}
camera{ location <-2.5, 1.0, -2.8> 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{
   vertex_vectors { 6296,
      <0.1,0.1,-0.058611978>
      <0.1,0.07704518,-0.05>
      <0.07704518,0.1,-0.05>
          ・・・（途中省略）
      <0.8040384,0.925,1.675>
      <0.775,0.92816883,1.675>
      <0.775,0.925,1.6781456>
   }
   face_indices { 3143,
      <3,2,1>
      <6,5,4>
      <6,4,7>
        ・・・（途中省略）
      <6287,6286,6285>
      <6290,6289,6288>
     <6293,6292,6291>
   }

   rotate x*(-90)
   texture { T_Glass3 } interior{ I_Glass ior 1.33 }
}

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

function writeData(i, mesh)
    filename = string("out", lpad(i,4,"0"), ".dat")
    out = open(filename, "w")
    # vertices
    println(out, "vertex_vectors { ", length(mesh.vertices), ",")
    map(mesh.vertices) do vertex
        println(out, "<", Float32(vertex[1]), ",", Float32(vertex[2]), ",", Float32(vertex[3]), ">")
    end
    println(out, "}")
    # face indices
    println(out, "face_indices { ", length(mesh.faces), ",")
    map(mesh.faces) do face
        println(out, "<", face[1], ",", face[2], ",", face[3], ">")
    end
    println(out, "}")
    close(out)
end