任意の項を計算する - 低水準の汎用的な構築手順 (非推奨)¶
このセクションでは, 現在では非推奨と見なされている GetFEM で実装された汎用構築プロシージャの最初のバージョンを説明します.これにより,線形の場合に任意の行列の集合を作成可能になります.非線形の場合,特殊な "non_linear_term" オブジェクトを実装する必要があります.これは複雑で, GetFEM という非常に低水準の内部ツールを使用する必要があります.汎用弱形式言語(GWFL)は,この困難を回避するために開発されました( 任意の項を計算する - 高水準の汎用的な構築手順 - 弱形式言語 (GWFL) 参照).
ファイル getfem/getfem_assembling.h
で見ることができるように,以前の構築手順はすべて getfem::generic_assembly
オブジェクトを使用しています.例えば,スカラーFEMのためのボリュームソース項の構築は,コードの以下の抜粋を用いて行われます.
getfem::generic_assembly assem;
assem.push_im(mim);
assem.push_mf(mf);
assem.push_mf(mfdata);
assem.push_data(F);
assem.push_vec(B);
assem.set("Z=data(#2);"
"V(#1)+=comp(Base(#1).Base(#2))(:,j).Z(j);");
assem.assembly();
最初の命令はオブジェクトを宣言し,使用するデータを設定します.積分法を保持する mesh_im オブジェクト,2つの mesh_fem オブジェクト,入力データ F
,および宛先ベクトル B
を含みます.
入力データは mfdata
で定義されているベクトル \(F\) です.1つは \(\sum_{j} f_j (\int_\Omega \phi^i \psi^j)\) を評価することです.命令はメッシュの凸の cv
ごとに実行されます.(項 #1
と項 #2
は最初の mesh_fem 2つ目は mf
と mfdata
です). Z=data(#2);
はそれぞれの凸について, push_data
で与えられた最初のデータ引数の値を push_data
に対応するインデックスで受け取ることを示します. 2番目の( #2
)の凸に付けられた自由度 mesh_fem (ここでは Z = F[mfdata.ind_dof_of_element(cv)]
となります).
V(#1)+=...
の部分は,次の式の結果が出力ベクトルに蓄積されることを意味します( push_vec
で提供されます).ここでもまた, #1
は,最初の( #1
) mesh_fem に関して,現在の凸の自由度に対応するインデックスに結果を書き込むことを意味します.
右側の comp(Base(#1).Base(#2))(:,j).Z(j)
には2つの操作が含まれています.最初のものは,凸のテンソルの計算です: comp(Base(#1).Base(#2))
は2次元のテンソルとして評価されます \(\int\phi^i \psi^j\) をすべての自由度について計算します mf
と現在の凸面に付けられた mfdata
の \(j\) .次の部分はリダクション演算で, C(:,j).Z(j)
:それぞれの指数(ここでは \(j\) )を合計します.結果は \(\sum_j c_{i,j} z_j\) となります.
comp(Base(#1).Base(#2))
の内部で使用される積分法は mim
から取得されます.別の mesh_im オブジェクトから積分法を使用する必要がある場合は,comp
の最初の引数として指定することができます,例えば, comp(\%2, Base(#1).Grad(#2))
は2番目の mesh_im オブジェクトを使用します(getfem++-2.0の新機能).
他の例として,ベクトルLaplacianの剛性行列を構築します.
getfem::generic_assembly assem;
assem.push_im(mim);
assem.push_mf(mf);
assem.push_mf(mfdata);
assem.push_data(A);
assem.push_mat(SM);
assem.set("a=data$1(#2);"
"M$1(#1,#1)+=sym(comp(vGrad(#1).vGrad(#1).Base(#2))(:,j,k,:,j,k,p).a(p))");
assem.assembly();
出力は, assem.push_mat(SM)
で挿入された疎行列に書き出されます. M$1(#1,#1)
の $1
は最初の行列を "push" していることを示しています(オプションですが,2つの行列を構築するならば,2つ目はこのように参照する必要があります). sym
関数は,結果が対称であることを保証します(これが行われないと,いくつかの丸め誤差が対称性を取り消し,構築が少し遅くなることに留意してください).次に, comp
部分は7次テンソルを評価し,
ここで, \(\varphi^i_j\) は mf
の \(i\) 次の基底関数の \(j\) 次成分で, \(\psi^p\) は(スカラー) 2番目の mesh_fem の基底関数を構築したいので
換算は次のとおりです.
comp
関数で Grad
の代わりに vGrad
が使われました.なぜなら, Laplacian ベクトル をアセンブルしているためです.そのため,それぞれの vGrad
部分に3次元の(自由度番号,成分番号,および派生番号)があります.スカラLaplacianの場合, comp(Grad(#1).Grad(#1).Base(#2))(:,k,:,k,p).a(p)
を使うことができます.しかし,ベクトル形式はベクトルとスカラーの両方の場合に機能する利点があります.
最後の命令である assem.assembly()
は各凸の式を評価します.境界を越えて構築を呼び出すには, assem.assembly(rg)
を呼び出します.ここで rg
は getfem::mesh_region
オブジェクト. rg
も数値であるかもしれません.その場合メッシュ領域は mim.linked_mesh().region(rg)
です.
3番目の例は,メッシュ境界上のスカラーまたはベクトルフィールドの \(L^2\) ノルムを計算する方法を示します.
assem.push_im(mim);
assem.push_mf(mf);
assem.push_data(U);
std::vector<scalar_type> v(1);
assem.push_vec(v);
assem.set("u=data(#1);"
"V()+=u(i).u(j).comp(vBase(#1).vBase(#1))(i,k,j,k)");
assem.assembly(boundary_number);
これは読みやすいです. assembly
が返ってくるとき, v[0]
は次の通りです.
4番目と最後の例は,完全なHookeテンソルによる線形弾性問題の(準最適な)構築を示しています.
assem.set("h=data$1(qdim(#1),qdim(#1),qdim(#1),qdim(#1),#2);"
"t=comp(vGrad(#1).vGrad(#1).Base(#2));"
"e=(t{:,2,3,:,5,6,:}+t{:,3,2,:,5,6,:}+t{:,2,3,:,6,5,:}+t{:,3,2,:,6,5,:})/4;"
"M(#1,#1)+= sym(e(:,j,k,:,m,n,p).h(j,k,m,n,p))");
元の方程式は次のとおりです.
\(h\) はHookeテンソル,そして \(:\) は行列間のスカラ積を意味します.定数ではないと仮定しているので,第2の mesh_fem: \(h_{ijkl}(x)=\sum_p h_{ijkl}^p \psi^p\) に \(h\) が与えられます.したがって,最初の行は,最初のデータが実際に5次元のテンソルであることを宣言します.最初の4番目のものはすべて最初の mesh_fem のターゲット次元に等しく,最後のものは次の第2の mesh_fem の自由度数に等しくなり, comp
部分はベクトルLaplacianの場合と同じ7次テンソルを計算します.このテンソルから,順列を使って \(\varepsilon(\varphi^i)_{jk}\varepsilon(\phi^l)_{mn}\psi^p\) を評価し最後に表現はHookeテンソルに換算されます.
comp
コマンドの中で利用可能な操作¶
Base(#i)
: i 次 mesh_fem の基底関数の値を評価します.Grad(#i)
: i 次 mesh_fem の基底関数の勾配の値を評価します.Hess(#i)
: i 次 mesh_fem の基底関数のHessianの値を評価します.Normal()
: 単位法線を評価します(体積積分には使用しないでください).NonLin$x(#mf1,... #mfn)
:push_nonlinear_term(pnonlinear_elem_term)
で挿入された,リストされた mesh_fem オブジェクトを使用した x 次非線形項の計算.GradGT()
,GradGTInv()
:現在の凸の幾何学的変換の勾配(およびその逆数)を評価します.
注釈
comp
コマンドの中の任意のデータオブジェクトを参照し, comp()
の中で縮約を実行することができます.この機能は非線形項の構築を高速化するために有用です(ファイル getfem/getfem_nonlinear_elasticity.h
の使用例を参照).
他の操作¶
スライスは縮約演算と混在することがあります. t(:,4,i,i)
は第2次元のインデックス4でスライスをとり,次元3と4の対角を縮約します. スライスのインデックス番号は1ではなく0で始まることに留意してください!!
mdim(#2)
は第2の mesh_fem に関連するメッシュ次元として評価され, qdim(#2)
は mesh_fem の対象次元です.
テンソルの対角は t{:,:,3,3}
で得ることができます.(厳密には t{1,2,3,3}
と同等です: コロンは読みやすさを向上させるためだけです).これは,並べ替え操作と同じ演算子です. t{:,:,1,1}
や t{:,:,4,4}
は有効な操作ではないことに注意してください.
print
コマンドは,テンソルを見るために使うことができます: "print comp(Base(#1));"
は各凸の基底関数の積分を出力します.
If there is more than one data array, output array or output sparse
matrix, one can use data$2
, data$3
, V$2
, M$2
,...