SPHによる巻き波のシミュレーション第２回

前回のエントリでは、可動壁によって波を起こすことで巻き波のシミュレーションをしようとしたが、うまく巻き波が起きなかった。

そこで、巻き波を起こすにはどうしたらいいのかいろいろと調べてみたところ、RadovitzkyとOritzの"Lagrangian finite element analytics of newtonian fluid flows."という論文を見つけた。

この論文は、有限要素法によってラグランジュ的に自由表面をシミュレーションする手法を提示している。その適用例の１つとして砕波を扱っており、そこに孤立波の理論解が示されていた。これを適用すれば波の岸に打ち寄せる波を再現することができそう。

孤立波の初期条件は、水面の高さをη、x方向の初速をu、y方向の初速をvとすると、



のようになる（x方向のオフセットを省略）。ここで、gは重力加速度で9.8[m/s^2]、Hとdはそれぞれ基準となる水面の高さと波の高さである。

また、tanhとsechは双曲線関数で、その定義は

である。

この初期条件を適用して、またシミュレーションしてみると、



お！巻き波が起きた！

シミュレーション結果を見ると、左から進んできた孤立波が、斜面にさしかかったところでせり上がって、巻き波になってくずれる様子が見てとれる。また、巻き波が打ち寄せたあとにさーっと引いていく様子も、実際の砂浜とよく似ている。

いい感じに巻き波が起きるようになったので、次は、このエントリのように綺麗にレンダリングしてみます。

■関連するエントリ
SPHによる巻き波のシミュレーション第１回
SPHによる巻き波のシミュレーション第３回
-

SPHによる巻き波のシミュレーション第１回

波の一種に巻き波というものがある。浅海域に進入した波が対称性を失って、前面が切り立った壁のようになり、頂部の水が前に飛び出す形で巻くように砕けるものをいう。

粒子法では、しぶきがあがるような激しい変形を表現することができる。そこで、巻き波をシミュレーションしたいと思った。



左側の境界条件を周期的に動かすことで波を起こしている。また、真ん中から右は斜面になっている。

ところが、結果をみてみるとうまくいかない…

波が減衰して斜面まで届かないのと、斜面の底部分で不自然に左向きの流れが起きているのが気になる。

波が減衰してしまう点については、容器の形や可動壁の動きを調整しようと思う。

斜面の底部分で不自然な流れが起きている点については、理由がよくわからない…境界条件にはペナルティ法を使っているが、そのせいで斜面方向に変に力がかかっているんだろうか？

もうすこし調整が必要そうだ。

■関連するエントリ
SPHによる巻き波のシミュレーション第２回
SPHによる巻き波のシミュレーション第３回
-

【論文】p = k ( ρ - ρ0 )でのρ0の意味

Mullerらの"Particle-Based Fluid Simulation for Interactive Applications"では、圧力pを密度ρから求める際に、定常密度ρ0を引いて

p = k(ρ - ρ0)

としている（kは定数）。ρ0を引く理由としてMullerらは、

... Since pressure forces depend on the gradient of the pressure field, the offset mathematically has not effect on pressure forces. However, the offset does influence the gradient of a field smoothed by SPH and makes the simulation numerically more stable. ...

と述べている。

なんのこっちゃ？

オフセット（ρ0）が数学的に圧力に影響を与えないといいつつ、それによってシミュレーションがより安定化するといっている。しかし、数学的に影響を与えないのであれば、安定化するもなにもないではないか。

そこで、引用元である、Desbrunの"Smoothed Particels: A new paradigm for animating highly deformable bodies"の中で、なぜρ0が導入されているのかを見てみた。

Desbrunらによると、

... In astrophysics applications, pressure forces were often combined with gravitational forces balancing the expansion phenomenon.
  In contrast, we would like to animate materials with constant density at rest. Consequently, the material sould exhibit some internal cohesion, resulting in attraction-repulsino forces as in the Lennard-Jones model. ...

とある。

なるほど。

もともとSPHは天体のシミュレーションで用いられていた手法なので、真空中での粒子の運動を対象としている。一方で、Desbrunらが行おうとしている大変形体のシミュレーションは、大気中での運動が対象となる。

大気中での運動なので、本来であれば、大変形体の周りに空気があって、その気圧と大変形体の圧力が釣り合っていないといけない。

それを論文でどう扱っているかというと、周りの空気も粒子としてシミュレーションに含めるのではなく、その代わりに、ρからρ0を引いて気圧の影響を補正した圧力を用いるという方法をとっている。つまり、定常状態での圧力を仮想的に0気圧となるようにオフセットをとって計算している。

つまり、定常状態での圧力を仮想的に0気圧とすることで、気圧と大変形体の圧力を釣り合わせているということになる。

実際に、ρ0を変えてシミュレーションをしてみた。

【ρ0=600】


【ρ0=300】


【ρ0=0】


確かに、ρ0が小さくなるにつれて、粒子が空気中に広がっていっている様子がわかる。ρ0=0に至っては、まさに真空中に拡散しているように見える結果が得られた。

結論
圧力pを密度ρから求める際に定常密度ρ0を引いてp = k (ρ - ρ0)とするのは、気圧の影響を加味して粒子が拡散してしまわないようにするためだった。SPHを宇宙から地球上に持ってくるわけですな。

■関連する記事
粒子法のプログラム第１回（概要）
粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）
【粒子法】粒子を流体としてレンダリング
３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…

ClickTaleの無料利用分を使い切ってしまった

おととい入れてみたClickTale。

２日で1週間分のレコーディング数を使い切ってしまった…。これじゃ使えないなぁ。

１人１人のユーザの動きを追える素晴らしいツールなんだけど、月$99は払えないや。

自分で作ろうか？？

ClickTaleの再生がAnalyticsに計測されてしまわないか？

このブログにClickTaleを入れてみた。

ClickTaleはユーザのマウスの動きをレコーディングしてくれるアクセス解析ツールで、レコーディングされたマウスの動きを再生できる。

