粒子法による流体解析

termat:http://termat.sakura.ne.jp/
*
*　　MPS法（moving particle semi-implicit method、粒子法の一種）　による流体解析の習作
*
* 　参考文献：粒子法(計算力学レクチャーシリーズ) 越塚誠一著丸善出版
* url=http://pub.maruzen.co.jp/book_magazine/book_data/search/4621075225.html
*
* 　1.注意事項
* 　　　計算にとても時間がかかります。
* OS=Windows7、CPU=Core i7、 Memory=4GBの環境で、計算に、
* 問題No1=80秒、問題No2=247秒、問題No3＝335秒を要しました。
*
* 　2.操作
* 　　　1)チェックボックスで問題を選ぶと、計算を開始します。
* 　　　2)実時間で2秒分の計算を終えると、計算結果を表示します。
* 　　　3)計算が終了した段階で画面をクリックすると、計算結果を再表示します。
*
* 　3.備考
* 　　　計算条件ファイル（CSVファイル形式）は、以下の書式でデータを保存している。
* 　　　1行目-------平均粒子間距離
* 　　　2行目以降----格子タイプ,x座標,y座標
*

termat

forked:0
favorite:8
lines:771
license : see code comments
modified : 2010-03-21 03:37:26

Embed: Tweet

/*
 *  termat:http://termat.sakura.ne.jp/ 
 * 
 *　　MPS法（moving particle semi-implicit method、粒子法の一種）　による流体解析の習作
 * 
 * 　参考文献：粒子法(計算力学レクチャーシリーズ) 越塚 誠一 著 丸善出版
 *          url=http://pub.maruzen.co.jp/book_magazine/book_data/search/4621075225.html
 * 
 * 　1.注意事項
 * 　　　計算にとても時間がかかります。
 *    OS=Windows7、CPU=Core i7、 Memory=4GBの環境で、計算に、
 *    問題No1=80秒、問題No2=247秒、 問題No3＝335秒を要しました。
 *  
 * 　2.操作
 * 　　　1)チェックボックスで問題を選ぶと、計算を開始します。
 * 　　　2)実時間で2秒分の計算を終えると、計算結果を表示します。
 * 　　　3)計算が終了した段階で画面をクリックすると、計算結果を再表示します。
 * 
 * 　3.備考
 * 　　　計算条件ファイル（CSVファイル形式）は、以下の書式でデータを保存している。
 * 　　　1行目-------平均粒子間距離
 * 　　　2行目以降----格子タイプ,x座標,y座標
 * 
 */
