DelFEMでH面導波管伝達問題を解く
DelFEMでH面導波管伝達問題を解いてみました。
サンプルとしてDelFEM4Netに追加しています。
wg2d
本当は固有値解析も含めてDelFEMでやりたかったのですがちょっとハードルが高かったので固有値解析はKrdLab.Lisysを使用しました。
ん？メモリーが足りないといわれる？ちょっと分かりませんでした。ごめんなさい。
自作ラッパーでメモリーリークがありました。delete忘れです。ああ、ごめんなさい。
苦労した点：
DelFEMのリニア方程式ソルバーがFEMに特化している（だから計算は早い）ため、FEMの剛性行列、質量行列に変更を加えるような解析方法を選択するとかなり強引なことをする必要がありました。
計算した結果が最初は合わず、どこか間違っているのか調べたところ隣接要素以外の要素行列が全体行列へマージされないようでした。
一方H面導波管の伝達問題で使用される固有モード展開を境界条件に適用する方法は、境界に接する要素すべてについて関連性を持たせるので隣接要素以外へもマージしたい訳です。
幸いマージ対象の設定を行う低レベルのI/Fがあったのでそれを利用しました。
下記が該当箇所です。結構はまったので参考になればうれしいです。
/// <summary>
/// 境界のパターン追加(要素アレイ単位)
/// </summary>
/// <param name="ls">リニアシステム</param>
/// <param name="world">ワールド座標系</param>
/// <param name="fieldValId">フィールド値ID</param>
/// <param name="no_c_all">全体節点番号配列</param>
/// <param name="to_no_boundary">全体節点番号→境界上節点番号マップ</param>
/// <param name="eaId">要素アレイID</param>
/// <returns></returns>
private static bool makeWaveguidePortBCPattern_EachElementAry(CZLinearSystem ls, CFieldWorld world, uint fieldValId, uint[] no_c_all, Dictionary<uint, uint> to_no_boundary, uint eaId)
{
// フィールドを取得
CField valField = world.GetField(fieldValId);
System.Diagnostics.Debug.Assert(valField.GetInterpolationType(eaId, world) == INTERPOLATION_TYPE.LINE11);
// 要素アレイを取得
System.Diagnostics.Debug.Assert(world.IsIdEA(eaId));
CElemAry ea = world.GetEA(eaId);
System.Diagnostics.Debug.Assert(ea.ElemType() == ELEM_TYPE.LINE);
if (!world.IsIdField(fieldValId))
{
return false;
}
// 境界上の節点数
uint node_cnt = (uint)no_c_all.Length;
// 線要素の節点数
uint nno = 2;
// 要素節点の全体節点番号
uint[] no_c = new uint[nno];
// 要素剛性行列(コーナ-コーナー)
CZMatDia_BlkCrs_Ptr mat_cc = ls.GetMatrixPtr(fieldValId, ELSEG_TYPE.CORNER, world);
// 要素残差ベクトル(コーナー)
CZVector_Blk_Ptr res_c = ls.GetResidualPtr(fieldValId, ELSEG_TYPE.CORNER, world);
//System.Diagnostics.Debug.WriteLine("NBlkMatCol:" + mat_cc.NBlkMatCol()); // （境界でなく領域の総節点数と同じ?)
//System.Diagnostics.Debug.WriteLine("NBlkMatRow:" + mat_cc.NBlkMatRow());
//System.Diagnostics.Debug.WriteLine("LenBlkCol:" + mat_cc.LenBlkCol());
//System.Diagnostics.Debug.WriteLine("LenBlkRow:" + mat_cc.LenBlkRow());
//System.Diagnostics.Debug.WriteLine("node_cnt:" + node_cnt);
// どうもマトリクスにマージされていない 節点を含む隣接要素以外のマージができないようになっている？
// パターンを追加してみる
using (CIndexedArray crs = new CIndexedArray())
{
crs.InitializeSize(mat_cc.NBlkMatCol());
//crs.Fill(mat_cc.NBlkMatCol(), mat_cc.NBlkMatRow());
using (UIntVectorIndexer index = crs.index)
using (UIntVectorIndexer ary = crs.array)
{
for (int iblk = 0; iblk < (int)crs.Size(); iblk++)
{
index[iblk] = 0;
}
for (int iblk = 0; iblk < mat_cc.NBlkMatCol(); iblk++)
{
// 現在のマトリクスのインデックス設定をコピーする
uint npsup = 0;
ConstUIntArrayIndexer cur_rows = mat_cc.GetPtrIndPSuP( (uint)iblk, out npsup);
foreach (uint row_index in cur_rows)
{
ary.Add(row_index);
}
index[iblk + 1] = (uint)ary.Count;
// 境界の節点に関して列を追加
int ino_boundary = -1;
int col = -1;
for (uint ino_boundary_tmp = 0; ino_boundary_tmp < node_cnt; ino_boundary_tmp++)
{
int tmp_col = (int)no_c_all[ino_boundary_tmp];
if (tmp_col == iblk)
{
col = tmp_col;
ino_boundary = (int)ino_boundary_tmp;
break;
}
}
if (col != -1 && ino_boundary != -1)
{
// 関連付けられていない節点をその行の列データの最後に追加
int last_index = (int)index[col + 1] - 1;
int add_cnt = 0;
for (uint jno_boundary = 0; jno_boundary < node_cnt; jno_boundary++)
{
uint row = no_c_all[jno_boundary];
if (ino_boundary != jno_boundary) // 対角要素は除く
{
if (!cur_rows.Contains(row))
{
ary.Insert(last_index + 1 + add_cnt, row);
add_cnt++;
//System.Diagnostics.Debug.WriteLine("added:" + col + " " + row);
}
}
}
if (add_cnt > 0)
{
//ndex[col + 1] += (uint)add_cnt;
index[col + 1] = (uint)ary.Count;
}
}
}
}
crs.Sort();
System.Diagnostics.Debug.Assert(crs.CheckValid());
mat_cc.AddPattern(crs);
}
return true;
}
