例題

ステップごとの基本的な例

この例は，以下のような über-canonical 問題に対するgetfemの基本的な使い方を示しています．これは Laplacian, \(-\Delta u = f\) を正方形上で，Dirichlet条件 \(u = g(x)\) をドメイン境界上で解きます．この例の m-file は demo_step_by_step.m という名前の下にあります． GetFEM ディストリビューションの interface/tests/matlab-octave/ ディレクトリにあります．

最初のステップは meshを作成する ことです．単純なジオメトリ( gf_mesh('generate', mesher_object mo, scalar h)) を参照)用に単純な構造化メッシュまたは非構造化メッシュを作成することも，外部メッシャ( gf_mesh('import', string FORMAT, string FILENAME) を参照)に依存します．この例では，ノードが \(\{x_{i=0\ldots10,j=0..10}=(i/10,j/10)\}\) である通常の cartesian mesh を考えます．

>> % creation of a simple cartesian mesh
>> m = gf_mesh('cartesian',[0:.1:1],[0:.1:1]);
m =
     id: 0
    cid: 0

m の値を調べれば，それが2つの整数を含む構造体であることに気づくでしょう．最初のものはその識別子であり，2番目のものはそのクラスID，すなわちその型の識別子です．この小さな構造は，実際のオブジェクトに対する単なる "handle" または "descriptor" であり， GetFEM メモリに格納され， Octave と MatLab データ構造では表現できません．いずれにしても， gf_workspace('stats') コマンドを使って|gf|オブジェクトを検査することはできます．

次に，頂点の番号と凸の番号が付いた メッシュを見てみましょう ．

>> % we enable vertices and convexes labels
>> gf_plot_mesh(m, 'vertices', 'on', 'convexes', 'on');

このように，メッシュは規則的であり，節点と凸の番号付けも規則的です(これは直交メッシュでは保証されていますが，自由度に同様の番号付けは望めません)．

次のステップでは mesh_fem オブジェクト を作成します．これは，メッシュをFEMの集合とリンクします．

>> % create a mesh_fem of for a field of dimension 1 (i.e. a scalar field)
>> mf = gf_mesh_fem(m,1);
>> gf_mesh_fem_set(mf,'fem',gf_fem('FEM_QK(2,2)'));

最初の命令は新しい gfMeshFem オブジェクトを構築し，2番目の引数はこのオブジェクトがスカラーフィールドの補間に使用されることを指定します(未知変数 \(u\) はスカラーフィールドです)．第2の命令は，すべての凸(各基底関数は4次の多項式であり \(k\) 次 \(P_k\) 多項式であることを覚えておいてください，ここで \(2k\) 次 \(Q_k\) 多項式です)に \(Q_2\) 有限要素法を割り当てます． \(Q_2\) は多項式有限要素法なので，参照要素でその基底関数の式を見ることができます．

>> gf_fem_get(gf_fem('FEM_QK(2,2)'), 'poly_str');
ans =
    '1 - 3*x - 3*y + 2*x^2 + 9*x*y + 2*y^2 - 6*x^2*y - 6*x*y^2 + 4*x^2*y^2'
    '4*x - 4*x^2 - 12*x*y + 12*x^2*y + 8*x*y^2 - 8*x^2*y^2'
    '-x + 2*x^2 + 3*x*y - 6*x^2*y - 2*x*y^2 + 4*x^2*y^2'
    '4*y - 12*x*y - 4*y^2 + 8*x^2*y + 12*x*y^2 - 8*x^2*y^2'
    '16*x*y - 16*x^2*y - 16*x*y^2 + 16*x^2*y^2'
    '-4*x*y + 8*x^2*y + 4*x*y^2 - 8*x^2*y^2'
    '-y + 3*x*y + 2*y^2 - 2*x^2*y - 6*x*y^2 + 4*x^2*y^2'
    '-4*x*y + 4*x^2*y + 8*x*y^2 - 8*x^2*y^2'
    'x*y - 2*x^2*y - 2*x*y^2 + 4*x^2*y^2'