package
{
	import flash.display.Bitmap;
	import flash.display.BitmapData;
	import flash.display.MovieClip;
	import flash.display.Sprite;
	import flash.events.Event;
	import flash.events.MouseEvent;
	import flash.net.URLLoader;
	import flash.net.URLRequest;
	import flash.system.Security;
	import flash.text.TextField;
	import flash.text.TextFormat;
	import com.bit101.components.*;
	
	[SWF(width = "480", height = "480", backgroundColor = "0x000000", fps = "60")] 
	public class MPS3 extends Sprite {
		private var loader:URLLoader;
		private var cont:MPSController;
		private var scale:Number;
		private var colors:Array = [0x0000ff, 0x555555, 0xaaaaaa];
		private var radius:Number;
		private var time:TextField;
		private var bitmap:BitmapData;
		private var mc:MovieClip;
		private var mc2:MovieClip;
		private var ball:Vector.<Metaball>;
		private var rec:Vector.<Array>;
		private var times:Vector.<Number>;
		private var progress:ProgressBar;
		private var endTime:Number = 2.0;
		private var isThinking:Boolean = true;
		private var counter:int = 0;
		private var nextTime:Number = 0.02;
		private var isRendering:Boolean = true;

public function MPS3() {
			mc = new MovieClip();
			bitmap = new BitmapData(stage.stageWidth, stage.stageHeight, false, 0x000000);
			mc.addChild(new Bitmap(bitmap));
			addChild(mc);
			mc2 = new MovieClip();
			addChild(mc2);
			cont = new MPSController();
			scale = 480.0;
			time = new TextField();
			var tf:TextFormat = new TextFormat();
			tf.color = 0xffffff;
			time.defaultTextFormat = tf;
			time.text = "time=0.0";
			time.x = 10;	
			time.y = 10;
			mc2.addChild(time);
			progress=new ProgressBar();
			progress.x = 140;
			progress.y = 60;
			addChild(progress);
			rec = new Vector.<Array>();
			times = new Vector.<Number>();
			var prob1 : RadioButton = new RadioButton(this, 420, 5, "No.1", true, function(e:Event):void { init("prob01.csv"); } );
			var prob2 : RadioButton = new RadioButton(this, 420, 20, "No.2",false, function(e:Event):void { init("prob02.csv"); } );
			var prob3 : RadioButton = new RadioButton(this, 420, 35, "No.3",false, function(e:Event):void { init("prob03.csv"); } );
			addChild(prob1); addChild(prob2); addChild(prob3);
			addEventListener(Event.ENTER_FRAME, update);
			stage.addEventListener(MouseEvent.MOUSE_DOWN, onDown);
			init("prob01.csv");
		}
		
		private function init(file:String):void {
			while (cont.array.length > 0) cont.array.pop();
			while (rec.length > 0) rec.pop();
			while (times.length > 0) times.pop();
			isRendering = true
			load("http://termat.sakura.ne.jp/air/"+file, "http://termat.sakura.ne.jp/crossdomain.xml");
		}
		
		private function load(url:String, sec:String = null):void {
			loader = new URLLoader();
			if (sec != null)Security.loadPolicyFile(sec);
			loader.addEventListener(Event.COMPLETE, onCompleted);
			loader.load(new URLRequest(url));
		}
	
		private function onDown(e:MouseEvent):void {
			counter = 0.0;
			isRendering = !isRendering;
		}
		
		private function onCompleted(e:Event):void {
			var str:String = loader.data as String;
			str = (str.split("\r\n")).join("\n");
			str = (str.split("\r")).join("\n");
			var data:Array = new Array();
			var tmp:Array = str.split("\n");
			var l:Number = tmp.length;
			for (var i:Number = 0; i < l; i++) data.push(tmp[i].split(","));
			if (data[data.length - 1].length == 1 && data[data.length - 1][0] == "") { data.pop(); }
			cont.aveParticleDis = parseFloat(data[0][0]);
			radius = cont.aveParticleDis / 2.0;
			ball = new Vector.<Metaball>();
			var cs:Vector.<Array> = createPallet();
			for (i = 1; i < data.length; i++) {
				tmp = data[i];
				var p:Particle = new Particle(parseInt(tmp[0]), parseFloat(tmp[1]), parseFloat(tmp[2]), 0, 0, 0, 0);
				cont.array.push(p);
				var m:Metaball = new Metaball(p, cs, radius, scale);
				m.init();
				ball.push(m);
			}
			cont.init();
			isThinking = true;
			counter = 0;
			nextTime = 0.02;
			draw();
		}
		
		private function update(e:Event):void {
			if (isThinking) {
				if (cont.totalTime > nextTime) record();
				cont.solve();
				if (cont.totalTime > endTime) {
					record();
					isThinking = false;
					progress.clear();
				}else {
					progress.update(cont.totalTime / endTime);
				}
			}else if (rec.length > counter) {
				time.text = "time=" + times[counter].toString();
				time.width = time.textWidth + 10;
				var i:int = 0, it:int = 0;
				var pos:Array = rec[counter];
				for each(var p:Particle in cont.array) {
					p.x = pos[i++];	p.y = pos[i++];
				}
				draw();
				counter++;
			}
		}
		
		private function draw():void {
			var it:int = 0;
			mc2.graphics.clear();
			bitmap.fillRect(bitmap.rect, 0x000000);
			for each(var p:Particle in cont.array) {
				if (p.type == 0 && isRendering) {
					ball[it].draw(bitmap);
				}else {
					mc2.graphics.beginFill(colors[p.type]);
					mc2.graphics.drawCircle(p.x * scale + 120, -(p.y * scale-400), radius * scale);
					mc2.graphics.endFill();
				}
				it++;
			}
		}
		
		private function record():void {
			var tmp:Array = new Array();
			for each(var p:Particle in cont.array) {
				tmp.push(p.x);	tmp.push(p.y);
			}
			rec.push(tmp);
			times.push(cont.totalTime);
			nextTime += 0.02;
		}
		
		private function createPallet():Vector.<Array>{
			var ret:Vector.<Array> = new Vector.<Array>();
			var r:int, g:int, b:int;
			for (var i:int=0;i<256;i++){
				r = g = b = 0;
				if (i >= 0) b = Math.min(255, 4 * i * 2);
				if (i >= 2) g = Math.min(255, 4 * (i / 8));
				if (i >= 4) r = Math.min(255, 4 * (i / 16));
				ret[i]=[r,g,b];
			}
			return ret;
		}
	}
}