再生機能で１つ気になったのが、再生したときのページ読み込みが解析対象ページに埋まっているGoogle Analyticsで計測されてしまわないかという点。

Analyticsで確認してみると、ClickTaleからのアクセスは計測されないみたいでよかった。

OneApp

MicrosoftがOneAppという携帯電話向けのアプリケーションフレームワークを発表した。

OneAppは、途上国の人々がTwitterやFacebookのようなサービスに携帯電話からアクセスできるようにするサービスで、プロセッサが非力でメモリが少なくてもJavaが動きさえすれば使うことができる。

大型Webアプリケーションを動かすために必要な処理の一部をMicrosoftのサーバがクラウドとして代行する仕組みで、デベロッパはJavaScriptとXMLによる比較的簡単な作業でアプリをOneAppに対応させることができる。

Microsoftが途上国向けにサービスを提供する理由には、人道的理由もあるが、途上国が潜在的に巨大な市場であるというビジネス上のものもあるだろう。

iTunesでムービーのサムネイルを変えるには

iTunesのムービーのサムネイルは好きなタイムラインのものに変えることができる。

▼変え方

1. サムネイルにしたいタイミングでムービーを停止する
2. ムービーを右クリック
3. 「ポスターフレームを設定」を選ぶ

これで好きなカットをサムネイルにできる♪

粒子法のシーンを２倍のサイズにしてみた

「粒子法のシーンを２倍のサイズにしてみたが…」で失敗したものをもう一度計算した。

【２倍のサイズ】


【元のサイズ】


初期配置や境界の形はまったく同じにして、xyz各方向に２倍している。一見同じように見えるけど、床のチェック模様の細かさが違うところに注目して欲しい。

２倍サイズにすると、より高いところからさらに多くの流体が落下するため、盛大なしぶきがあがった。上のレンダリング結果だと黒くなってしまってわかりづらいが、陰関数曲面をそのままレンダリングするとその様子がよくわかる。

粒子法(SPH)のプログラム

粒子法(SPH)のプログラムを解説したシリーズです。ソースコードも公開しています。

粒子法(SPH)のプログラム一覧

動画

シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。

動画の一覧

--

【論文】Particle-Based Fluid Simulation for Interactive Applications

Matthias Mullerらの「Particle-Based Fluid Simulation for Interactive Applications」を読んだ。

この論文では、SPHをベースに、自由表面を持つ流体をインタラクティブにシミュレーションしレンダリングする方法を提案している。

この論文におけるMullerらの貢献は、

• Desbrunらの大変形体向けの手法を流体に拡張
• インタラクティビティの実現

の２点。

Desbrunの手法を流体に拡張するにあたり、Mullerらは流体が受ける力をナビエ=ストークス方程式から直接導出するようにし、さらに、表面張力をモデル化した項を加えることを提案している。

また、インタラクティブに流体を描画できるように、高速に計算できるカーネルを設計すると共に、キャッシュヒット率が向上するようデータ構造を工夫している。

レンダリングにはPoint SplattingとMarching Cubesの２通りの手法を用いており、2003年当時のPC(Pentium4+GeForce4)で2000個程度の粒子を計算した例では、Point Splattingで20fps、Marching Cubesで5fpsでの描画を実現している。

ハンニバル戦記

塩野七生の「ローマ人の物語」を読んでいる。

今日は、２巻の「ハンニバル戦記」を読み始めた。ハンニバル戦記は、そのタイトルの通りハンニバルのアルプス越えで有名なポエニ戦役について書かれている。

ポエニ戦役は紀元前264年〜同133年までローマとカルタゴの間で行われた戦争で、今日読んだのはそのうちポエニ戦役のきっかけとなったメッシーナでの戦い。

紀元前265年、シチリア半島の東端に位置するメッシーナは、南に位置するシラクサから攻め込まれていた。

メッシーナは、海峡を挟んで対岸にあるローマへ救援を求める。頼られたローマは迷ったあげく、クラウディウス率いる軍を派遣した。

ローマはシラクサ、カルタゴを撃破。シラクサとは同盟関係を結ぶことでシチリア島に影響力を伸ばすことに成功する。

しかし、これによってカルタゴの勢力圏と直接接することになり、その後のポエニ戦役へと繋がっていく…

次はいよいよポエニ戦役の始まり！

由比ヶ浜

夕方に鎌倉の由比ケ浜までふらっと行ってきた。

鎌倉には何度かいったことがあるけど、由比ヶ浜にいったのは初めて。

行ってみると、広いし、水もきれいだし、設備もあるし、駅から近い。すごく良い！

海というと子供のころの印象が強くて、車で何時間もかけて頑張っていかないといけないというイメージがあったけど、由比ケ浜は全然そんなことがない。

すごくいいところを見つけた♪

粒子法のシーンを２倍のサイズにしてみたが…

３次元の粒子法シミュレーションを、さらにシーンを２倍のサイズにして計算してみた。



あ…カメラの位置を調節していなかった…
再計算します。

追記：再計算しました。こちら


■関連する記事
粒子法のプログラム第１回（概要）
粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）
【粒子法】粒子を流体としてレンダリング
３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…
粒子法のシーンを２倍のサイズにしてみた

３次元の粒子法シミュレーション

粒子法のプログラム第１回（概要）や【粒子法】粒子を流体としてレンダリングでは粒子を２次元で計算していたので、今回はそれを３次元に拡張して計算してみた。平面だったものを、そのままz方向にも粒子を配置している。


【粒子表現】


【流体表現】



■流体シミュレーションに関するエントリ
粒子法のプログラム第１回（概要）
粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）
【粒子法】粒子を流体としてレンダリング
３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…
粒子法のシーンを２倍のサイズにしてみた
SPHによる巻き波のシミュレーション２
Haskell、OCamlでSPH法
このあとやりたいこと
カメラ位置を変えて流体をレンダリング
固液連成シミュレーション