MatLab の "object oriented" 機能を利用することもできます．お気付きかもしれませんが，クラス "foo" が getfem-interface によって提供される場合，それは関数 gf_foo でビルドされ，関数 gf_foo_get および gf_foo_set で操作されます．しかし， gfFoo コンストラクタでオブジェクトを作成し， get(..) および set(..) メソッドで操作することもできます．たとえば，前の手順は次のようになります．

>> gfFem('FEM_QK(2,2)');
gfFem object ID=0 dim=2, target_dim=1, nbdof=9,[EQUIV, POLY, LAGR], est.degree=4
  -> FEM_QK(2,2)
>> m=gfMesh('cartesian', [0:.1:1], [0:.1:1]);
gfMesh object ID=0 [16512 bytes], dim=2, nbpts=121, nbcvs=100
>> mf=gfMeshFem(m,1);
gfMeshFem object: ID=1 [804 bytes], qdim=1, nbdof=0,
  linked gfMesh object: dim=2, nbpts=121, nbcvs=100
>> set(mf, 'fem', gfFem('FEM_QK(2,2)'));
>> mf
gfMeshFem object: ID=1 [1316 bytes], qdim=1, nbdof=441,
  linked gfMesh object: dim=2, nbpts=121, nbcvs=100

さて， mf で数値積分を行うには， mesh_im オブジェクトを構築する 必要があります．

>> % assign the same integration method on all convexes
>> mim = gf_mesh_im(m, gf_integ('IM_EXACT_PARALLELEPIPED(2)'));

積分法を使用して各要素の各種積分を計算します．ここでは，正確な計算(非 求積公式)を実行することを選択します．これは，参照凸からのこれらの凸の幾何変換が線形(これはすべての単体に当てはまり，通常のメッシュの平行六面体にも当てはまりますが，一般的な四角形には当てはまりません)であり，選択したFEMが多項式であるために可能です．このため，すべての基底関数/基底関数の積/傾き/などを解析的に積分することが可能であり，FEMの手法や積分法にも多くの選択肢があります( ユーザドキュメント を参照)．

ただし，一般的なケースでは，厳密な積分法よりも近似積分法の方が適しています．

次に，Dirichlet条件を設定するために，領域の "boundary" を 検索する 必要があります．メッシュオブジェクトには，凸面と凸面の集合を保存する機能があります．これらの集合("regions" と呼ばれます)には，整数 #id を使用してアクセスします．

>> % detect the border of the mesh
>> border = gf_mesh_get(m,'outer faces');
>> % mark it as boundary #42
>> gf_mesh_set(m, 'region', 42, border);
>> gf_plot_mesh(m, 'regions', [42]); % the boundary edges appears in red

ここでは，メッシュの境界上にある凸の面(つまり，2つの凸によって共有されていない面)を見つけます．

備考:

インターフェースは大文字小文字を区別せず，空白はアンダースコアで置き換えることができるので， gf_mesh_get(m, 'OuTEr_faCes') を使うこともできます．

配列 border には二つの行があり，最初の行は凸の番号で，次の行は面の番号です(これは凸に対してローカルであり，面のグローバルな番号付けはありません)．次に，この面の集合を領域番号42に割り当てます．

この時点では，モデルを記述し，ソルバを実行して解を取得するだけです． "model" は gf_model (または gfModel )コンストラクタで作成されます．モデルは基本的に，全体線形システム(非線形問題の接線行列)とそれに関連するRHSを構築するオブジェクトです．代表的な修正は，考慮される問題(線形弾性，Laplacian等)に対する剛性マトリックスの挿入，拘束の集合の処理，Dirichlet条件，ソース項のRHSへの追加などです．全体正接行列とその右辺は， "model" 構造に格納されます．

簡単な方法で問題を作成してみましょう． \(u=x(x-1)y(y-1)+x^5\) とすると， \(\Delta u=2(x^2+y^2)-2(x+y)+20x^3\) ( \(Q^2\) 法を使っているので，FEMは正確な解を捕らえることができません)となります．

まず空の実数モデルから始めます．

>> % empty real model
>> md = gf_model('real');

(モデルは 'real' か 'complex' のどちらかです)．そして有限要素法 mf のシステムでは， u は未知であると宣言します．