import flash.display.MovieClip;
import flash.text.TextField;
import flash.text.TextFormat;
import flash.display.BitmapData;
class ProgressBar extends MovieClip {
	private var mes:TextField;
	
	public function ProgressBar() {
		super();
		mes = new TextField();
		var tf:TextFormat = new TextFormat();
		tf.color = 0xffffff;
		tf.size = 16;
		mes.defaultTextFormat = tf;
		mes.x = 90;
		addChild(mes);
	}
	
	public function update(r:Number):void {
		graphics.clear();
		graphics.lineStyle(1, 0x0000ff);
		graphics.drawRect(0, 0, 200, 30);
		graphics.beginFill(0x0000ff);
		graphics.drawRect(0, 0, 200 * r, 30);
		graphics.endFill();
		r = Math.round(r * 100.0);
		mes.text = r.toString() + "%";
		mes.width = mes.textWidth + 5;
	}
	
	public function clear():void {
		graphics.clear();
		mes.text = "";
	}
}

class Metaball {
	public var colors:Vector.<Array> ;
	private var p:Particle;
	private var pixels:Vector.<Pixel>;
	private var RADIUS:int;
	private var NUM_PIXELS:int;
	private var radius:Number;
	private var scale:Number;
		
	public function Metaball(p:Particle,c:Vector.<Array>,radius:Number,scale:Number){
		this.p = p;
		this.radius = radius;
		this.scale = scale;
		RADIUS = Math.floor(radius * scale * 3) as int;
		NUM_PIXELS=(2 * RADIUS) * (2 * RADIUS);
		pixels = new Vector.<Pixel>();
		for (var i:int = 0; i < NUM_PIXELS; i++) pixels[i] = new Pixel();
		colors=c;
	}
		
	public function init():void {
		var no:int;
		var c:int = 0;
		for (var i:int = -RADIUS; i < RADIUS; i++) {
			for (var j:int = -RADIUS; j < RADIUS; j++) {
				var z:Number = RADIUS * RADIUS - i * i - j * j;
				if (z < 0) {
					no = 0;
				}else{
					z = Math.sqrt(z);
					var t:Number = z / RADIUS;
					no = (int) (255 * (t * t * t * t));
					if (no > 255)no = 255;
					if (no < 0)no = 0;
				}
				pixels[c].dx = i;
				pixels[c].dy = j;
				pixels[c].no = no;
				c++;
			}
		}
	}
		
	public function draw(bitmap:BitmapData):void {
		var x:int = Math.round(p.x * scale + 120) as int;
		var y:int = Math.round( -(p.y * scale-400)) as int;
		for (var i:int = 0; i < NUM_PIXELS; i++) {
			var sx:int = x + pixels[i].dx;
			if (sx < 0 || sx > bitmap.width-1)continue;
			var sy:int = y + pixels[i].dy;
			if (sy < 0 || sy > bitmap.height-1)continue;
			var no:int = pixels[i].no;
			var color:Array = colors[no];
			var c:uint = bitmap.getPixel(sx, sy);
			var p:Array = getRgb(c);
			for (var j:int = 0; j < 3; j++) {
				p[j] += color[j];
				if (p[j] > 255)p[j] = 255;
			}
			bitmap.setPixel(sx, sy, get16(p[0], p[1], p[2]));
		}
	}
	
	private function get16(rr:int,gg:int,bb:int):uint {
		return (rr << 16) + (gg << 8) + bb;
	}
	
	private function getRgb(c:uint):Array {
		var r:int = (c & 0xff0000) >> 16;
		var g:int = (c & 0xff00) >> 8;
		var b:int = (c & 0xff);
		return [r, g, b];
	}
}

class Pixel {
	public var dx:int;
	public var dy:int;
	public var no:int;
}

class Particle {
	private static const WATER:int = 0;
	private static const WALL:int = 2;
	private static const PWALL:int = 1;
	private static const RHO:Array = [1000.0, 1000.0, 1000.0];	//密度
	private static const ALPHA:Array = [0.0000001, 0.0000001, 0.0000001];//比圧縮率	
	private static var DIM:Number = 2; 
	private static var neighSize:int = 100 ; 
	public var press:Number=0.0; 
	public var dencity:Number=0.0; 
	public var bcon:int=0 , check:int=0; 
	public var x:Number=0.0 , y:Number=0.0; 
	public var vx:Number=0.0 , vy:Number=0.0; 
	private var dx:Number=0.0 , dy:Number=0.0; 
	public var poiss:Array , ic:Array;
	public  var dp:Number=0.0, source:Number=0.0, pmin:Number=0.0; 
	public var type:int=0 ;

public function Particle( type:int , x:Number , y:Number , vx:Number , vy:Number , p:Number , d:Number ) { 
		this.x = x ;	this.y = y ;
		this.type = type ;
		press = p ;	dencity = d ;
		dp = 0.0 ;
		this.vx = vx ;	this.vy = vy ;
		dx = 0.0 ;	dy = 0.0 ;
		poiss = createArray(neighSize + 1);
		ic = createArray(neighSize + 1);
	}
	