■剛体シミュレーションに関するエントリ
粒子ベース剛体シミュレーション（プレビュー）
粒子ベース多体衝突シミュレーション
引き続き、粒子ベース剛体シミュレーション

粒子法(SPH)のプログラム

粒子法(SPH)のプログラムを解説したシリーズです。ソースコードも公開しています。

粒子法(SPH)のプログラム一覧

動画

シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。

動画の一覧

--

Yahoo! BOSS

米Yahoo!のサービスに、BOSSというのがある。BOSSはウェブ検索サービスを構築するためのオープンプラットフォームで、APIを通じてサードパーティーが自由にYahoo! Searchのインデックスやランキングアルゴリズムなどの技術を使うことができる。

TechCrunchによると、今年の５月現在でBOSSは1日あたり3000万の検索クエリを扱っており、さらに利用量が伸びている。この規模は、Bingと互角かやや抜いているレベルの規模だ。

BOSSを採用している企業には、セマンティック検索エンジンを手掛けるhaikaやソーシャルサービスのリアルタイム検索エンジンOneRiot、より理解しやすい検索結果の提供を目指すCluuzといったものがある。TechCrunchのpower searchにも使われているようだ。

BOSSを使えば大きな投資資金がなくても検索エンジンを作ることができる。その計算資源を使って、何かビジネスができないだろうか？

Haskellでレイトレーシング（第５回）〜タプルを返せる

Haskellでは、関数の戻り値としてタプルを返すことができます。

今回は、タプルを返す場合についてCommon LispとHaskellでの書き方の違いを書きます。

Common Lispでタプルを返す場合は、以下のようにvalues関数で値を返して、multiple-value-bindでその値を受け取ります。

(defun sendray (pt xr yr zr)
  (multiple-value-bind (s int) (first-hit pt xr yr zr)
    (if s
        (* (lambert s int xr yr zr) (surface-color s))
        0)))

(defun first-hit (pt xr yr zr)
  (let ...略...
    (values surface hit)))

first-hitがレイと交差した面とその交点の組を返す関数です。それをsendrayから呼び出して、multiple-value-bindで受け取っています。

一方、Haskellではこのように書きます。

sendray :: [Surface] -> Light -> Point ->
           Double -> Double -> Double -> Color
sendray world light pt xr yr zr =
  case first_hit world pt xr yr zr of
    Just (s@(Sphere (Color col) _ _ _), int)
      -> Color ((lambert s int light) * col)
    Nothing
      -> Color 0.0

first_hit :: [Surface] -> Point ->
             Double -> Double -> Double ->
             Maybe (Surface, Point)
first_hit world pt xr yr zr = nearest (map hit world)
  where nearest :: [Maybe (Double, Surface, Point)]
                   -> Maybe (Surface, Point)
        nearest hits = case foldr cmp Nothing hits of
                         Just (d, s, h) -> Just (s, h)
                         otherwise      -> Nothing
     
        ...略...

first_hitでMaybe (Surface, Point)型を返し、それをsendrayでパターンマッチで受け取っています。Haskellでは、タプルを返すのに関数をはさむ必要はなく、そのまま書くことができます。

次回は、「5.衝突点の計算」について書きます。


■他の記事
Haskellでレイトレーシング（第１回）〜導入
Haskellでレイトレーシング（第２回）〜ループは使わない
Haskellでレイトレーシング（第３回）〜型による分岐にはパターンマッチを用いる
Haskellでレイトレーシング（第４回）〜NilはMaybeモナドで表現する
Haskellでレイトレーシング（第５回）〜タプルを返せる

Haskellでレイトレーシング（第４回）〜NilはMaybeモナドで表現する

今回は、関数の戻り値に正規の値とそうでない値がある場合に、Haskellではどのように表現するのかを書きます。

関数が正規の結果とそうでない結果を返す場合の例として、視点から投げられたレイと空間中の物体の交点を返す関数を考えます。

この関数ではレイと物体の交点があるときにはその交点座標を返します。では、レイと物体が交差しない場合には何を返すのでしょうか？

Common Lispのコードでは、レイと物体が交差しない場合には、（暗黙的に）Nilを返すようになっています。

(defun sphere-intersect (s pt xr yr zr)
  (let* ...略...
    (if n
        (make-point :x (+ (x pt) (* n xr))
                    :y (+ (y pt) (* n yr))
                    :z (+ (z pt) (* n zr))))))

レイと物体の交点を求めるための方程式の解nがある場合はmake-pointで作った交点を返し、解nがない場合には暗黙的にNilを返しています。

実用上はこのようにNilを返す形でも問題ないのですが、HaskellではこれをMaybeモナドを使ってより適切に表現することができます。

intersect :: Surface -> Point -> Double -> Double -> Double -> Maybe Point
intersect (Sphere _ r (Point cx cy cz) _) (Point px py pz) xr yr zr =
    let ...略...
    in case minroot a b c of
         Just n  -> Just (Point (px+xr*n) (py+yr*n) (pz+zr*n))
         Nothing -> Nothing

minrootで求めた方程式の解nがある場合にはJust Pointを返し、解nがない場合にはNothingを返しています。

つまり、関数の戻り値をMaybe Point(交点があるかもしれないよ)型と定義し、交点がある場合にはJust Point(交点はコレ！)を返し、交点がない場合にはNothing(交点がなかった…)を返すのです。

いかがでしょう？やりたいことをすごく自然に表現できているのではないでしょうか？

次回は、「4.タプルを返せる」について書きます。

■他の記事
Haskellでレイトレーシング（第１回）〜導入
Haskellでレイトレーシング（第２回）〜ループは使わない
Haskellでレイトレーシング（第３回）〜型による分岐にはパターンマッチを用いる
Haskellでレイトレーシング（第４回）〜NilはMaybeモナドで表現する
Haskellでレイトレーシング（第５回）〜タプルを返せる

Haskellでレイトレーシング（第３回）〜型による分岐にはパターンマッチを用いる