>> % declare that "u" is an unknown of the system
>> % on the finite element method `mf`
>> gf_model_set(md, 'add fem variable', 'u', mf);

ここで，以下の問題を処理する "generic elliptic" ブリックを追加します．ここで \(-\nabla\cdot(A:\nabla u) = \ldots\) 問題をハンドルします．ここで \(A\) はスカラー場，行列場，4次テンソル場の場合があります．デフォルトでは， \(A=1\) です．これを変数 u に追加します．

>> % add generic elliptic brick on "u"
>> gf_model_set(md, 'add Laplacian brick', mim, 'u');

次に，領域の境界にDirichlet条件を追加します．

>> % add Dirichlet condition
>> Uexact = gf_mesh_fem_get(mf, 'eval', {'(x-.5).^2 + (y-.5).^2 + x/5 - y/3'});
>> gf_model_set(md, 'add initialized fem data', 'DirichletData', mf, Uexact);
>> gf_model_set(md, 'add Dirichlet condition with multipliers', mim, 'u', mf, 42, 'DirichletData');

最初の2行は，Dirichlet条件の値を表すモデルのデータを定義します．3つ目は，境界番号 42 の変数 u にDirichlet条件を追加します．Dirichlet条件はLagrange乗数を用いて課しました．他にはペナルティ法があります． gfMeshFem 引数も必要です．これはDirichlet条件 \(u=g\) が弱形式 \(\int_\Gamma u(x)v(x) = \int_\Gamma g(x)v(x)\ \forall v\) で課されるからです．ここで \(v\) はここでは mf で与えられる乗数空間で使われます．

備考:

多項式は mf で補間されます．これは mf がLagrange型の場合のみ可能です．この最初の例では，未知変数と g のようなデータに同じ gfMeshFem を使用していますが，一般的なケースでは mf はLagrangianではなく，Dirichlet条件やソース項などの記述には別の(Lagrangian) mesh_fem が使用されます．

ソース項は以下の行で追加できます．

>> % add source term
>> f = gf_mesh_fem_get(mf, 'eval', { '2(x^2+y^2)-2(x+y)+20x^3' });
>> gf_model_set(md, 'add initialized fem data', 'VolumicData', mf, f);
>> gf_model_set(md, 'add source term brick', mim, 'u', 'VolumicData');

ソルバを起動するだけです．線形システムは次の命令で組み立てられ，解かれます．

>> % solve the linear system
>> gf_model_get(md, 'solve');

これで，(解かれた線形系のような他のものと同様に)モデルに解が含まれました．抽出され， Octave か MatLab 図に表示されます．

>> % extracted solution
>> u = gf_model_get(md, 'variable', 'u');
>> % display
>> gf_plot(mf, u, 'mesh','on');

別のLaplacianによる厳密解

これは tests/matlab-octave/demo_laplacian.m の例です．

% trace on;
gf_workspace('clear all');
m = gf_mesh('cartesian',[0:.1:1],[0:.1:1]);
%m=gf_mesh('import','structured','GT="GT_QK(2,1)";SIZES=[1,1];NOISED=1;NSUBDIV=[1,1];')

% create a mesh_fem of for a field of dimension 1 (i.e. a scalar field)
mf = gf_mesh_fem(m,1);
% assign the Q2 fem to all convexes of the mesh_fem,
gf_mesh_fem_set(mf,'fem',gf_fem('FEM_QK(2,2)'));

% Integration which will be used
mim = gf_mesh_im(m, gf_integ('IM_GAUSS_PARALLELEPIPED(2,4)'));
%mim = gf_mesh_im(m, gf_integ('IM_STRUCTURED_COMPOSITE(IM_GAUSS_PARALLELEPIPED(2,5),4)'));
% detect the border of the mesh
border = gf_mesh_get(m,'outer faces');
% mark it as boundary #1
gf_mesh_set(m, 'boundary', 1, border);
gf_plot_mesh(m, 'regions', [1]); % the boundary edges appears in red
pause(1);

