今回は制御工学ではなく有限要素法に関係する内容です。有限要素法についても全く素人なのですが、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] )

結果 f:id:Otepipi:20190415222048p:plain

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] )

結果 f:id:Otepipi:20190415222207p:plain

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] )

結果 f:id:Otepipi:20190415222326p:plain

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] )

結果 f:id:Otepipi:20190415222521p:plain

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] )

結果 f:id:Otepipi:20190415222657p:plain

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] )

結果 f:id:Otepipi:20190415222950p:plain

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] )

結果 f:id:Otepipi:20190415223122p:plain

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] )

結果 f:id:Otepipi:20190415223329p:plain

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)

結果 f:id:Otepipi:20190415223503p:plain

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)

結果 f:id:Otepipi:20190415223628p:plain