lattice Boltzmann method 20100914 上手く行かない
/**
* Copyright kuma360 ( http://wonderfl.net/user/kuma360 )
* MIT License ( http://www.opensource.org/licenses/mit-license.php )
* Downloaded from: http://wonderfl.net/c/rvYo
*/
//lattice Boltzmann method
//参考
//http://homepage.mac.com/mike1336/md/app/j351_400/393_latticeBoltzmannMethodP/latticeBoltzmannMethod.html
//http://en.wikipedia.org/wiki/Lattice_Boltzmann_methods
//http://www.jstage.jst.go.jp/article/kaigan/65/1/65_56/_article
//http://hexadrive.jp/index.php?id=25
//http://www.jstage.jst.go.jp/article/kaigan/65/1/256/_pdf
//たぶんいろいろ間違ってます。すぐ発散するし。
package
{
import flash.display.Bitmap;
import flash.display.BitmapData;
import flash.display.Sprite;
import flash.events.Event;
import flash.events.MouseEvent;
[SWF(frameRate='20')];
public class Main extends Sprite {
public static var _canvas:BitmapData ;
private var L:Lattice = new Lattice ;
private var P:ParticleManager = new ParticleManager ;
public function Main():void {
_canvas = new BitmapData ( 465 , 465 , false , 0 ) ;
addChild ( new Bitmap ( _canvas ) ) ;
stage.addEventListener (
MouseEvent.CLICK ,
function () :void {
L.reset () ;
}
);
addEventListener (
Event.ENTER_FRAME ,
function ( e:Event ):void {
L.run () ;
P.run ( L ) ;
var x:int = ( Math.floor ( mouseX / _LS ) + _size ) % _size ;
var y:int = ( Math.floor ( mouseY / _LS ) + _size ) % _size ;
L._v[ x ][ y ].i[1].f = 8 ;
}
) ;
}
}
}
import flash.geom.Rectangle;
import flash.text.engine.BreakOpportunity;
//枝
class iDirc {
//f 粒子の数
public var f:Number = 0 ;
//ft 計算後の粒子の数
public var ft:Number = 0 ;
//feq 粘性
public var feq:Number = 0 ;
}
//格子上の点
class Cell {
//ρ 密度
public var p:Number = 0 ;
//u 流速
public var ux:Number = 0 ;
public var uy:Number = 0 ;
//方向
public var i:Vector.<iDirc> = new Vector.<iDirc> ;
public function Cell () {
i.push ( null ) ;
i.push ( new iDirc ) ;
i.push ( new iDirc ) ;
i.push ( new iDirc ) ;
i.push ( new iDirc ) ;
i.push ( new iDirc ) ;
i.push ( new iDirc ) ;
i.push ( new iDirc ) ;
i.push ( new iDirc ) ;
i.push ( new iDirc ) ;
}
//並進過程
public function step1 ( L:Lattice , X:int , Y:int ) :void {
var I:int = 0 ;
var X1:int = X - 1 ;
var X2:int = X ;
var X3:int = X + 1 ;
var Y1:int = Y - 1 ;
var Y2:int = Y ;
var Y3:int = Y + 1 ;
X1 = ( X1 + _size ) % _size ;
Y1 = ( Y1 + _size ) % _size ;
X3 = ( X3 + _size ) % _size ;
Y3 = ( Y3 + _size ) % _size ;
i[1].ft = L._v[X2][Y2].i[1].f;
i[2].ft = L._v[X3][Y2].i[2].f;
i[3].ft = L._v[X1][Y2].i[3].f;
i[4].ft = L._v[X2][Y1].i[4].f;
i[5].ft = L._v[X2][Y3].i[5].f;
i[6].ft = L._v[X3][Y3].i[6].f;
i[7].ft = L._v[X1][Y3].i[7].f;
i[8].ft = L._v[X3][Y1].i[8].f;
i[9].ft = L._v[X1][Y1].i[9].f;
}
//衝突過程
public function step2 () :void {
var I:int = 0 ;
//ρ 密度
p = 0 ;
for ( I = 1 ; I < 10 ; ++ I ) {
p += i[I].ft ;
}
if ( p > 0 ) {
//u 流速
ux = 0 ;
ux += i[2].ft ;
ux += i[6].ft ;
ux += i[8].ft ;
ux -= i[3].ft ;
ux -= i[7].ft ;
ux -= i[9].ft ;
ux /= p ;
//u 流速
uy = 0 ;
uy += i[5].ft ;
uy += i[6].ft ;
uy += i[7].ft ;
uy -= i[4].ft ;
uy -= i[8].ft ;
uy -= i[9].ft ;
uy /= p ;
//
var vv:Number = 1.5 * ( ux * ux + uy * uy ) ;
//※よくわからない
var v2:Number = ux;
var v3:Number = -ux;
var v4:Number = -uy;
var v5:Number = uy;
var v6:Number = ux+uy;
var v7:Number = -ux+uy;
var v8:Number = ux-uy;
var v9:Number = -ux-uy;
//衝突
i[1].feq = (4 / 9) * p * ( 1 - vv ) ;
i[2].feq = (1 / 9) * p * ( 1 - vv + 3 * v2 + 4.5 * v2 * v2 ) ;
i[3].feq = (1 / 9) * p * ( 1 - vv + 3 * v3 + 4.5 * v3 * v3 ) ;
i[4].feq = (1 / 9) * p * ( 1 - vv + 3 * v4 + 4.5 * v4 * v4 ) ;
i[5].feq = (1 / 9) * p * ( 1 - vv + 3 * v5 + 4.5 * v5 * v5 ) ;
i[6].feq = (1 / 36) * p * ( 1 - vv + 3 * v6 + 4.5 * v6 * v6 ) ;
i[7].feq = (1 / 36) * p * ( 1 - vv + 3 * v7 + 4.5 * v7 * v7 ) ;
i[8].feq = (1 / 36) * p * ( 1 - vv + 3 * v8 + 4.5 * v8 * v8 ) ;
i[9].feq = (1 / 36) * p * ( 1 - vv + 3 * v9 + 4.5 * v9 * v9 ) ;
} else {
i[1].feq = 0 ;
i[2].feq = 0 ;
i[3].feq = 0 ;
i[4].feq = 0 ;
i[5].feq = 0 ;
i[6].feq = 0 ;
i[7].feq = 0 ;
i[8].feq = 0 ;
i[9].feq = 0 ;
}
}
//緩和
public function step3 () :void {
//※1/Γ * Fieqの所だと思うけど よくわからない
var I:int = 0 ;
for ( I = 1 ; I < 10 ; ++ I ) {
i[I].f = i[I].ft - _OMEGA * ( i[I].ft - i[I].feq ) ;
}
}
}
//格子
class Lattice {
//ループ用
private var X:int = 0 ;
private var Y:int = 0 ;
//格子
public var _v:Vector.< Vector.<Cell> > ;
public function Lattice () {
_v = new Vector.< Vector.<Cell > > ;
for ( X = 0 ; X < _size ; ++ X ) {
var V:Vector.<Cell> = new Vector.<Cell> ;
for ( Y = 0 ; Y < _size ; ++ Y ) {
V.push ( new Cell ) ;
}
_v.push ( V ) ;
}
reset ( ) ;
}
//動作
public function run () :void {
for ( X = 0 ; X < _size ; ++ X ) {
for ( Y = 0 ; Y < _size ; ++ Y ) {
_v[X][Y].step1( this , X , Y ) ;
}
}
for ( X = 0 ; X < _size ; ++ X ) {
for ( Y = 0 ; Y < _size ; ++ Y ) {
_v[X][Y].step2() ;
}
}
for ( X = 0 ; X < _size ; ++ X ) {
for ( Y = 0 ; Y < _size ; ++ Y ) {
_v[X][Y].step3() ;
}
}
Main._canvas.fillRect ( Main._canvas.rect , 0 ) ;
for ( X = 0 ; X < _size ; ++ X ) {
for ( Y = 0 ; Y < _size ; ++ Y ) {
Main._canvas.fillRect ( new Rectangle ( X * _LS , Y * _LS , _LS , _LS ) , ( ( _v[X][Y].ux * 70 + 120 ) << 16 ) | ( ( _v[X][Y].p * 70 + 120 ) << 8 ) | ( ( _v[X][Y].uy * 70 + 120 ) ) ) ;
}
}
}
public function reset ( ):void {
for ( X = 0 ; X < _size ; ++ X ) {
for ( Y = 0 ; Y < _size ; ++ Y ) {
_v[X][Y].i[1].f = 4/9 ;
_v[X][Y].i[2].f = 1/9 ;
_v[X][Y].i[3].f = 1/9 ;
_v[X][Y].i[4].f = 1/9 ;
_v[X][Y].i[5].f = 1/9 ;
_v[X][Y].i[6].f = 1/36 ;
_v[X][Y].i[7].f = 1/36 ;
_v[X][Y].i[8].f = 1/36 ;
_v[X][Y].i[9].f = 1/36 ;
_v[X][Y].i[1].ft = 4/9 ;
_v[X][Y].i[2].ft = 1/9 ;
_v[X][Y].i[3].ft = 1/9 ;
_v[X][Y].i[4].ft = 1/9 ;
_v[X][Y].i[5].ft = 1/9 ;
_v[X][Y].i[6].ft = 1/36 ;
_v[X][Y].i[7].ft = 1/36 ;
_v[X][Y].i[8].ft = 1/36 ;
_v[X][Y].i[9].ft = 1/36 ;
}
}
}
}
class Particle {
public var x:Number = Math.random() * 465 ;
public var y:Number = Math.random() * 465 ;
}
class ParticleManager {
private var v:Vector.<Particle> = new Vector.<Particle> ;
public function run ( L:Lattice ):void {
while ( v.length < 5000 ) {
v.push ( new Particle ) ;
}
for each ( var P1:Particle in v ) {
var x:int = ( Math.floor ( P1.x / _LS ) + _size * 9999 ) % _size ;
var y:int = ( Math.floor ( P1.y / _LS ) + _size * 9999 ) % _size ;
P1.x -= L._v[x][y].ux * 60 ;
P1.y -= L._v[x][y].uy * 60 ;
}
for each ( var P2:Particle in v ) {
Main._canvas.setPixel ( P2.x , P2.y , 0xFFFFFF ) ;
}
for ( var I:int = 0 ; I < v.length ; ++ I ) {
if ( v[I].x < 0 ) {
v.splice ( I , 1 ) ;
continue ;
}
if ( v[I].y < 0 ) {
v.splice ( I , 1 ) ;
continue ;
}
if ( 460 < v[I].x ) {
v.splice ( I , 1 ) ;
continue ;
}
if ( 460 < v[I].y ) {
v.splice ( I , 1 ) ;
continue ;
}
}
}
}
//格子の数
const _size:int = 46 ;
//格子のサイズ
const _LS:int = 10 ;
//※たぶん粘性
const _OMEGA:Number = 0.999 ;