% interpolate the exact solution
Uexact = gf_mesh_fem_get(mf, 'eval', { 'y.*(y-1).*x.*(x-1)+x.^5' });
% its second derivative
F      = gf_mesh_fem_get(mf, 'eval', { '-(2*(x.^2+y.^2)-2*x-2*y+20*x.^3)' });


md=gf_model('real');
gf_model_set(md, 'add fem variable', 'u', mf);
gf_model_set(md, 'add Laplacian brick', mim, 'u');
gf_model_set(md, 'add initialized fem data', 'VolumicData', mf, F);
gf_model_set(md, 'add source term brick', mim, 'u', 'VolumicData');
gf_model_set(md, 'add initialized fem data', 'DirichletData', mf, Uexact);
gf_model_set(md, 'add Dirichlet condition with multipliers', mim, 'u', mf, 1, 'DirichletData');

gf_model_get(md, 'solve');
U = gf_model_get(md, 'variable', 'u');

% Version with old bricks
% b0=gf_mdbrick('generic elliptic',mim,mf);
% b1=gf_mdbrick('dirichlet', b0, 1, mf, 'penalized');
% gf_mdbrick_set(b1, 'param', 'R', mf, Uexact); 
% b2=gf_mdbrick('source term',b1);
% gf_mdbrick_set(b2, 'param', 'source_term', mf, F);
% mds=gf_mdstate(b1);
% gf_mdbrick_get(b2, 'solve', mds)
% U=gf_mdstate_get(mds, 'state');

disp(sprintf('H1 norm of error: %g', gf_compute(mf,U-Uexact,'H1 norm',mim)));

subplot(2,1,1); gf_plot(mf,U,'mesh','on','contour',.01:.01:.1); 
colorbar; title('computed solution');

subplot(2,1,2); gf_plot(mf,U-Uexact,'mesh','on'); 
colorbar;title('difference with exact solution');

線形および非線形弾性

この例では， GiD で生成されたメッシュを使用しています．オブジェクトは2次4面体でメッシュ分割されます．この例の m-file は GetFEM ディストリビューションの tests/matlab-octave ディレクトリの demo_tripod.m にあります．

disp('This demo is an adaption of the original tripod demo')
disp('which uses the new "brick" framework of getfem')
disp('The code is shorter, faster and much more powerful')
disp('You can easily switch between linear/non linear')
disp('compressible/incompressible elasticity!')

linear = 1
incompressible = 0


gf_workspace('clear all');
% import the mesh
m=gfMesh('import','gid','../meshes/tripod.GiD.msh');
mfu=gfMeshFem(m,3);     % mesh-fem supporting a 3D-vector field
mfd=gfMeshFem(m,1);     % scalar mesh_fem, for data fields.
% the mesh_im stores the integration methods for each tetrahedron
mim=gfMeshIm(m,gf_integ('IM_TETRAHEDRON(5)'));
% we choose a P2 fem for the main unknown
gf_mesh_fem_set(mfu,'fem',gf_fem('FEM_PK(3,2)'));
% the material is homogeneous, hence we use a P0 fem for the data
gf_mesh_fem_set(mfd,'fem',gf_fem('FEM_PK(3,0)'));
% display some informations about the mesh
disp(sprintf('nbcvs=%d, nbpts=%d, nbdof=%d',gf_mesh_get(m,'nbcvs'),...
             gf_mesh_get(m,'nbpts'),gf_mesh_fem_get(mfu,'nbdof')));
P=gf_mesh_get(m,'pts'); % get list of mesh points coordinates
pidtop=find(abs(P(2,:)-13)<1e-6); % find those on top of the object
pidbot=find(abs(P(2,:)+10)<1e-6); % find those on the bottom
% build the list of faces from the list of points
ftop=gf_mesh_get(m,'faces from pid',pidtop); 
fbot=gf_mesh_get(m,'faces from pid',pidbot);
% assign boundary numbers
gf_mesh_set(m,'boundary',1,ftop);
gf_mesh_set(m,'boundary',2,fbot);

E = 1e3; Nu = 0.3;
% set the Lame coefficients
lambda = E*Nu/((1+Nu)*(1-2*Nu));
mu = E/(2*(1+Nu));