今回は、パターンマッチを使った分岐について書きます。

Common Lispのコードでは、レイとSurfaceの交点を計算するときに、typecase関数を使ってSufaceの型によって処理を分岐させています。分岐用の関数を用意しなければなりません。

(defun intersect (s pr xr yr zr)
 (funcall (typecase s (sphere #'sphere-intersect))
          s pt xr yr zr))

(defun sphere-intersect (s pr xr yr zr)
 ...)

一方、Haskellではパターンマッチを使うことで、型による分岐を自然に表現することができます。

intersect :: Surface -> Point -> Double -> Double -> Double -> Maybe Point
intersect (Sphere _ r (Point cx cy cz)) (Point px py pz) xr yr zr =
   let a = (sq xr) + (sq yr) + (sq zr)
       b = 2 * ((px-cx)*xr + (py-cy)*yr + (pz-cz)*zr)
       c = (sq (px-cx)) + (sq (py-cy)) + (sq (pz-cz)) - (sq r)
       n = minroot a b c
   in case n of
        Just n  -> Just (Point (px+xr*n) (py+yr*n) (pz+zr*n))
        Nothing -> Nothing

今回はSphereのみですが、その他の型を追加した場合にはマッチさせるパターンを追記するだけです。

このように、Haskellでは引数のパターンによって分岐することができるのです。

次回は、「3.NilはMaybeモナドで表現する」について書きます。

■他の記事
Haskellでレイトレーシング（第１回）〜導入
Haskellでレイトレーシング（第２回）〜ループは使わない
Haskellでレイトレーシング（第３回）〜型による分岐にはパターンマッチを用いる
Haskellでレイトレーシング（第４回）〜NilはMaybeモナドで表現する
Haskellでレイトレーシング（第５回）〜タプルを返せる

Haskellでレイトレーシング（第２回）〜ループは使わない

いま、Haskellでレイトレーサを書いています。

第１回のブログで、Common LispのコードをHaskellで書くときにポイントとなったところをあげました。

今回は、そのポイントのうち、「1.ループを使わない」について書きます。

ループを使わない

各ピクセルごとに色を計算するために、Common Lispのコードでは以下のように手続き的にループを回しています。

(defun tracer (pathname &optional (res 1))
 (with-open-file (p pathname :direction :output)
   (format p "P2 ~A ~A 255" (* res 100) (* res 100))
     (let ((inc (/ res)))
       (do ((y -50 (+ y inc)))
           ((< (- 50 y) inc))
         (do ((x -50 (+ x inc)))
             ((< (- 50 x) inc))
           (print (color-at x y) p)))

一方、Haskellでは手続き的なループは使わないので、このような書き方はしません。

それではどのように書くのかというと、以下のようにリストで表現するのです。

tracer :: FilePath -> [Surface] -> Int -> IO()
tracer pathname world res =
   let
       cols = [color_at world (pos i) (pos j)  | j <- [0..(res*100)-1],                                                   
                                                 i <- [0..(res*100)-1]]
   in
       do putPGM pathname (res*100) (res*100) cols
   where
         putPGM :: FilePath -> Int -> Int -> [Int] -> IO()
         putPGM pathname w h cols =
             let header = "P2 " ++ show w ++ " " ++ show h ++ " 255\n"
                 body   = foldr (++) "" [(show c) ++ "\n" | c <- cols]
             in do writeFile pathname (header ++ body)
         
         pos :: Int -> Double
         pos i = (fromIntegral i) / (fromIntegral res) - 50.0

putPGMがpgm画像を出力するアクションで、それにピクセルの色のリストcolsを渡すのです。

このような書き方をすると、「すべてのピクセルをリストにしてから処理するなんて、メモリを無駄に使いすぎるのではないか」と思ってしまいます。しかし、Haskellは遅延評価なので、必要な要素から順々に計算されていくため、リスト全体をメモリに持っておくということは起きません。

このように、リストと遅延評価を使うことで、手続き型のループと同様の処理をHaskellで表現することができるのです。

次回は、「2.型による分岐にはパターンマッチを用いる」について書きます。

■他の記事
Haskellでレイトレーシング（第１回）〜導入
Haskellでレイトレーシング（第２回）〜ループは使わない
Haskellでレイトレーシング（第３回）〜型による分岐にはパターンマッチを用いる
Haskellでレイトレーシング（第４回）〜NilはMaybeモナドで表現する
Haskellでレイトレーシング（第５回）〜タプルを返せる

Haskellでレイトレーシング（第１回）〜導入

Haskellでレイトレーサを書き始めました。ベースにしているのは、「ANSI Common Lisp」という本に載っていたレイトレーサのサンプルです。Common Lispで書かれているものをHaskellで書いています。

Haskellで書いたコードはこれで、実行すると以下のようなpgm画像が得られます。

Common LispのコードをHaskellで書くにあたっては、次のようなところがポイントでした。

1.ループを使わない

Common Lispのコードでは、各ピクセルごとに手続き的にループを回すようになっていますが、Haskellではそれをループを使わずに記述します。どのように書くのかというと、全ピクセルをリストで表現し、それを関数に渡すのです。Haskellの遅延評価が活きてくるポイントです。

2.型による分岐にはパターンマッチを用いる

Common Lispのコードでは、物体の型を判別して処理を分岐しています。Haskellでは、引数のパターンマッチでこれをもっと自然に表現できます。

3.NilをMaybeモナドで表現する

レイが物体と衝突しない場合に、Commpn Lispのコードでは暗黙的にNilを返すように記述されています。それを、Maybeモナドを使いNothingで表現するようにしました。レイが物体に衝突するときにはその交点を返し、衝突しないときにはNothingを返します。

4.タプルを返せる

Common Lispのコードではmultiple-value-bindを使って値の組を返す部分があります。Haskellではそのままタプルを返すことができます。

5.衝突点の計算

Common Lispのコードでは、もっとも視点に近い衝突点を探すためにループが用いられています。Haskellでは、この衝突点の計算もループを使わずに表現します。foldrを使って、Haskellらしい形で表現しました。

次回以降の投稿で、これらのポイントを１つずつ説明していきます。

■他の記事
Haskellでレイトレーシング（第１回）〜導入
Haskellでレイトレーシング（第２回）〜ループは使わない
Haskellでレイトレーシング（第３回）〜型による分岐にはパターンマッチを用いる
Haskellでレイトレーシング（第４回）〜NilはMaybeモナドで表現する
Haskellでレイトレーシング（第５回）〜タプルを返せる

メアリー・ブレア展

現代美術館(MOT)で開かれている、メアリー・ブレア展に行ってきました。メアリー・ブレアは、ディズニー映画のコンセプト・アートを手がけた女性です。

「シンデレラ」や「不思議の国のアリス」、「ピーターパン」などのコンセプト・アートが展示されて、映画の雰囲気がギュッと凝縮されていました。メアリーの絵の特徴は鮮やかで強烈な色使いで、マティスやデュフィを連想させます。

こういう絵が好き♪

世界のへそ

今日のBingトップはバリ島のアグン山。そのHotspotの中に、

バリ島北部に位置するアグン山。ヒンドゥー教で世界の中心と考えられていることから「世界のへそ」と呼ばれています。

とあった。

これを見て「なんでヒンドゥー教で世界の中心と考えられている山がバリ島にあるの？」と思い、調べてみた。

検索してみたところ、ここによると、

バリ島では、インドネシアのほかの島々とはちがって、島民の93.2％が、ヒンドゥー教を信仰している。これに、アグン山信仰に代表されるような、むかしから伝わる土着の自然信仰がくわわって、独特のバリ・ヒンドゥー教世界がつくられた。

とある。

また、Wikipediaによると、

バリ島では、クディリ朝の支配下に入った11世紀初めごろからヒンドゥー・ジャワ文化の影響が及び始め、その後しばらくジャワの支配を離れるが、マジャパヒト王国がバリを征服した1343年以後、16世紀初めにジャワのイスラム化によって同王国が滅亡するまでにヒンドゥー化が広く浸透した。

とある。

どうやら、バリ島には11世紀初めにヒンドゥー教が渡来してきて、それが土着の自然信仰と混ざって、インドとは異なるバリ・ヒンドゥー教ができたらしい。そして、その中にアグン山信仰も受け継がれているそう。

バリ島が独特なのは、こういうことだったのね。

Electronics Manufacturing Service

Appleは、自分で工場を持ってせっせとiPodを作っているわけではなくて生産を外部に委託している。そういった生産設備を持たない会社が製品の生産を委託する先の会社って何というのだろうと思っていた。

シリコンチップの場合はファブレス(fubless)とファウンダリ(foundary)という単語で調べられるが、電子機器の場合はどうも呼び方が違うようで検索しても見つからない。

では、そういう会社は何と呼ばれているかというと、EMS(Electronics Manufacturing Service)というらしい。

大手では、

• ホンハイ（台湾）
• フレクストロニクス（シンガポール）
• ソレクトロン（アメリカ）
• ウィストロン（台湾）
• セレスティカ（カナダ）
• シークス（日本）

といった会社がある。やっぱり台湾が強いのか、２社が台湾にある。あとはアメリカ、カナダ、シンガポール。日本にもあるみたい。

発注元と受託先の役割分担としてはおおまかに２種類あって、

• 設計は発注元が行い生産を受託するOEM
• 上流工程である設計工程も含めて受託するODM(original design manufacturer)

があるらしい。ODMでは設計もやる。

EMSについて疑問に思ったことは、

• 製品の仕様を用意すれば内部設計なんかはやってくれるのだろうか？
• こういう会社と協力すればベンチャーがプロダクトを売ることができるのでは？
• ロットってどれくらいで受け付けるんだろう？
• たとえば毎年新しいのがたくさんでてくるおもちゃなんかも、こういう会社で生産されているんだろうか？

粒子法のプログラム最終回（粒子の出力）

最終回となる今回は、計算した粒子の出力について説明します。

このプログラムでは、計算した粒子をpov-rayのフォーマットで出力しています。

第１回に載せたムービーは、このプログラムでタイムステップごとに出力した.povファイルをレンダリングし、さらにffmpegで動画にしたものです。

粒子を出力するコードは以下のようになっています。

void output_particles( Particles* ps )
{
static int num = 0;
std::string file_name = \
(boost::format( "result%08d.pov" ) % num++ ).str();
std::cout << "processing " << file_name << " ..." << std::endl;
std::ofstream f( file_name.c_str() );
if ( ! f )
{
std::cerr << "cannot open " << file_name << std::endl;
exit(1);
}

f << "#include ¥"colors.inc¥"¥n"
<< "camera {¥n"
<< "  location <10,>¥n"
<< "  look_at <10,>¥n"
<< "}¥n"
<< "light_source { <0,> color White }¥n";
FOR_EACH_PARTICLES( ps, p )
{
f << "sphere {¥n"
<< "  <" << r(0) << ", " << r(1) << ", " << r(2) << ">, 0.5¥n"
<< "  texture {¥n"
<< "    pigment { color Yellow }¥n"
<< "  }¥n"
<< "}¥n";
}
f << std::endl;
}

タイムステップごとに.povファイルを出力しています。

さらに、この投稿に載せたムービーでは、計算した粒子から求めた陰曲面をマーチングキューブでポリゴン化し、反射や屈折も考慮してレンダリングしています。

これでひととおりの説明が終わりました。このあとの発展としては、3次元にしたり、物体との相互作用を入れたり、GPUで処理したりといったことを考えています。

粒子法のプログラム

粒子法のプログラム第１回（概要）
粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）

その他の解説エントリ

SPHによる巻き波のシミュレーション１
SPHによる巻き波のシミュレーション２
SPHによる巻き波のシミュレーション３
このあとやりたいこと

固液連成シミュレーションに関するエントリ

粒子法による固液連成シミュレーション

流体シミュレーションに関するエントリ

【粒子法】粒子を流体としてレンダリング
３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…
粒子法のシーンを２倍のサイズにしてみた
Haskell、OCamlでSPH法
カメラ位置を変えて流体をレンダリング
Bunny-shaped fluid simulation

剛体シミュレーションに関するエントリ

粒子ベース剛体シミュレーション（プレビュー）
粒子ベース多体衝突シミュレーション
引き続き、粒子ベース剛体シミュレーション
Falling Rigid Bunnies

動画

シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。

動画の一覧

--

粒子法のプログラム第６回（境界条件と粒子位置の更新）

今回は、粒子の境界条件と粒子位置の更新について説明します。

境界条件と粒子位置の更新のコードは以下のようになります。

void calculate_position( Particles* ps )
{
vec accel(3), g(3), norm(3);
double speed, diff, adj;

g(0) = 0.0; g(1) = -9.8; g(2) = 0.0;
FOR_EACH_PARTICLES( ps, p )
{
accel = p->f * SPH_PMASS;

speed = norm2( accel );
if ( speed > SPH_LIMIT*SPH_LIMIT ) {
accel *= SPH_LIMIT / sqrt(speed);
}

// Z-axis walls
diff = 2.0 * RADIUS - ( p->r(2) - MIN(2) ) * SPH_SIMSCALE;
if ( diff > EPSILON )
{
norm(0) = 0.0; norm(1) = 0.0; norm(2) = 1.0;
adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP * inner_prod( norm, p->v );
accel += adj * norm;
}
diff = 2.0 * RADIUS - ( MAX(2) - p->r(2) ) * SPH_SIMSCALE;
if ( diff > EPSILON )
{
norm(0) = 0.0; norm(1) = 0.0; norm(2) = -1.0;
adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP * inner_prod( norm, p->v );
accel += adj * norm;
}

// X-axis walls
diff = 2.0 * RADIUS - ( p->r(0) - MIN(0) ) * SPH_SIMSCALE;
if ( diff > EPSILON )
{
norm(0) = 1.0; norm(1) = 0.0; norm(2) = 0.0;
adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP * inner_prod( norm, p->v );
accel += adj * norm;
}
diff = 2.0 * RADIUS - ( MAX(0) - p->r(0) ) * SPH_SIMSCALE;
if ( diff > EPSILON )
{
norm(0) = -1.0; norm(1) = 0.0; norm(2) = 0.0;
adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP * inner_prod( norm, p->v );
accel += adj * norm;
}

// Y-axis walls
diff = 2.0 * RADIUS - ( p->r(1) - MIN(1) ) * SPH_SIMSCALE;
if ( diff > EPSILON )
{
norm(0) = 0.0; norm(1) = 1.0; norm(2) = 0.0;
adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP * inner_prod( norm, p->v );
accel += adj * norm;
}
diff = 2.0 * RADIUS - ( MAX(1) - p->r(1) ) * SPH_SIMSCALE;
if ( diff > EPSILON )
{
norm(0) = 0.0; norm(1) = -1.0; norm(2) = 0.0;
adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP * inner_prod( norm, p->v );
accel += adj * norm;
}

accel += g;
p->v += accel * DT;
p->r += p->v * DT / SPH_SIMSCALE;

// stick on x-y plane
p->r(2) = 0.0;

}
}

それぞれの粒子について、粒子にかかる力から加速度を計算したあと、各壁面について境界条件を適用し、重力加速度を加えます。そうやって求まった加速度から、粒子の速度を求め、粒子の位置を更新します。

境界条件には、ペナルティ法を用いています。ペナルティ法とは、粒子が壁面へめりこんだ場合に、そのめり込みの程度によって壁面との衝突による反発力を与える方法です。

このようにして、境界条件を処理し、粒子の位置を更新しています。

次回は計算結果の出力です。最終回です。

粒子法のプログラム

粒子法のプログラム第１回（概要）
粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）

その他の解説エントリ

SPHによる巻き波のシミュレーション１
SPHによる巻き波のシミュレーション２
SPHによる巻き波のシミュレーション３
このあとやりたいこと

固液連成シミュレーションに関するエントリ

粒子法による固液連成シミュレーション

流体シミュレーションに関するエントリ

【粒子法】粒子を流体としてレンダリング
３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…
粒子法のシーンを２倍のサイズにしてみた
Haskell、OCamlでSPH法
カメラ位置を変えて流体をレンダリング
Bunny-shaped fluid simulation

剛体シミュレーションに関するエントリ

粒子ベース剛体シミュレーション（プレビュー）
粒子ベース多体衝突シミュレーション
引き続き、粒子ベース剛体シミュレーション
Falling Rigid Bunnies

動画

シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。

動画の一覧

--

Vivienne Westwoodの長財布

財布って、どんなのを使ってます？

私は父親からアメリカ土産にもらった革の三つ折り財布を10年以上使っているのですが、最近、長財布もいいなぁと思っていました。

それで今日、表参道に行って、Vivienne Westwoodの長財布を買ってきました。

買ったのは、写真にあるADVANのネイビー。表面に味のあるアドバン加工を施したレザーをベースに、全面に大胆なORB型の型押しを施してあります。ORBはVivienneのいろんなアイテムに施されていて、そこには「伝統をもって未来を創る」という彼女のコンセプトが込められています。

あと、外側のネイビーと内側のイエローのコンビネーションがとってもクール！この絶妙なコンビネーションが素晴らしい！！

ブティックで買い物をしたのが初めてだったので、ドキドキしたー

大事に長く使おう！

粒子法のプログラム第５回（力の計算）

今回は、粒子法のプログラムでどのように力を計算しているかを説明します。

力の計算は次のコードになります。

void calculate_force( Particles* ps )
{
double pterm, vterm, r, c;
vec dr(3), force(3), fcurr(3);

FOR_EACH_PARTICLES( ps, p )
{
force(0) = 0.0; force(1) = 0.0; force(2) = 0.0;
FOR_EACH_PARTICLES( ps, pj )
{
if ( p == pj ) continue;
dr = ( p->r - pj->r ) * SPH_SIMSCALE;
r = norm( dr );
if ( H > r )
{
c = H - r;
pterm = -0.5 * c * SpikyKern * ( p->p + pj->p ) / r;
vterm = LapKern * SPH_VISC;
fcurr = pterm * dr + vterm * ( pj->v - p->v );
fcurr *= c * p->rho * pj->rho;
force += fcurr;
}
}
p->f = force;
}
}

それぞれの粒子について、ナビエ=ストークス方程式の圧力項、粘性項にあたる力を計算しています。

影響半径内にある自分以外の粒子に対して、距離に応じて重みを足し合わせています。SpikyKernとLapKernがカーネルの係数です。

力の計算も、密度と圧力の計算と同様に、近傍粒子探索を工夫したり対称性を利用することで、高速化できます。

次回は、どのように粒子の位置を更新するかという説明になります。

粒子法のプログラム

粒子法のプログラム第１回（概要）
粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）

その他の解説エントリ

SPHによる巻き波のシミュレーション１
SPHによる巻き波のシミュレーション２
SPHによる巻き波のシミュレーション３
このあとやりたいこと

固液連成シミュレーションに関するエントリ

粒子法による固液連成シミュレーション

流体シミュレーションに関するエントリ

【粒子法】粒子を流体としてレンダリング
３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…
粒子法のシーンを２倍のサイズにしてみた
Haskell、OCamlでSPH法
カメラ位置を変えて流体をレンダリング
Bunny-shaped fluid simulation

剛体シミュレーションに関するエントリ

粒子ベース剛体シミュレーション（プレビュー）
粒子ベース多体衝突シミュレーション
引き続き、粒子ベース剛体シミュレーション
Falling Rigid Bunnies

動画

シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。

動画の一覧

--

目覚めを良くするための最良の方法

最近、枕を変えてみた。

これまではふつうのやわらかい枕。それを、四角くて固めのい草枕に変えた。

すると、目覚めがすごくすっきりするようになった。

枕でここまで変わるのかー

びっくり。これから蒸し暑くなるので、その点でもグッド！

粒子法のプログラム第４回（密度と圧力の計算）

今回は、粒子法のプログラムでどのように密度と圧力を計算しているかを説明します。

粒子の密度と圧力の計算は次のコードになります。

void calculate_density_and_pressure( Particles* ps )
{
double r2, c, sum, H2;
vec dr;

H2 = H*H;

FOR_EACH_PARTICLES( ps, p )
{
sum = 0.0;
FOR_EACH_PARTICLES( ps, pj )
{
if ( p == pj ) continue;
dr = (p->r - pj->r) * SPH_SIMSCALE;
r2 = norm2( dr );
if ( H2 > r2 )
{
c = H2 - r2;
sum += c * c * c;
}
}
p->rho = sum * SPH_PMASS * Poly6Kern;
p->p   = ( p->rho - SPH_RESTDENSITY ) * SPH_INTSTIFF;
p->rho = 1.0 / p->rho;  // take inverse for later calculation
}
}

それぞれの粒子について、影響半径内にある自分以外の粒子に対して、距離に応じて重みを足し合わせています。そのあと、足し合わせた重みに、粒子質量についての係数とカーネル係数をかけています。密度が計算された後、圧力を密度から計算します。

上のやり方だと、計算量のオーダーは、各粒子について他の粒子それぞれとの計算になるので、粒子の数をnとするとO(n^2)となります。しかし、影響半径外の粒子との計算は不要なので、他粒子すべてと計算する必要は本来ありません。なので、近傍粒子をあらかじめマップに入れておき、それを参照するようにすれば、計算量をO(n)にすることができます。

また、粒子pと粒子pjについて対称なので、その対称性を利用すれば、計算量をさらに半分にすることが可能です。

密度と圧力は、このように計算されています。いったん素直な実装になっていますが、上の２つの点で高速化することができます。

次回は、粒子にかかる力をどのように計算しているかという説明です。

粒子法のプログラム

粒子法のプログラム第１回（概要）
粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）

その他の解説エントリ

SPHによる巻き波のシミュレーション１
SPHによる巻き波のシミュレーション２
SPHによる巻き波のシミュレーション３
このあとやりたいこと

固液連成シミュレーションに関するエントリ

粒子法による固液連成シミュレーション

流体シミュレーションに関するエントリ

【粒子法】粒子を流体としてレンダリング
３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…
粒子法のシーンを２倍のサイズにしてみた
Haskell、OCamlでSPH法
カメラ位置を変えて流体をレンダリング
Bunny-shaped fluid simulation

剛体シミュレーションに関するエントリ

粒子ベース剛体シミュレーション（プレビュー）
粒子ベース多体衝突シミュレーション
引き続き、粒子ベース剛体シミュレーション
Falling Rigid Bunnies

動画

シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。

動画の一覧

--

【粒子法】粒子を流体としてレンダリング

粒子法では流体を粒子として表現します。

そのため、そのまま描画すると、この投稿のように球のあつまりとして描画することになります。

これを、流体として見えるようにレンダリングしてみました。

マーチングキューブで陰曲面をポリゴン化し、レイトレーシングでレンダリングしています。

フォトンマップも使っているのでコースティクスも表現されています。

屈折率を高めにしたのでキラキラ。




■固液連成シミュレーションに関するエントリ
粒子法による固液連成シミュレーション

■関連する記事
粒子法のプログラム第１回（概要）
粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）
【粒子法】粒子を流体としてレンダリング
３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…
粒子法のシーンを２倍のサイズにしてみた

粒子法(SPH)のプログラム

粒子法(SPH)のプログラムを解説したシリーズです。ソースコードも公開しています。

粒子法(SPH)のプログラム一覧

動画

シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。

動画の一覧

--

粒子法のプログラム第３回（データ構造）

今回は、粒子法のプログラムのデータ構造を説明します。

使用したデータ構造は、Particle構造体とParticlesベクタの２つです。

Particles構造体は以下です。

typedef boost::numeric::ublas::vector<double> vec;
typedef struct
{
vec    r;
vec    v;
double rho;
double p;
vec    f;
} Particle;

boostライブラリのublasを使い、３次元ベクトルを表現しています。
rは粒子座標、vは粒子速度、rhoは粒子地点での密度、pは粒子地点での圧力、fは粒子が受ける力です。

Particlesベクタは以下です。

typedef std::vector<particle> Particles;

STLのvectorコンテナでParticle構造体のベクタを表現しています。

使用しているデータ構造は、この２つのみです。このように、粒子法はデータ構造が非常に簡単なのが特徴です。

次回からは、粒子法のアルゴリズムの説明になります。

粒子法のプログラム

粒子法のプログラム第１回（概要）
粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）

その他の解説エントリ

SPHによる巻き波のシミュレーション１
SPHによる巻き波のシミュレーション２
SPHによる巻き波のシミュレーション３
このあとやりたいこと

固液連成シミュレーションに関するエントリ

粒子法による固液連成シミュレーション

流体シミュレーションに関するエントリ

【粒子法】粒子を流体としてレンダリング
３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…
粒子法のシーンを２倍のサイズにしてみた
Haskell、OCamlでSPH法
カメラ位置を変えて流体をレンダリング
Bunny-shaped fluid simulation

剛体シミュレーションに関するエントリ

粒子ベース剛体シミュレーション（プレビュー）
粒子ベース多体衝突シミュレーション
引き続き、粒子ベース剛体シミュレーション
Falling Rigid Bunnies

動画

シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。

動画の一覧

--

さくらの根元に・・・

今日は日本橋にいってきました。

通りの両側にさくらが並んでいて、３分咲きくらいになっていたのですが、その根元に、今の時期絶対ないはずのものが。

なんと、雪があったんです。

もちろん降った雪が積もっていたわけではなく、たぶん、桜の開花を遅らせるためなんだろうなと思います。

聞くこと

最近、「聞く」ことに注意している。

これまでも、もちろん人の話を聞いて相手が何を言おうとしているのか推測して返すようにはしていたのだけれど、どうも上手くいっていなかった。

なんで上手くいかないんだろうと思っていると、その原因は推測がずれていることだということに気づいた。

それ以来、人の話を聞くときには、自分の推測があっているか、「それってこういうこと？」と確認するようにしている。

そうすると、びっくりするくらい会話がスムーズになった！
すばらしい。

リボンが流行り

今日、新宿に行ったのでついでにルミネを見てきた。

ルミネでお店を見てみると、リボン付きの服をよくみかけた。

つい最近日経MJでリボンが流行っているよという読んだことがあったので、「確かに流行ってる！」と納得。

フィール・ザ・カライ・センセーション！

シンガポールのマクドナルドのテレビCM。
http://www.mcdonalds.com.sg/tvcs/wasabi_tvc.mov

台詞が日本語なんだけど、そのあとのメッセージがまるっきりカタカナ英語で、めっちゃウケる。ワサビ＝日本＝カタカナ英語。ちょっと！バカにされてるよ！

あと「カライ・センセーション」というのも謎だけど、確かに、ワサビの辛さって、hotでもspicyとも違って、"karai"んだなぁと。

チンチョウゲ

オフィスの前にチンチョウゲの花が咲いています。
チンチョウゲが子供のころから好きです。

子供のころ住んでいたマンションの下に植わっており、２月３月になるとあの特徴的な香りといっしょに咲いて、

けど、１年に１度しか咲かないので、
「ちんちょうげ」という名前がなかなか覚えられず、

何か不思議なにおい、として記憶に残っています。

はるみ

はるみというみかんが売ってたので、食べてみた。

薄皮をむかなくても食べられて、ふつうのみかんより大きいから食べごたえもあって、ちゃんと甘いので、すごくおいしい！

３月８日 東京都美術館〜上野公園〜大崎

昨日は昼から東京都美術館へ行って、上野公園を散歩したあと、大崎のスタバで新卒採用のブログを書いた。

東京都美術館では、「アート＆クラフツ」展をみた。展示のテーマは、１９席後半にイギリスで起こった「アーツ＆クラフツ」というデザイン運動。ウィリアム・モリスが主導した「アーツ＆クラフツ」というのは、手仕事の良さを見直し、自然や伝統から美を再発見し、シンプルなライフスタイルを提案する、という運動で、その背景には、産業革命以降の工業化によって粗悪な大量生産品が増えたことがあったそう。

東京都美術館から出たあと、上野公園を散歩した。あいにくの曇り空だったが、ふと見上げると、さくらが咲いていて、鳥が花の周りを飛び回っていた。そうか、３月か、と思った。空のくもり色と、逆光の枝の黒、鳥の影、そのなかに、咲き始めのさくらのほのかなピンク色が印象的で、ノートにそれを描いた。

そのあとは、大崎へいって、新卒採用のブログを書いた。就活中の学生に、うちの会社へ興味を持ってもらうためのブログ。ターゲットとシナリオを考えた結果、どんな人間がこの会社にいるのかを知ってもらうのがベストだと考え、徹底的に自分のこれまで歩んできた道を話すことにした。ちょっとした自叙伝？

充実した１日だった。感謝！

白いたいやき

白いたいやきが流行っているらしい。

白い理由は、タピオカでんぷんを使っているから。

大福みたいなもちもちな食感だって。

そういえば、白いイチゴというのも、こないだみたなー。

白いイチゴは、突然変異でたまたまできたのを、たくさん作れるようにしたんだって。

赤いイチゴとならべて、紅白で売ってたりするのがいい感じ。