	private function createArray(n:int):Array {
		var ret:Array = new Array();
		for (var i:int = 0; i < n; i++) ret[i] = 0;
		return ret;
	}
	
	public function distanceSq( p:Particle ):Number { 
		var xx:Number = x - p.x;
		var yy:Number = y - p.y;
		return ( xx * xx + yy * yy) ;
	}
	
	public function getVerocity():Number { 
		return Math.sqrt(vx * vx + vy * vy) ;
	}
	
	public function calcBuoyancy( g:Number , dt:Number ):void { 
		vy -= g * dt ;
	}
	
	public function calcPNDensity( re:Number , array:Array , neigh:Array , num:int ):void {
		dencity = calcDensity(re , array , neigh , num) ;
	}
	
	public function calcDensity( re:Number , array:Array , neigh:Array , num:int ):Number {
		var ret:Number = 0 ;
		for( var i:int = 0 ; i < num ; i++ ) {
			var r:Number = Math.sqrt(distanceSq(array[neigh[i]]));
			var w:Number = getWeight(r , re) ;
			ret += w ;
		}
		return ret ;
	}
	
	public static function getWeight( r:Number , re:Number ):Number { 
		var rr:Number = r / re ;
		if( rr < 1.0 ) {
			return 1.0 / rr - 1.0 ;
		}
		else {
			return 0.0 ;
		}
	}
	
	public function calcConvection( dt:Number ):void { 
		if( isWall() || isPWall() ) return;
		x += vx * dt ;
		y += vy * dt ;
	}
	
	public function calcConvectionDv( dt:Number ):void { 
		if( isWall() || isPWall() ) return;
		x += dx * dt ;
		y += dy * dt ;
	}
	
	public function isWall():Boolean { 	return (type == WALL);}
	public function isPWall():Boolean { return (type == PWALL);}
	public function isWater():Boolean { return (type == WATER);}
	
	public function calcPressureGradientTerm( dt:Number , re:Number , n0p:Number , array:Array , neigh:Array , num:int ):void {
		for( var i:int = 0 ; i < num ; i++ ) {
			var p:Particle = array[neigh[i]] ;
			if( p.isWall()) continue ;
			var rr:Number = Math.sqrt(distanceSq(p));
			var ddv:Number = dt * ( p.dp - pmin ) / rr ;
			ddv = ddv * getWeight(rr , re) / RHO[type];
			ddv = ddv / n0p * DIM;
			dx += ddv * (x - p.x) / rr ;
			dy += ddv * (y - p.y) / rr ;
		}
	}
	
	public function calcCollision( dt:Number , r2lim:Number , col_rat:Number , id:int , array:Array , neigh:Array , num:int ):void {
		var m1:Number = RHO[type] ;
		var gx:Number , gy:Number , gz:Number ;
		for( var i:int = 0 ; i < num ; i++ ) {
			if( neigh[i] <= id ) continue ;
			var p:Particle = array[neigh[i]] ;
			var xx:Number = p.x - x ;
			var yy:Number = p.y - y ;
			var rr2:Number = distanceSq(p) ;
			if( rr2 < r2lim ) {
				var rr:Number = Math.sqrt(rr2) ;
				var m2:Number = RHO[p.type] ;
				var mm:Number = m1 + m2 ;
				gx = ( m1 * vx + m2 * p.vx ) / mm ;
				gy = ( m1 * vy + m2 * p.vy ) / mm ;
				gx = m1 * ( vx - gx ) ;
				gy = m1 * ( vy - gy ) ;
				var vabs:Number = ( gx * xx + gy * yy) / rr ;
				if( vabs < 0.0 ) continue ;
				var vrat:Number = 1.0 + col_rat ;
				gx = vrat * vabs * xx / rr ;
				gy = vrat * vabs * yy / rr ;
				if( !isWall() && !isPWall() ) {
					vx -= gx / m1 ;
					vy -= gy / m1 ;
					x -= dt * gx / m1 ;
					y -= dt * gy / m1 ;
				}
				if( !p.isWall() && !p.isPWall() ) {
					p.vx += gx / m2 ;
					p.vy += gy / m2 ;
					p.x += dt * gx / m2 ;
					p.y += dt * gy / m2 ;
				}
			}
		}
	}
	