% create a meshfem for the pressure field (used if incompressible ~= 0)
mfp=gfMeshFem(m); set(mfp, 'fem',gfFem('FEM_PK_DISCONTINUOUS(3,0)'));
if (linear)
  % the linearized elasticity , for small displacements
  b0 = gfMdBrick('isotropic_linearized_elasticity',mim,mfu)
  set(b0, 'param','lambda', lambda);
  set(b0, 'param','mu', mu);
  if (incompressible)
    b1 = gfMdBrick('linear incompressibility term', b0, mfp);
  else
    b1 = b0;
  end;
else
  % See also demo_nonlinear_elasticity for a better example
  if (incompressible)
    b0 = gfMdBrick('nonlinear elasticity',mim, mfu, 'Mooney Rivlin');
    b1 = gfMdBrick('nonlinear elasticity incompressibility term',b0,mfp);
    set(b0, 'param','params',[lambda;mu]);
  else
    % large deformation with a linearized material law.. not
    % a very good choice!
    b0 = gfMdBrick('nonlinear elasticity',mim, mfu, 'SaintVenant Kirchhoff');
    set(b0, 'param','params',[lambda;mu]);
    %b0 = gfMdBrick('nonlinear elasticity',mim, mfu, 'Ciarlet Geymonat');
    b1 = b0;
  end;
end

% set a vertical force on the top of the tripod
b2 = gfMdBrick('source term', b1, 1);
set(b2, 'param', 'source_term', mfd, get(mfd, 'eval', {0;-10;0}));

% attach the tripod to the ground
b3 = gfMdBrick('dirichlet', b2, 2, mfu, 'penalized');

mds=gfMdState(b3)

disp('running solve...')

t0=cputime; 

get(b3, 'solve', mds, 'noisy', 'max_iter', 1000, 'max_res', 1e-6, 'lsolver', 'superlu');
disp(sprintf('solve done in %.2f sec', cputime-t0));

mfdu=gf_mesh_fem(m,1);
% the P2 fem is not derivable across elements, hence we use a discontinuous
% fem for the derivative of U.
gf_mesh_fem_set(mfdu,'fem',gf_fem('FEM_PK_DISCONTINUOUS(3,1)'));
VM=get(b0, 'von mises',mds,mfdu);

U=get(mds, 'state'); U=U(1:get(mfu, 'nbdof'));

disp('plotting ... can also take some minutes!');

% we plot the von mises on the deformed object, in superposition
% with the initial mesh.
if (linear),
  gf_plot(mfdu,VM,'mesh','on', 'cvlst', get(m, 'outer faces'),...
	  'deformation',U,'deformation_mf',mfu);
else
  gf_plot(mfdu,VM,'mesh','on', 'cvlst', get(m, 'outer faces'),...
	  'deformation',U,'deformation_mf',mfu,'deformation_scale',1);
end;

caxis([0 100]);
colorbar; view(180,-50); camlight;
gf_colormap('tripod');

% the von mises stress is exported into a VTK file
% (which can be viewed with 'mayavi -d tripod.vtk -m BandedSurfaceMap')
% see http://mayavi.sourceforge.net/
gf_mesh_fem_get(mfdu,'export to vtk','tripod.vtk','ascii',VM,'vm')

次の図は Von Mises 応力を示しています．

変形三脚

ブリックフレームワークを使わない

モデルブリックは，最終的な線形システムのアセンブリのほとんどの詳細を隠蔽するため，非常に便利です．しかし，より低い水準で，線形システムのアセンブリと解析を直接 Octave か MatLab で処理することも可能です．例えば， demo_tripod_alt.m はアセンブリが明示的であることを除き， demo_tripod.m によく似ています．

nbd=get(mfd, 'nbdof');
F = gf_asm('boundary_source', 1, mim, mfu, mfd, repmat([0;-10;0],1,nbd));
K = gf_asm('linear_elasticity', mim, mfu, mfd, ...
           lambda*ones(1,nbd),mu*ones(1,nbd));

