/* * 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) ; } } } } 粒子法による流体解析