	public function setMatrix( dt:Number , lambda:Number , re:Number , n0:Number , array:Array , iccg:Array , num:int ):void {
		poiss[0] = 0 ;
		for( var i:int = 0 ; i < num ; i++ ) {
			var p:Particle = array[iccg[i]] ;
			if( p.bcon == -1 ) {
				poiss[i + 1] = 0.0 ;
			}else {
				var rr:Number = Math.sqrt(distanceSq(p));
				var val:Number = 2.0 * DIM / lambda * getWeight(rr , re) / n0 / ((RHO[type]+RHO[p.type])/2.0);
				poiss[i + 1] = -val ;
				poiss[0] += val ;
			}
		}
		poiss[0] += ALPHA[type] / dt / dt ;
	}
	
	public function incom_decomp( id:int , array:Array , iccg:Array , num:Array ):void {
		if( bcon != 0 ) return  ;
		var sum:Number = poiss[0] ;
		for( var i:int = 0 ; i < num[id][1] ; i++ ) {
			if( iccg[id][i] > id ) continue ;
			var p:Particle = array[iccg[id][i]] ;
			if( p.bcon != 0 ) continue ;
			sum = sum - ic[i + 1] * ic[i + 1] ;
		}
		ic[0] = Math.sqrt(sum) ;
		for(i = 0 ; i < num[id][1] ; i++ ) {
			if( iccg[id][i] < id ) continue ;
			p = array[iccg[id][i]] ;
			if( p.bcon != 0 ) continue ;
			sum = poiss[i + 1] ;
			var id2:int = iccg[id][i] ;
			for( var mj:int = 0 ; mj < num[id2][1] ; mj++ ) {
				if( iccg[id2][mj] >= id ) continue ;
				var q:Particle = array[iccg[id2][mj]] ;
				if( q.bcon != 0 ) continue ;
				for( var mi:int = 0 ; mi < num[id][1] ; mi++ ) {
					if( iccg[id2][mj] == iccg[id][mi] ) {
						sum = sum - ic[mi + 1] * p.ic[mj + 1] ;
						break ;
					}
				}
			}
			ic[i + 1] = sum / ic[0] ;
			for(mj = 0 ; mj < num[id2][1] ; mj++ ) {
				if( id == iccg[id2][mj] ) {
					p.ic[mj + 1] = ic[i + 1] ;
					break ;
				}
			}
		}
	}
	
	public function setBcon( nop:Number , cutoff:Number ):void { 
		if(isWall() ) {
			bcon = -1 ;	check = -1 ;
		}else if( dencity / nop < cutoff ) {
			bcon = 1 ;	check = 1 ;
		}else {
			bcon = 0 ;	check = 0 ;
		}
	}
	
	public function checkBcon( array:Array , neigh:Array , num:int ):void {
		for( var i:int = 0 ; i < num ; i++ ) {
			var p:Particle = array[neigh[i]] ;
			if( p.check == 0 ) p.check = 1 ;
		}
		check = 2 ;
	}
	
	public function setSource( dt:Number , nop:Number ):void { 
		source = 1.0 / dt / dt * ( dencity - nop ) / nop ;
	}
	
	public function coeffCheck( diag:Number ):void { 
		if( check == 0 )poiss[0] *= diag ;
	}

public function calcDv( px:Particle , dt:Number , rep:Number , nOp:Number ):void { 
		var rr:Number = Math.sqrt(distanceSq(px));
		var ddv:Number = dt * ( px.dp - pmin ) / rr ;
		ddv = ddv * getWeight(rr , rep) / RHO[type];
		ddv = ddv / nOp * DIM;
		dx += ddv * ( x - px.x ) / rr ;
		dy += ddv * ( y - px.y ) / rr ;
	}
	
	public function updatePDv():void { 
		press += dp ;
		vx += dx ;	vy += dy ;
		dp = 0 ;
		dx = 0.0 ;	dy = 0.0 ;
	}
}
class MPSController {
	public var aveParticleDis:Number;	//初期配置における近接粒子間距離
	public const GRAV:Number=9.8;		//重力加速度
	public const DT_MAX:Number=0.01;	//時間増分上限
	public const CR_MAX:Number=0.2;		//クーラン数上限
	public const EPSP:Number=0.000000001;	//圧力計算における収束条件
	public const IMPX:int=50;		//圧力計算における反復数の上限
	public const IMNP:int=10;		//圧力計算における反復数の下限
	public const RE_PND:Number=2.1;	//粒子数密度の計算に用いる重み関数の半径
	public const RE_ICCG:Number=4.0;	//圧力のラプラシアン計算に用いる重み関数の半径
	public const MIN_LEN:Number=0.5;	//粒子間の再接近距離
	public const DIRRICHLET:Number=0.97;	//自由表面の条件を与えるパラメータ
	public const COL_RATE:Number=0.2;	//接近した粒子の反発率
	public const LAM_RATE:Number=0.2;	//圧力計算の緩和係数
	public const INCREASE_DT:Number = 1.1 ; 
	