% handle Dirichlet condition
[H,R]=gf_asm('dirichlet', 2, mim, mfu, mfd, repmat(eye(3),[1,1,nbd]), zeros(3, nbd));
[N,U0]=gf_spmat_get(H, 'dirichlet_nullspace', R);
KK=N'*K*N;
FF=N'*F;
% solve ...
disp('solving...'); t0 = cputime;
lsolver = 1 % change this to compare the different solvers
if (lsolver == 1),     % conjugate gradient
  P=gfPrecond('ildlt',KK);
  UU=gf_linsolve('cg',KK,FF,P,'noisy','res',1e-9);
elseif (lsolver == 2), % superlu
  UU=gf_linsolve('superlu',KK,FF);
else                   % the matlab "slash" operator
  UU=KK \ FF;
end;
disp(sprintf('linear system solved in \%.2f sec', cputime-t0));
U=(N*UU).'+U0;

getfem-interface では，ベクトルと行列のアセンブリは gf_asm 関数により行われます．Dirichlet条件 \(u(x) = r(x)\) は，弱形式 \(\int (h(x)u(x)).v(x) = \int r(x).v(x)\quad\forall v\) で処理されます(ここで， \(h(x)\) は \(3\times 3\) 行列です．これは定数であり，単位と等価です)．元のシステムからDirichlet条件を除去することにより，縮約システム KK UU = FF を構築します．ペナルティ法を使用してDirichlet条件を処理する方が効率的(よりシンプル)であることに注意してください．

その他の例

• demo_refine.m スクリプトは，端部がクランプされた単純な2次元または3次元バーを示します．応力が特異な領域(クランプ領域とNeumann境界の間の遷移)でより良い近似を得るために適合リファインメントを用いました．

• demo_nonlinear_elasticity.m スクリプトは，曲げ捻りの3次元バーを示しています．これは，変形が多数のステップで適用されるため，準静的な問題です．各ステップで非線形(有限変形)弾性問題を解きます．

• demo_stokes_3D_tank.m スクリプトはタンク内のStokes(粘性流体)問題を示します． demo_stokes_3D_tank_draw.m はメッシュスライスとストリームラインを使用して，解の適切なプロットを描画する方法を示します． demo_stokes_3D_tank_alt.m は古い例で，廃止された gf_solve 関数を使っています．

• demo_bilaplacian.m スクリプトは， GetFEM の例である tests/bilaplacian.cc を改作したものです．bilaplacian(かKirchhoff-Loveのプレートモデル)を正方形で求解します．

• demo_plasticity.m スクリプトは， GetFEM の例である tests/plasticity.cc を応用したものです．2次元または3次元バーは多数のステップで曲げられ，材料の塑性が考慮されます(材料の Von Mises が所定の閾値を超えたときに塑化が起こります)．

• demo_wave2D.m は高次の幾何学的変換と高次のFEMを持つ2次元スカラー波動方程式の例(円柱による平面波の回折)です．

Octave/Matlabオブジェクト指向フィーチャーの使用

GetFEM ツールボックスの基本機能では，拡張 Octave か MatLab 機能は使用しません(ただし，getfemオブジェクトのハンドルは小さな構造体に格納されます)．しかしツールボックスには，ハンドルをカプセル化して本当の Octave / MatLab オブジェクトのように見せるオブジェクトの集合が付属しています．目的は特別な機能を提供することではなく，ツールボックスの一貫性を改善することです．

使用例を次に示します．

>> m=gf_mesh('cartesian',0:.1:1,0:.1:1)
m =
     id: 0
    cid: 0

>> m2=gfMesh('cartesian',0:.1:1,0:.1:1)
gfMesh object ID=1 [17512 bytes], dim=2, nbpts=121, nbcvs=100
% while \kw{m} is a simple structure, \kw{m2} has been flagged
% as  an object of class gfMesh.  Since the \texttt{display} method for
% these  objects  have  been  overloaded,  the  toolbox  displays  some
% information about the mesh instead of the content of the structure.
>> gf_mesh_get(m,'nbpts')
ans =
   121
% pseudo member access (which calls ##gf_mesh_get(m2,'nbpts'))
>> m2.nbpts
ans =
   121

詳細は，OOコマンドのリファレンス GetFEM OOコマンド を参照してください．