Octaveでメッシュを作成する DistMesh
今回は制御工学ではなく有限要素法に関係する内容です。 有限要素法についても全く素人なのですが、Octaveでメッシュを作成するプログラムが公開されていましたので、それのチュートリアルを試してみたいと思います。
DistMesh
DistMeshはMATLABとGNU Octaveで動作する自動メッシュ生成プログラムです。2Dのジオメトリでは三角形、3Dでは四面体のメッシュを生成します。
example
例題を行う準備として、上記githubに公開されているdistmesh.m
のコピーを、Octaveのカレントディレクトリに作成します。その後、githubのReadMeに記載してある例を次々に実施してみます。
Example 1: Uniform mesh on unit circle
fd = @(p) sqrt(sum(p.^2,2)) - 1; fh = @(p) ones(size(p,1),1); [p,t] = distmesh( fd, fh, 0.2, [-1,-1;1,1] ); patch( 'vertices', p, 'faces', t, 'facecolor', [.9, .9, .9] )
結果
Example 2: Uniform mesh on ellipse
fd = @(p) p(:,1).^2/2^2 + p(:,2).^2/1^2 - 1; fh = @(p) ones(size(p,1),1); [p,t] = distmesh( fd, fh, 0.2, [-2,-1;2,1] ); patch( 'vertices', p, 'faces', t, 'facecolor', [.9, .9, .9] )
結果
Example 3: Uniform mesh on unit square
fd = @(p) -min(min(min(1+p(:,2),1-p(:,2)),1+p(:,1)),1-p(:,1)); fh = @(p) ones(size(p,1),1); [p,t] = distmesh( fd, fh, 0.2, [-1,-1;1,1], [-1,-1;-1,1;1,-1;1,1] ); patch( 'vertices', p, 'faces', t, 'facecolor', [.9, .9, .9] )
結果
Example 4: Uniform mesh on complex polygon
pv = [-0.4 -0.5;0.4 -0.2;0.4 -0.7;1.5 -0.4;0.9 0.1; 1.6 0.8;0.5 0.5;0.2 1;0.1 0.4;-0.7 0.7;-0.4 -0.5]; fd = { 'l_dpolygon', [], pv }; fh = @(p) ones(size(p,1),1); [p,t] = distmesh( fd, fh, 0.1, [-1,-1; 2,1], pv ); patch( 'vertices', p, 'faces', t, 'facecolor', [.9, .9, .9] )
結果
Example 5: Rectangle with circular hole, refined at circle boundary
drectangle = @(p,x1,x2,y1,y2) -min(min(min(-y1+p(:,2),y2-p(:,2)),-x1+p(:,1)),x2-p(:,1)); fd = @(p) max( drectangle(p,-1,1,-1,1), -(sqrt(sum(p.^2,2))-0.5) ); fh = @(p) 0.05 + 0.3*(sqrt(sum(p.^2,2))-0.5); [p,t] = distmesh( fd, fh, 0.05, [-1,-1;1,1], [-1,-1;-1,1;1,-1;1,1] ); patch( 'vertices', p, 'faces', t, 'facecolor', [.9, .9, .9] )
結果
Example 6: Square, with size function point and line sources
dcircle = @(p,xc,yc,r) sqrt((p(:,1)-xc).^2+(p(:,2)-yc).^2)-r; fd = @(p) -min(min(min(p(:,2),1-p(:,2)),p(:,1)),1-p(:,1)); dpolygon = @(p,v) feval('l_dpolygon',p,v); fh = @(p) min(min(0.01+0.3*abs(dcircle(p,0,0,0)), ... 0.025+0.3*abs(dpolygon(p,[0.3,0.7;0.7,0.5;0.3,0.7]))),0.15); [p,t] = distmesh( fd, fh, 0.01, [0,0;1,1], [0,0;1,0;0,1;1,1] ); patch( 'vertices', p, 'faces', t, 'facecolor', [.9, .9, .9] )
結果
Example 7: NACA0012 airfoil
hlead = 0.01; htrail = 0.04; hmax = 2; circx = 2; circr = 4; a = 0.12/0.2*[0.2969,-0.126,-0.3516,0.2843,-0.1036]; fd = @(p) max( dcircle(p,circx,0,circr), ... -((abs(p(:,2))-polyval([a(5:-1:2),0],p(:,1))).^2-a(1)^2*p(:,1)) ); fh = @(p) min(min(hlead+0.3*dcircle(p,0,0,0),htrail+0.3*dcircle(p,1,0,0)),hmax); fixx = 1 - htrail*cumsum(1.3.^(0:4)'); fixy = a(1)*sqrt(fixx) + polyval([a(5:-1:2),0],fixx); pfix = [[circx+[-1,1,0,0]*circr; 0,0,circr*[-1,1]]'; 0,0; 1,0; fixx,fixy; fixx,-fixy]; bbox = [circx-circr,-circr; circx+circr,circr]; h0 = min([hlead,htrail,hmax]); [p,t] = distmesh( fd, fh, h0, bbox, pfix ); patch( 'vertices', p, 'faces', t, 'facecolor', [.9, .9, .9] )
結果
Example 8: Uniform mesh on unit sphere
fd = @(p) sqrt(sum(p.^2,2)) - 1; fh = @(p) ones(size(p,1),1); [p,t] = distmesh( fd, fh, 0.2, [-1,-1,-1;1,1,1] ); f = [t(:,[1:3]); t(:,[1,2,4]); t(:,[2,3,4]); t(:,[3,1,4])]; patch( 'vertices', p, 'faces', f, 'facecolor', [.9, .9, .9] )
結果
Example 9: Uniform mesh on unit cube
fd = @(p) -min(min(min(min(min(p(:,3),1-p(:,3) ),p(:,2)),1-p(:,2)),p(:,1)),1-p(:,1)); fh = @(p) ones(size(p,1),1); pfix = [-1,-1,-1;-1,1,-1;1,-1,-1;1,1,-1; -1,-1,1;-1,1,1;1,-1,1;1,1,1]; [p,t] = distmesh( fd, fh, 0.2, [-1,-1,-1;1,1,1], pfix ); f = [t(:,[1:3]); t(:,[1,2,4]); t(:,[2,3,4]); t(:,[3,1,4])]; patch( 'vertices', p, 'faces', f, 'facecolor', [.9, .9, .9] ), view(3)
結果
Example 10: Uniform mesh on cylinder
fd = @(p) -min(min(p(:,3),4-p(:,3)),1-sqrt(sum(p(:,1:2).^2,2))); fh = @(p) ones(size(p,1),1); pfix = [-1,-1,-1;-1,1,-1;1,-1,-1;1,1,-1; -1,-1,1;-1,1,1;1,-1,1;1,1,1]; [p,t] = distmesh( fd, fh, 0.5, [-1,-1,0;1,1,4], [] ); f = [t(:,[1:3]); t(:,[1,2,4]); t(:,[2,3,4]); t(:,[3,1,4])]; patch( 'vertices', p, 'faces', f, 'facecolor', [.9, .9, .9] ), view(3)
結果
今回は以上になります。 チュートリアルを実行しただけでしたが、また何かに活用できないか考えてみたいです。