	public var n0p:Number=0.0;		//repに対応する粒子数密度
	public var n0iccg:Number = 0.0;	//reiccgに対応する粒子数密度
	public var n0pperticle:int=0;		//check
	public var array:Array;
	public var totalTime:Number ; 
	private var rep:Number ,rep2:Number, reiccg:Number , reiccg2:Number ; 
	private var r2lim:Number , ra:Number , lambda:Number ; 
	private var dt:Number ; 
	private var DIAG:Number = 2.0;
	private var maxNeigh:int = 100; 
	private var neigh:Array, iccg:Array, num:Array ; 
	private var r:Array, q:Array;
	
	public function MPSController(){ 
		array = new Array() ;
	}
	
	public function init():void { 
		if ( array.length == 0 ) return ;
		num = new Array();
		neigh = new Array();
		iccg = new Array();
		r = new Array();
		q = new Array();
		for (var i:int = 0; i < array.length; i++) {
			num[i] = [0, 0];
			var t0:Array = [];
			var t1:Array = [];
			for (var j:int = 0; j < 100; j++) {
				t0[j] = 0.0;	t1[j] = 0.0;
			}
			neigh[i] = t0;	iccg[i] = t1;
			r[i] = 0.0;	q[i] = 0.0;
		}
		totalTime = 0 ;
		rep = aveParticleDis * RE_PND ;
		reiccg = aveParticleDis * RE_ICCG ;
		reiccg2 = reiccg * reiccg ;
		rep2 = rep * rep;
		if( n0p == 0.0 ) {
			setNeighbor(rep , reiccg) ;
			var p:Particle = array[n0pperticle] ;
			p.calcPNDensity(rep , array , neigh[n0pperticle] , num[n0pperticle][0]) ;
			n0p = p.dencity;
			n0iccg = p.calcDensity(reiccg , array , iccg[n0pperticle] , num[n0pperticle][1]) ;
		}
		r2lim = Math.pow(aveParticleDis , 2) * Math.pow(MIN_LEN , 2) ;
		ra = aveParticleDis / 1.7724539 ;
		lambda = LAM_RATE * ( reiccg2 * reiccg2 / 12.0 - reiccg * ra * ra * ra / 3.0 + ra * ra * ra * ra / 4.0 ) / ( reiccg2 / 2.0 - reiccg * ra + ra * ra / 2.0 ) ;
		dt = DT_MAX;
		setNeighbor(rep , reiccg);
		for (i = 0 ; i < array.length ; i++) array[i].calcPNDensity(rep , array , neigh[i] , num[i][0]);
	}
	
	public function solve():void {
		dt = calcDt(CR_MAX , aveParticleDis , DT_MAX , dt) ;
		totalTime += dt ;
		for( var i:int = 0 ; i < array.length ; i++)array[i].calcBuoyancy(GRAV , dt);
		for (i = 0 ; i < array.length ; i++) array[i].calcConvection(dt);
		for (i = 0 ; i < array.length ; i++) array[i].calcCollision(dt, r2lim, COL_RATE, i, array, neigh[i], num[i][0]);
		setNeighbor(rep , reiccg) ;
		for (i = 0 ; i < array.length ; i++ ) array[i].calcPNDensity(rep , array , neigh[i] , num[i][0]);
		for(i = 0 ; i < array.length ; i++ )array[i].setBcon(n0p , DIRRICHLET);
		setSource(n0p);
		checkBcon();
		setMatrix();
		for(i = 0 ; i < array.length ; i++ )array[i].coeffCheck(DIAG) ;
		for(i = 0 ; i < array.length ; i++ )array[i].incom_decomp(i , array , iccg , num) ;
		pcg_solver(EPSP, IMNP, IMPX, dt);
		for(i = 0 ; i < array.length ; i++ ) if( array[i].dp < 0.0 ) array[i].dp = 0.0 ;
		set_pmin() ;
		rev_pressgrad(dt , rep , n0p) ;
		for(i = 0 ; i < array.length ; i++ )array[i].calcConvectionDv(dt) ;
		for(i = 0 ; i < array.length ; i++ )array[i].updatePDv() ;
		setNeighbor(rep , reiccg) ;
	}
	
	private function setNeighbor( rep:Number , reiccg:Number ):void { 
		for( var i:int = 0 ; i < array.length ; i++ ) {
			num[i][0] = 0;	num[i][1] = 0;
		}
		for(i = 0 ; i < array.length ; i++ ) {
			for( var j:int = i + 1 ; j < array.length ; j++ ) {
				if( array[i].isWall() && array[j].isWall() ) continue ;
				var d:Number = array[i].distanceSq(array[j]) ;
				if( d <= rep2 ) {
					neigh[i][num[i][0]++] = j ;
					neigh[j][num[j][0]++] = i ;
				}
				if( d <= reiccg2 ) {
					iccg[i][num[i][1]++] = j ;
					iccg[j][num[j][1]++] = i ;
				}
			}
		}
	}
	
	private function setMatrix():void {
		for(var i:int = 0 ; i < array.length ; i++ ) {
			if( array[i].bcon != 0 ) continue ;
			array[i].setMatrix(dt , lambda , reiccg , n0iccg , array , iccg[i] , num[i][1]) ;
		}
	}
	
	private function calcDt( cr:Number , aveParticleDis:Number , dtmax:Number , dtold:Number ):Number { 
		var vmax:Number = 0.0 ;
		for( var i:int = 0 ; i < array.length ; i++ ) {
			if( array[i].isWall() || array[i].isPWall() ) continue ;
			var vv:Number = array[i].getVerocity() ;
			if( vv > vmax ) vmax = vv ;
		}
		var dtlimit:Number = dtold * INCREASE_DT ;
		if( vmax == 0 ) return dtlimit ;
		var dt:Number = cr * aveParticleDis / vmax ;
		if( dt > dtlimit ) dt = dtlimit ;
		if( dt > dtmax ) dt = dtmax ;
		return dt ;
	}
	
	private function checkBcon():void { 
		var count:int,i:int ;
		do{
			count = 0 ;
			for(i = 0 ; i < array.length ; i++ ) {
				if( array[i].check == 1 ) {
					array[i].checkBcon(array , neigh[i] , num[i][0]) ;
					count++ ;
				}
			}
		}while( count != 0 ) ;
	}
	
	private function setSource(nop:Number ):void { 
		var num_dirichlet:int = 0 ;
		for( var i:int = 0 ; i < array.length ; i++ ) {
			if( array[i].isWall() ) continue ;
			if( array[i].bcon == 0 ) {
				array[i].setSource(dt , nop) ;
			}else if( array[i].bcon == 1 ) {
				array[i].source = 0.0;
				if( array[i].isWater()) num_dirichlet++ ;
			}
		}
	}
	
	private function pcg_solver( eps:Number , imin:Number , imax:Number , dt:Number ):Boolean { 
		for( var i:int = 0 ; i < array.length ; i++ ) {
			array[i].dp = 0.0 ;
			r[i] = 0.0;	q[i] = 0.0;
		}
		var s:Array = mul_matrix_vector() ;
		for(i = 0 ; i < array.length ; i++ ) r[i] = array[i].source - s[i] ;
		solver_ll(q , r) ;
		for(i = 0 ; i < array.length ; i++ )array[i].press = q[i];
		var rqo:Number = mul_vector(r , q) ;
		for ( var k:int = 0 ; k < imax ; k++ ) {
			mul_matrix_vector_press(s) ;
			var ps:Number = mul_vector_press(s) ;
			var aa:Number = rqo / ps ;
			for(i = 0 ; i < array.length ; i++ ) {
				var p0:Particle = array[i] ;
				p0.dp = p0.dp + aa * p0.press ;
			}
			for(i = 0 ; i < r.length ; i++ ) r[i] = r[i] - aa * s[i] ;
			solver_ll(q , r) ;
			var rqn:Number = mul_vector(r , q) ;
			var bb:Number = rqn / rqo ;
			rqo = rqn ;
			for (i = 0 ; i < array.length ; i++ ) array[i].press = (q[i] + bb * array[i].press);
			if( checkconv(r , eps , k , imin , imax , dt) ) return true ;	
		}
		return false ;
	}
	
	private function mul_matrix_vector():Array {
		var y:Array = new Array();
		for ( var i:int = 0 ; i < array.length ; i++ ) {
			y[i] = 0.0;
			if( array[i].bcon != 0 ) continue ;
			y[i] = array[i].poiss[0] * array[i].dp ;
			for( var j:int = 0 ; j < num[i][1] ; j++ ) {
				var p:Particle = array[iccg[i][j]] ;
				if( p.bcon != -1 ) continue ;
				y[i] = y[i] + p.poiss[j + 1] * p.dp ;
			}
		}
		return y ;
	}
	
	private function solver_ll( y:Array , xx:Array ):void {
		var x:Array = new Array();
		for( var i:int = 0 ; i < xx.length ; i++ ) x[i] = xx[i] ;
		for(i = 0 ; i < array.length ; i++ ) {
			var p:Particle = array[i] ;
			if( p.bcon != 0 ) continue ;
			for( var j:int = 0 ; j < num[i][1] ; j++ ) {
				var px:Particle = array[iccg[i][j]] ;
				if( iccg[i][j] > i || px.bcon != 0 ) continue ;
				x[i] = x[i] - p.ic[j + 1] * y[iccg[i][j]] ;
			}
			y[i] = x[i] / p.ic[0] ;
		}
		for(i = 0 ; i < y.length ; i++ ) x[i] = y[i] ;
		for(i = array.length - 1 ; i >= 0 ; i-- ) {
			p = array[i] ;
			if( p.bcon != 0 ) continue ;
			for(j = 0 ; j < num[i][1] ; j++ ) {
				px= array[iccg[i][j]] ;
				if( iccg[i][j] < i || px.bcon != 0 ) continue ;
				x[i] = x[i] - p.ic[j + 1] * y[iccg[i][j]] ;
			}
			y[i] = x[i] / p.ic[0] ;
		}
	}
	
	private function mul_vector( x1:Array , x2:Array ):Number {
		var ans:Number = 0.0 ;
		for( var i:int = 0 ; i < array.length ; i++ )if( array[i].bcon == 0 ) ans += x1[i] * x2[i] ;
		return ans ;
	}
	
	private function mul_matrix_vector_press( y:Array ):void {
		for ( var i:int = 0 ; i < array.length ; i++ ) {
			var p:Particle = array[i] ;
			if( p.bcon != 0 ) continue ;
			y[i] = p.poiss[0] * p.press ;
			for( var j:int = 0 ; j < num[i][1] ; j++ ) {
				var px:Particle = array[iccg[i][j]] ;
				if( px.bcon == -1 ) continue ;
				y[i] = y[i] + p.poiss[j + 1] * px.press;
			}
		}
	}
	
	private function mul_vector_press( x1:Array ):Number {
		var ans:Number = 0.0 ;
		for ( var i:int = 0 ; i < array.length ; i++ ) if( array[i].bcon == 0 ) ans += x1[i] * array[i].press;
		return ans ;
	}

private function checkconv( r:Array , eps:Number , k:Number , imin:Number , imax:Number , dt:Number ):Boolean {
		var checker:int = 0 ;
		var err1sum:Number = 0.0 ;
		for( var i:int = 0 ; i < array.length ; i++ ) {
			if( array[i].bcon != 0 ) continue ;
			var err1:Number = 0.0 ;
			var rr:Number = r[i] ;
			if ( rr > 0 ) { err1 = rr * dt * dt ; } else { err1 = -rr * dt * dt ; }
			if ( err1 > eps ) { checker++; }
			err1sum += err1 ;
		}
		if( checker == 0 && k >= imin ) {
			return true ;
		}else {
			return false ;
		}
	}
	
	private function set_pmin():void { 
		for( var i:int = 0 ; i < array.length ; i++ ) {
			var p:Particle = array[i] ;
			if( p.bcon == -1 || p.isWall() ) continue ;
			p.pmin = p.dp ;
			for( var j:int = 0 ; j < num[i][1] ; j++ ) {
				var px:Particle = array[iccg[i][j]] ;
				if( px.isWall() ) continue ;
				if( px.dp < p.pmin ) p.pmin = px.dp ;
			}
		}
	}
	
	private function rev_pressgrad( dt:Number , rep:Number , n0p:Number ):void { 
		for( var i:int = 0 ; i < array.length ; i++ ) {
			var p:Particle = array[i] ;
			if( p.bcon == -1 || p.isPWall() ) continue ;
			for( var j:int = 0 ; j < num[i][1] ; j++ ) {
				var px:Particle = array[iccg[i][j]] ;
				if( px.isWall() ) continue ;
				p.calcDv(px , dt , rep , n0p) ;
			}
		}
	}
}

Code Fullscreen

Preview Fullscreen

bradsedito

shimo3651

kidaipu

tkinjo

flasher_bbb

fabure

demouth

Mae_ITR