要素メッシュの生成
メッシュ生成の機能を使うためには,有限要素法(FEM)のパッケージをロードしなければならない.
Needs["NDSolve`FEM`"]はじめに
数値解法の多くでは,関心領域をその領域の近似で置き換えて使う.この近似は,離散領域と呼ばれる.離散領域は,より小さい要素の集まりに分割され,これらの要素を合計すると,離散領域全体になる.この分割された離散領域は,メッシュと呼ばれる.そうすると,数値解を求めることは,より小さな要素についての解を計算してから,部分的な解をメッシュ全体の解になるように組み合せることに基づく.
NDSolveは例えば,内部で領域
をElementMeshオブジェクトに変換する.このElementMeshは,数値解析が行われる領域
の離散的な近似バージョンである.NDSolve等の数値関数は,記号的な領域記述の代りに,ElementMeshを入力として受け取ることができる.このことは,メッシュ生成のプロセスにおいて大きな柔軟性を与える.実際要素メッシュは,例えば外部ツールによって生成することも可能になる.
いろいろな方法でElementMeshを作成することができ,さまざまな関数がメッシュ作成のプロセス中の手助けをするために提供されている.ElementMeshを作成する主な関数として,ToElementMeshがある.ToElementMeshは,概念的に異なるさまざまなメソッドでElementMeshを作成することを可能にする.
陰的な要素メッシュの作成は,ImplicitRegion等の陰関数をElementMeshに変換することに基づく.それに対して,明示的なメッシュ生成は,GraphicsComplex等の明示的な表現を要素メッシュに変換することに基づく.手作業でのメッシュ生成も,この明示的な生成の一種である.ここでは,メッシュ要素の明示的な集合を与えて,要素メッシュを形成する.
陰的な要素メッシュ生成も,明示的な要素メッシュ生成も,さらに小さな部分に分割することが可能である.関数ToBoundaryMesh は,陰的あるいは明示的な入力の境界表現を生成する.この境界表現を今度はToElementMeshに与えて完全なメッシュを形成する.この他の用法の中でも,数値メソッドが境界表現のみを必要とする場合には,ToBoundaryMeshが便利である.
一つ心に留めておかなければならない重要なことは,どのように要素メッシュが作成された場合でも,簡単な場合を除いて,要素メッシュは厳密な領域の近似に過ぎないという点である.要素メッシュが領域の重要な部分を捉える厳密さは,数値解の質に大きく係わってくる.
いったんメッシュが作成されると,NDSolve等の数値関数にこれを渡すことができる.
ElementMeshをNDSolveに渡す
最初の例として,ディリクレ(Dirichlet)境界条件を持つポアソン(Poisson)偏微分方程式を標準的な方法で解く.このプロセスの間に,メッシュが内部的に生成される.次のステップでは,このメッシュがNDSolveへの引数として事前定義され与えられる.
Ω = ImplicitRegion[True, {{x, 0, 2}, {y, 0, 1}}]op = -Laplacian[u[x, y], {x, y}] - 20Γ = DirichletCondition[u[x, y] == 0, x == 0 || x == 2]ufun = NDSolveValue[{op == 0, Γ}, u, {x, y}∈Ω]Show[
ContourPlot[ufun[x, y], {x, y}∈Ω, AspectRatio -> Automatic],
ufun["ElementMesh"]["Wireframe"]]NDSolveが有限要素法を利用した場合には,Interpolation関数がElementMeshを保存することに注意する.ElementMeshが補間関数とともに保存された場合には,それを抽出することができる.
ufun["ElementMesh"]陰的なパラメトリック領域を指定する代りに,明示的なElementMeshを指定することも可能である.これは,ToElementMeshを使って行うことができる.
mesh = ToElementMesh["Coordinates" -> {{0., 0.}, {1., 0.}, {2., 0.}, {2., 1.}, {1., 1.}, {0., 1.}}, "MeshElements" -> {TriangleElement[{{1, 2, 5}, {5, 6, 1}, {2, 3, 4}, {4, 5, 2}}]}]mesh["Wireframe"]要素メッシュの可視化についての詳細は,要素メッシュの可視化についてのチュートリアルを参照されたい.
次に,同じ偏微分方程式を解く.今回は,明示的なメッシュのみが定義されている.
ufun = NDSolveValue[{op == 0, Γ}, u[x, y], {x, y}∈mesh];Show[
ContourPlot[ufun, {x, y}∈mesh, AspectRatio -> Automatic],
mesh["Wireframe"]]ほとんどの場合,4つの要素だけでは,正確な解を表すのには十分ではないことを心に留めておくべきである.
MeshOptionsを通してElementMesh作成のオプションをNDSolveに渡す
ToElementMeshとToBoundaryMeshのオプションはすべて,NDSolveに直接与えることができる.
ufun = NDSolveValue[{op == 0, Γ}, u, {x, y}∈Ω, Method -> {"FiniteElement", "MeshOptions" -> {MaxCellMeasure -> 0.01}}];代りに,要素メッシュをシミュレーションの前に生成し,NDSolveに与えることも可能である.
mesh = ToElementMesh[Ω, MaxCellMeasure -> 0.01];
mesh["Wireframe"]ufun = NDSolveValue[{op == 0, Γ}, u, {x, y}∈mesh]ufun = NDSolveValue[{D[u[t, x, y], t] - Laplacian[u[t, x, y], {x, y}] - 20 == 0, DirichletCondition[u[t, x, y] == 0, x == 0 || x == 2], u[0, x, y] == 0}, u, {t, 0, 1}, {x, y}∈Ω, Method -> {"PDEDiscretization" -> {Automatic, "SpatialDiscretization" -> {"FiniteElement", "MeshOptions" -> {MaxCellMeasure -> 0.01}}}}]NDSolveのオプションの指定とその解法段階についての詳細は,NDSolveの「詳細とオプション」のセクションに記載されている.
NDSolveにメッシュが与えられた場合には,メッシュオプションは使っても効果がなくなる.
ElementMeshとMeshRegionを比較する
ElementMeshの機能についてさらに詳しく説明する前に,ElementMeshオブジェクトをMeshRegionオブジェクトと比べてみるとよいかも知れない.
ElementMeshのデフォルトの"MeshOrder"は,2である.
ToElementMesh[Disk[]]["MeshOrder"]メッシュ次数を通して表現されたより高次の要素を取り扱えるということには,2つの明らかな利点がある.
- 必要であれば,例えばMeshRegionへの変換が簡単にできる
- 簡単.1つずつのコンバータ(ToBoundaryMeshとToElementMesh)だけ
pde = D[-2u[x], x, x] + 3 D[u[x], x] + 10u[x] == 1;
exact = DSolveValue[{pde, 2 * u[-1] == -1, 3 * u[2] == 1 / 2}, u, x];
ufun1 = NDSolveValue[{pde, DirichletCondition[2 * u[x] == -1, x == -1],
DirichletCondition[3 * u[x] == 1 / 2, x == 2]}, u, {x}∈ToElementMesh[FullRegion[1], {{-1, 2}}, "MeshOrder" -> 1]];
ufun2 = NDSolveValue[{pde, DirichletCondition[2 * u[x] == -1, x == -1],
DirichletCondition[3 * u[x] == 1 / 2, x == 2]}, u, {x}∈ToElementMesh[FullRegion[1], {{-1, 2}}, "MeshOrder" -> 2]];
Plot[{exact[x] - ufun1[x], exact[x] - ufun2[x]}, {x, -1, 2}]mr = DiscretizeRegion[Disk[]];
Head[mr]em = ToElementMesh[Disk[]]π - RegionMeasure[mr]π - Total[em["MeshElementMeasure"], 2]ElementMeshは,二次近似を使うので,Diskをよりよく近似することができる.
MeshRegionとElementMeshの間の変換は簡単である.
ToElementMesh[mr]MeshRegion[em]もとのElementMeshが一次(デフォルト)よりも高い次数である場合には,変換されたMeshRegionを使ったときに期待通りの結果が得られないこともある.これが起こらないようにする1つの方法として,まずElementMeshを一次メッシュに変換してから,次にMeshRegionに変換するという方法がある.
em = ToElementMesh[Disk[]];
mr = MeshRegion[MeshOrderAlteration[em, 1]]GraphicsRow[{RegionPlot[mr], RegionPlot[em]}]MeshRegionのデータ構造と比べた利点を提供するためには,いくつかの注意点を考慮する必要がある.ElementMeshでは,境界メッシュと完全なメッシュの違いは,完全なメッシュ要素の有無によって示される.境界要素メッシュは,完全なメッシュ要素に対してAutomaticを設定する.
bmesh = ToBoundaryMesh[Disk[]]mesh = ToElementMesh[Disk[]]bmesh["MeshElements"]mesh["MeshElements"]//Short境界のElementMeshは,閉じられた曲線である必要はない.
ToBoundaryMesh["Coordinates" -> {{0, 0}, {1, 0}, {1, 1}}, "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}}]}]["Wireframe"]これは,完全メッシュを必要としない数値的アルゴリズム(例:境界積分)に必要である.閉じていない境界曲線は,しかしながら,完全な要素メッシュに変換することができない.このことがBoundaryMeshRegionとどのように違うかということに注意されたい.BoundaryMeshRegionは常に,境界を使ってその領域を囲むことによって,完全領域を表す.
ElementMeshとMeshRegionの間でさらに違う点は,MeshRegionには,完全な次元のメッシュ領域要素から切り離される低次元の成分が含まれる場合があるということである.
mr = MeshRegion[{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {1, 2 / 3}, {2, 2 / 3}}, {Polygon[{{1, 2, 3, 4}}], Line[{{5, 6}}]}]DimensionalMeshComponents[mr]ToElementMesh[Last[DimensionalMeshComponents[mr]]]["Wireframe"]境界要素メッシュには,内部構造(例えば,2つの物質領域を表すため)が含まれることがある.
bmesh = ToBoundaryMesh["Coordinates" -> {{0, 0}, {1, 0}, {1, 1}, {0, 1}, {1 / 6, 1 / 6}, {5 / 6, 1 / 6}, {5 / 6, 5 / 6}}, "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 1}}], LineElement[{{5, 6}, {6, 7}}]}];
bmesh["Wireframe"]mesh = ToElementMesh[bmesh, MaxCellMeasure -> Infinity];
mesh["Wireframe"]内部構造がどのように最終的なElementMeshにまだ存在しているかについて注意する.
以下は,ElementMeshを直接使用する場合の利点と問題点のリストである.残りの項目のいくつかについては,このチュートリアルでさらに説明する.
ElementMeshでの領域の近似
ElementMeshを作成するには,完全なメッシュにToElementMeshを使う,あるいは境界メッシュ表現にToBoundaryMeshを使う方法がある.
ToElementMeshが呼び出されると,ToBoundaryMeshがまず内部で呼び出される.ToElementMeshですべてを行うことも可能であるが,完全なメッシュを生成する前に,境界メッシュをチェックすると便利なこともある. ToBoundaryMeshは特に境界上のマーカーを使う場合に有用である.
手作業のメッシュ作成
最小のElementMeshは,座標と要素からなる.次の要素が利用できる.
- 1D:LineElement
要素は,そのタイプ(例えば,TriangleElement)と整数リストのリストによって指定される.
e1 = TriangleElement[{{1, 2, 3}, {3, 4, 1}}]整数は,座標に対応する指標(インシデントとも呼ばれる)である.上の場合では,TriangleElementには,インシデント{1,2,3}からなる要素とインシデント {3,4,1}からなる要素の2つの三角要素が含まれる.
メッシュを構築するには,座標を与えなければならない.座標は,1D,2D,3Dのいずれでもよいが,要素タイプに合ったものでなければならない.三角メッシュには,2D座標が与えられなければならない.
coordinates = {{0., 0.}, {1., 0.}, {1., 1.}, {0., 1.}}要素指標は,座標に対応する.2つ目の三角要素のインシデントは{3,4,1}であり,それぞれ座標{1.,1.},{1.,0.},{0.,0.}を参照する.
mesh = ToElementMesh["Coordinates" -> coordinates, "MeshElements" -> {e1}]要素インシデントの概念は,GraphicsComplexの概念に密接に関係している.
Show[
mesh["Wireframe"["MeshElementIDStyle" -> DarkGreen]],
mesh["Wireframe"["MeshElement" -> "PointElements", "MeshElementIDStyle" -> Red]]]線分メッシュ
1Dメッシュとしては,メッシュ要素 ei はLineElementである.境界要素 bi はPointElementである.
mesh = ToElementMesh["Coordinates" -> {{0.}, {0.5}, {1.}}, "MeshElements" -> {LineElement[{{1, 2}, {2, 3}}]}]Show[
mesh["Wireframe"["MeshElementIDStyle" -> DarkGreen]],
mesh["Wireframe"["MeshElement" -> "PointElements", "MeshElementIDStyle" -> Red]]]三角メッシュ
2Dメッシュとしては,メッシュ要素 ei はTriangleElementまたはQuadElementである.境界要素 bi はLineElementである.
三角メッシュを作成するには,座標と三角要素が必要である.線形三角形には3つのインシデントが含まれ,要素タイプはTriangleElementである.インシデント中の整数の数がメッシュの次数を決定する.三角形の場合は,1つの要素につき3つのインシデントが線形三角形に対応し,この三角形が一次メッシュに対応する.インシデントは反時計回りに与えられなければならない.
mesh = ToElementMesh["Coordinates" -> {{1.293, 0.228}, {1., 0.}, {0.94, 0.342}, {1.293, 0.}, {1.215, 0.442}, {2., 0.}, {1.879, 0.684}}, "MeshElements" -> {TriangleElement[{{1, 3, 2}, {1, 2, 4}, {1, 4, 6}, {1, 6, 7}, {1, 7, 5}, {1, 5, 3}}]}]mesh["MeshOrder"]Show[
mesh["Wireframe"],
mesh["Wireframe"["MeshElement" -> "PointElements", "MeshElementIDStyle" -> Red]]]Show[
mesh["Wireframe"["MeshElementIDStyle" -> DarkGreen]],
mesh["Wireframe"["MeshElement" -> "PointElements", "MeshElementIDStyle" -> Red]]]第1要素の第1の辺はノード{3,2}を繋ぐもので,近傍要素を持たないので,0である.ノード2を中央のノードに繋ぐ第2の辺は,第2要素に繋がれている.第3要素は,中央のノードをノード3に繋ぐ辺を持ち,この辺は第6要素に繋がれている.第4要素以降の要素すべてについても同じように取り扱われる.
mesh["ElementConnectivity"]TriangleElement等の要素コンテナは,マーカーも含むことができる.これらは,別の物質領域に印を付ける際に便利である.マーカーの数は,要素の数と同じでなければならない.
mesh = ToElementMesh["Coordinates" -> {{1.293, 0.228}, {1., 0.}, {0.94, 0.342}, {1.293, 0.}, {1.215, 0.442}, {2., 0.}, {1.879, 0.684}}, "MeshElements" -> {TriangleElement[{{1, 3, 2}, {1, 2, 4}, {1, 4, 6}, {1, 6, 7}, {1, 7, 5}, {1, 5, 3}}, {1, 1, 1, 2, 2, 2}]}]mesh["Wireframe"["MeshElementMarkerStyle" -> DarkGreen]]ここまでは,三角メッシュは一次であった.二次三角メッシュは,1つの要素について6つのインシデントを持つ.追加の座標は,中間のノードである.したがって,二次三角要素には,6つのインシデントが含まれる.最初の3つは線形インシデントであり,次の3つは二次インシデントである.TriangleElementを参照されたい.
coordinates = {{1.293, 0.228}, {1., 0.}, {0.94, 0.342}, {1.293, 0.}, {1.215, 0.442}, {2., 0.}, {1.879, 0.684}, {0.97, 0.171}, {1.0775, 0.392}, {1.1165, 0.285}, {1.1465, 0.}, {1.1465, 0.114}, {1.254, 0.335}, {1.293, 0.114}, {1.547, 0.563}, {1.586, 0.456}, {1.6465, 0.}, {1.6465, 0.114}, {1.9395, 0.342}};e1 = TriangleElement[{{1, 3, 2, 10, 8, 12}, {1, 2, 4, 12, 11, 14}, {1, 4, 6, 14, 17, 18}, {1, 6, 7, 18, 19, 16}, {1, 7, 5, 16, 15, 13}, {1, 5, 3, 13, 9, 10}}];mesh = ToElementMesh["Coordinates" -> coordinates, "MeshElements" -> {e1}]Show[
mesh["Wireframe"],
mesh["Wireframe"[ "MeshElement" -> "PointElements", "MeshElementIDStyle" -> Red]]]要素メッシュの可視化の詳細は,要素メッシュの可視化についてのチュートリアルを参照されたい.
クワッドメッシュ
QuadElementメッシュは,TriangleElementメッシュとまったく同じように働く.唯一の違いは,線形クワッド要素については,1つの要素につき4つのインシデントが必要で,二次要素については,1つの要素につき8つのインシデントが必要な点である.
nx = ny = 10;
coordinates = Flatten[ Table[{r Cos[θ], r Sin[θ]}, {r, 1., 2., 1 / (ny - 1)}, {θ, 0., 2 Pi / 3., (2 Pi / 3.) / (nx - 1)}], 1];incidents = Flatten[Table[{j * nx + i, j * nx + i + 1, (j - 1) * nx + i + 1, (j - 1) * nx + i}, {i, 1, nx - 1}, {j, 1, ny - 1}], 1];mesh = ToElementMesh["Coordinates" -> coordinates, "MeshElements" -> {QuadElement[incidents]}]mesh["Wireframe"]マーカーをTriangleElementの場合とまったく同じ方法で与えることもできる.
2Dにおける複合要素タイプのメッシュ
2Dメッシュについては,メッシュ要素 ei は,TriangleElementとQuadElementを組み合せたものでもよい.境界要素 bi はLineElementである.
coordinates = {{0., 0.}, {1., 0.}, {2., 0.}, {2.5, 0.5}, {0., 1.}, {1., 1.}, {2., 1.}, {3., 1.}, {2.5, 1.5}, {0., 2.}, {1., 2.}, {2., 2.}};e1 = QuadElement[{{1, 2, 6, 5}, {2, 3, 7, 6}, {5, 6, 11, 10}, {6, 7, 12, 11}}];
e2 = TriangleElement[{{3, 4, 7}, {4, 8, 7}, {7, 9, 12}, {7, 8, 9}}];mesh = ToElementMesh["Coordinates" -> coordinates, "MeshElements" -> {e1, e2}]mesh["Wireframe"]複合要素メッシュのすべての要素は,同じ次数でなければならない.一次三角要素と二次クワッド要素を同じメッシュに持つことは不可能である.
ele = {QuadElement[ElementIncidents[e1], {1, 1, 2, 2}], TriangleElement[ElementIncidents[e2], {1, 1, 2, 2}]};mesh = ToElementMesh["Coordinates" -> coordinates, "MeshElements" -> ele]mesh["Wireframe"["MeshElementMarkerStyle" -> DarkGreen]]要素の接続性は,マーカー境界についての情報を保つ.2つの繋がれている要素のマーカー値が異なる場合は常に,その要素の要素接続性エントリは負である.
mesh["Wireframe"["MeshElementMarkerStyle" -> Red, "MeshElementIDStyle" -> DarkGreen, "ContinuousElementID" -> True]]上のグラフィックスでは,マーカー1を持つ要素番号2が,マーカー2を持つ要素番号4と繋がっている.このマーカー地におけるジャンプの要素接続性は,負の符号に保存される.
mesh["ElementConnectivity"]境界層メッシュとしても使える2D複合要素メッシュの例は,偏微分方程式のモデルコレクションに記載されている.
2Dにおける境界メッシュ
境界メッシュは,完全なメッシュを生成する場合に有用である.2Dでは,境界要素 bi はLineElementである.
bmesh = ToBoundaryMesh["Coordinates" -> {{0., 0.}, {1., 0.}, {2., 0.}, {2., 1.}, {1., 1.}, {0., 1.}}, "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 5}, {5, 6}, {6, 1}}]}]bmesh["Wireframe"]ToElementMesh[bmesh]["Wireframe"]bmesh = ToBoundaryMesh["Coordinates" -> {{0., 0.}, {1., 0.}, {2., 0.}, {2., 1.}, {1., 1.}, {0., 1.}}, "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 5}, {5, 6}, {6, 1}, {2, 5}}]}]bmesh["Wireframe"]ToElementMesh[bmesh]["Wireframe"]境界メッシュには,メッシュの一部になるように,境界領域内に任意の点を加えることができる.
coordinates = Join[{{0., 0.}, {1., 0.}, {2., 0.}, {2., 1.}, {1., 1.}, {0., 1.}}, DeleteDuplicates[Table[{1. + 0.2 Sin[ϕ], 0.5 + 0.2 Cos[ϕ]}, {ϕ, 0, 2 π, π / 25}]]];bmesh = ToBoundaryMesh["Coordinates" -> coordinates, "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 5}, {5, 6}, {6, 1}}]}]bmesh["PointElements"]Show[
bmesh["Wireframe"],
bmesh["Wireframe"["MeshElement" -> "PointElements"]]]ToElementMesh[bmesh]["Wireframe"]四面体メッシュ
三次元の手作業によるメッシュ作成は,一次元や二次元の場合と同じような手順を踏む.
coordinates = N[Join[{{0, 0, 1}, {0, 0, 0}}, Table[{Sin[i], Cos[i], 0}, {i, 0, 2π - π / 4, π / 4}]]];Graphics3D[MapIndexed[Text[ToString[#2[[1]]], #1]&, coordinates]]mesh = ToElementMesh["Coordinates" -> coordinates, "MeshElements" -> {
TetrahedronElement[{{1, 2, 3, 4}, {1, 2, 4, 5}, {1, 2, 5, 6}, {1, 2, 6, 7}, {1, 2, 7, 8}, {1, 2, 8, 9}, {1, 2, 9, 10}, {1, 2, 10, 3}}, {1, 1, 1, 1, 2, 2, 2, 2}]
}]mesh["Wireframe"["MeshElement" -> "MeshElements", "MeshElementMarkerStyle" -> Red, Boxed -> False]]六面体メッシュ
nx = 13;ny = 6;nz = 7;
coordinates = Flatten[ Table[{r Cos[θ], r Sin[θ], h}, {h, -1., 1., 2 / (nz - 1)}, {r, 1., 2., 1 / (ny - 1)}, {θ, 0., 2 Pi / 3., (2 Pi / 3.) / (nx - 1)}], 2];incidents = Flatten[Table[Block[{p1 = (j - 1) * nx + i, p2 = j * nx + i, p3 = p2 + 1, p4 = p1 + 1, p5, p6, p7, p8},
{p5, p6, p7, p8} = {p1, p2, p3, p4} + k * nx * ny;
{p1, p2, p3, p4} += (k - 1) * nx * ny;
{p1, p2, p3, p4, p5, p6, p7, p8}], {i, 1, nx - 1}, {j, 1, ny - 1}, {k, 1, nz - 1}], 2];mesh = ToElementMesh["Coordinates" -> coordinates, "MeshElements" -> {HexahedronElement[incidents]}]mesh["Wireframe"]3Dにおける複合要素タイプのメッシュ
3Dメッシュについては,メッシュ要素 ei は,TetrahedronElementかHexahedronElementのいずれでもあり得るが,PrismElementでそれらが繋がれていない限り,どちらもであることはない.
3Dにおける境界メッシュ
境界メッシュは,完全なメッシュを生成するのに有用である.3Dでは,境界要素 bi はTriangleElementかQuadElementである.
coordinates = {{0., 0., 0.}, {2., 0., 0.}, {2., 1., 0.}, {0., 1., 0.}, {0., 0., 1.}, {2., 0., 1.}, {2., 1., 1.}, {0., 1., 1.}}boundaryElements = {QuadElement[{{1, 2, 3, 4}, {5, 6, 7, 8}, {2, 6, 7, 3}, {3, 7, 8, 4}, {4, 8, 5, 1}, {1, 5, 6, 2}}]}bmesh = ToBoundaryMesh["Coordinates" -> coordinates, "BoundaryElements" -> boundaryElements]bmesh["Wireframe"]ToElementMesh[bmesh, MaxCellMeasure -> Infinity]["Wireframe"]記号領域
記号領域の表現の使用と,その表現からElementMeshへの変換は,簡単である.任意のグラフィックスプリミティブを使って記号領域が表現できる.
Ω = Disk[]mesh = ToElementMesh[Ω]mesh["Wireframe"]ToElementMesh[Torus[]]["Wireframe"]ToElementMesh[RegionDifference[Cuboid[], Ball[]]]["Wireframe"]ToBoundaryMeshとToElementMeshにあるその他の多くの例を参照されたい.
数値領域
NumericalRegionは,記号的な領域表現と境界,および完全な要素メッシュを組み合せるためのものである.これは,現在自動化されているメッシュ生成を超える柔軟性を追加してくれる.例えば,ToBoundaryMeshの"MaxBoundaryCellMeasure"オプションを使うと,最大境界セルの大きさを制御することが可能であるが,このオプションでは境界の一部だけを細分化することはできない.そうすると,境界全体を細分化することによって,不要なメッシュ要素が数多く生成されてしまうこともある.この問題を回避する一つの方法に,手動で境界メッシュを生成する方法がある.さらに,領域の記号的表現があれば,最小要素を含む高品質な二次メッシュを生成することが可能である.
Ω = RegionDifference[Disk[{0, 0}, 1000], Disk[{0, 0}, 1]];mesh = ToElementMesh[Ω]{mesh["Wireframe"], mesh["Wireframe"[PlotRange -> {{-1.2, 1.2}, {-1.2, 1.2}}]]}この場合には,メッシュ生成過程において,領域中央の小さな穴を細分化した.PrecisionGoalを増やすことによって,さらに細分化することは可能であるが,必ずしも最適の解決法とは限らない.
mesh = ToElementMesh[Ω, PrecisionGoal -> 7]{mesh["Wireframe"], mesh["Wireframe"[PlotRange -> {{-1.2, 1.2}, {-1.2, 1.2}}]]}メッシュ要素の数が大幅に増えた.これを回避する一つの方法に,記号領域の記述を組み合せて,手動で境界メッシュを生成し, それらを最終的なメッシュのNumericalRegionに組み合せる方法がある.NumericalRegionは,ToNumericalRegionで生成される.
nr = ToNumericalRegion[Ω]手動で境界メッシュを生成するということは,領域の各部分を別々に生成してから,新しい境界メッシュでそれらを組み合せるということである.
bm1 = ToBoundaryMesh[Disk[{0, 0}, 1000]];
bm2 = ToBoundaryMesh[Disk[{0, 0}, 1], PrecisionGoal -> 5];bele1 = bm1["BoundaryElements"];
bele2 = bm2["BoundaryElements"];新しい境界メッシュは,両方のメッシュからの座標を結合してから新しい境界要素のリストを構築することにより生成する.新しい境界要素リストは,1つ目の境界メッシュからの境界要素と,2つ目の境界メッシュからの境界要素の指標を移動させたものを組み合せたものである.この場合,2つ目のメッシュの座標が1つ目のメッシュの座標の後ろに隠れてしまうので,指標は移動させる必要がある.
bmesh = ToBoundaryMesh["Coordinates" -> Join[bm1["Coordinates"], bm2["Coordinates"]], "BoundaryElements" -> Join[bele1, MapThread[#1[#2]&, {Head /@ bele2, Length[bm1["Coordinates"]] + ElementIncidents[bele2]}]], "RegionHoles" -> {{0, 0}}];{bmesh["Wireframe"], bmesh["Wireframe"[PlotRange -> {{-1.2, 1.2}, {-1.2, 1.2}}]]}SetNumericalRegionElementMesh[nr, bmesh]nr["SymbolicRegion"]nr["BoundaryMesh"]mesh = ToElementMesh[nr]{mesh["Wireframe"], mesh["Wireframe"[PlotRange -> {{-1.2, 1.2}, {-1.2, 1.2}}]]}mesh["MeshOrder"]境界要素メッシュを結合させることは,境界要素メッシュが分かれている場合に最もうまくいく.交点は通常大丈夫だが,2Dにおいて重複している辺あるいは3Dにおいて重複している面については,問題になることがある.このアプローチは,部分領域の大きさの違いが大きい場合に使うとよい.その他の場合には,部分領域の作成は, 次のセクションで説明するように,ブール演算を"RegionHoles"と"RegionMarker"オプションと組み合せた形で処理する方がよい.
特別な目的のメッシュ
領域の積
2Dおよび3Dのメッシュを作成するもう一つの方法に,領域の積を利用することがある.1Dメッシュに別の1Dメッシュを掛けると2Dメッシュが作成され,2Dメッシュに1Dメッシュを掛けると3Dメッシュが作成される.
mesh1D = ToElementMesh["Coordinates" -> {{0.}, {0.5}, {1.}}, "MeshElements" -> {LineElement[{{1, 2}, {2, 3}}]}]mesh2D = ElementMeshRegionProduct[mesh1D, mesh1D]mesh2D["Wireframe"]mesh = ToElementMesh[Annulus[], "MaxCellMeasure" -> 0.1];
mesh["Wireframe"]ElementMeshRegionProduct[mesh, ToElementMesh[Line[{{0}, {1 / 5}}], "MaxCellMeasure" -> 0.1]]["Wireframe"]グレーデッドメッシュ
関数ToGradedMeshは,1Dのグレーデッドメッシュを生成する.グレーデッドメッシュは,ノードが非一様に分布するメッシュである.
mesh = ToGradedMesh[Line[{{0}, {2}}], <|"Alignment" -> "Left"|>]MeshRegion[mesh]左にノードがより集中していることに注目されたい."Right","BothEnds","Central","Uniform"等,その他の形式のアラインメントも可能である.またカスタムのアラインメント関数を与えることもできる.
グレーデッドメッシュは,数多くの状況に役立つ.無限の範囲を模倣するために,比較的大きな領域をモデル化して,1つの部分にだけノードが集中するようにできる.これは,ToGradedMeshの関数ページの例にある通りである.グレーデッドメッシュが便利である別の分野には,材料の境界が不連続である場合がある.問題のある領域においてノードが集中しているグレーデッドメッシュを生成すると,解の近似の質が向上し,ときにはかなり向上することもある.
ElementMeshRegionProductを使うと,2Dおよび3Dのグレーデッドメッシュを作成することもできる.
meshX = ToGradedMesh[Line[{{0}, {2}}], <|"Alignment" -> "Central"|>];meshY = ToGradedMesh[Line[{{0}, {1}}], <|"Alignment" -> "Central", "ElementCount" -> 50|>];productMesh = ElementMeshRegionProduct[meshX, meshY]productMesh["Wireframe"]meshZ = ToGradedMesh[Line[{{0}, {1}}], <| "Alignment" -> "BothEnds", "ElementCount" -> 50, "MinimalDistance" -> 1 / 200|>]productMesh = ElementMeshRegionProduct[productMesh, meshZ]productMesh["Wireframe"]その他の例については,ToGradedMeshおよびElementMeshRegionProductのアプリケーションセクションを参照されたい.
完全整合層のメッシュ
完全整合層(PML) は,メッシュと対応する偏微分方程式を拡張し,PML内の信号の減衰を可能にする.これは,無限に拡張された領域を模倣するが,有限サイズの領域を使うという点で便利である.これらのタイプのメッシュの作成は,特定の例を見ることによって説明すると最も分かりやすい.音響学については,時間領域におけるPMLと周波数領域におけるPMLの使用と作成についての説明がある.3DのPMLの電磁気の例もある.
等高線からのメッシュ
等高線があるメッシュあるいは境界線を持つオブジェクトのメッシュが必要になることもある.そのようなメッシュは流動的なシミュレーションに使われることが多い.例えば,平面等のオブジェクトの周りの流体のシミュレーションを行う場合がそうである.このセクションでは,ピコ島の周りのメッシュがどのように作成されるかを示す.データは,島の地理オブジェクトから取り出したものである.
pico = GeoPosition[{38.5, -28.3}];
geoRange = GeoRange -> {{38.35, 38.6}, {-28.6, -28}};
geoGraphics = GeoGraphics[pico, geoRange]elevation = QuantityMagnitude[GeoElevationData[Options[geoGraphics, {GeoRange, GeoProjection, GeoZoomLevel}]]];hight = 248;length = 1500;
scaledElevation = Clip[elevation, {0, Max[elevation]}] / (Max[elevation] / hight);step = 10;
scaledElevationSample = scaledElevation[[1 ;; -1 ;; step, 1 ;; -1 ;; step]];listPlot = ListPlot3D[Reverse[scaledElevationSample], Boxed -> False, Axes -> False, DataRange -> {{0, length}, {0, length}}, Mesh -> None, PlotRange -> All, PlotStyle -> Texture[GeoImage[pico, "ReliefMap", geoRange, GeoZoomLevel -> 12]]]discretizedGraphics = DiscretizeGraphics[listPlot]boxhight = hight + 1 / 2hight;
openBox = RegionUnion[DiscretizeRegion[#, MaxCellMeasure -> Infinity]& /@ {
Polygon[{{0, 0, boxhight}, {length, 0, boxhight}, {length, length, boxhight}, {0, length, boxhight}}],
Polygon[{{0, 0, 0}, {length, 0, 0}, {length, 0, boxhight}, {0, 0, boxhight}}],
Polygon[{{0, length, 0}, {length, length, 0}, {length, length, boxhight}, {0, length, boxhight}}],
Polygon[{{0, 0, 0}, {0, length, 0}, {0, length, boxhight}, {0, 0, boxhight}}],
Polygon[{{length, 0, 0}, {length, length, 0}, {length, length, boxhight}, {length, 0, boxhight}}]
}]regionUnion = RegionUnion@@Flatten[{openBox, discretizedGraphics}];mesh = ToElementMesh[regionUnion, MaxCellMeasure -> Infinity];mesh["Wireframe"]次に,メッシュの一部が削除された起伏図を表す可視化を生成する.このためには,各要素の中心を計算し,最終的な可視化に表示したい要素だけを選ぶ.
eleCoords = NDSolve`FEM`GetElementCoordinates[mesh["Coordinates"], Flatten[ElementIncidents[mesh["MeshElements"]], 1]];
com = (Plus@@@eleCoords[[All, 1 ;; 4]]) / 4;theseEle = Pick[Range[Length[com]], com, _ ? ((0.25#[[1]] + #[[2]]) > 850&)];
theseEle = Transpose[{ConstantArray[1, Length[theseEle]], theseEle}];meshVisual = mesh["Wireframe"[theseEle, "MeshElement" -> "MeshElements", "ElementMeshDirective" -> Directive[EdgeForm[LightGray], FaceForm[Gray]]]];
Show[meshVisual, listPlot]領域の近似の質
Line,Polygon等のグラフィックスプリミティブ,あるいはMeshRegionでは,ElementMeshに変換しても無損失である.例えば,RectangleのElementMesh表現は,Rectangle自身が領域を表すのと同じぐらい厳密あるいは不正確である.近似の質を推測する1つの方法は,可能な場合に,問題の領域範囲を比べることである.
NIntegrate[1, {x, y}∈Rectangle[]]Total[ToElementMesh[Rectangle[]]["MeshElementMeasure"], 2]これを例えばDiskと比べてみよう.どれぐらい細かくメッシュが作成され,要素のメッシュ次数がどれほど高くても,離散化は近似にすぎない.
π - Total[ToElementMesh[Disk[]]["MeshElementMeasure"], 2]mesh = ToElementMesh[Disk[]];
Max[1 - Total[(Join@@GetElementCoordinates[mesh["Coordinates"], Join@@ElementIncidents[mesh["PointElements"]]]) ^ 2, {2}]]π - Total[ToElementMesh[Disk[], "MeshOrder" -> 1]["MeshElementMeasure"], 2]Diskの場合,ElementMeshが曲線の境界メッシュ要素を持つことができるので,メッシュ次数が重要な役割を果たす.
π - Total[ToElementMesh[Disk[], "MaxBoundaryCellMeasure" -> 0.01, "BoundaryMeshGenerator" -> "Continuation" ]["MeshElementMeasure"], 2]また,境界の粒度は,どのくらいうまく領域が近似できるかということにおいて影響を与える.領域変換の全体的な確度は,AccuracyGoalを通して制御される.
π - Total[ToElementMesh[Disk[], AccuracyGoal -> 8, "MeshOrder" -> 1]["MeshElementMeasure"], 2]π - Total[ToElementMesh[Disk[], AccuracyGoal -> 8, "MeshOrder" -> 2]["MeshElementMeasure"], 2]内部では,まずNumericalRegionが作成される.ToBoundaryMeshがその数値領域に対して呼び出されてから,完全なメッシュが作成される.これが境界上の新しいノードを取り込むことがある.二次メッシュでは,追加で中間のノードが挿入される.その後,これらの新しい境界ノードの位置を改善する関数が呼び出される.
π - Total[ToElementMesh[Disk[], "ImproveBoundaryPosition" -> False]["MeshElementMeasure"], 2]NumericalRegionには,領域の厳密な記号表現,境界ElementMeshおよび完全なElementMeshを含むことができる.記号領域が使用でき,境界ノードの厳密な位置と境界上のより高次のノードを求めることができるので,境界ノードの位置は改善される.この過程については,Numerical Regionsのセクションでさらに詳しく説明する.
境界メッシュの作成のために,2つの境界メッシュ生成器が使用できる.1つ目はRegionPlotに基づき,"RegionPlot"と呼ばれる.これは,高速の境界近似を提供する.デフォルトの境界メッシュ生成器は,"Continuation"と呼ばれ,継続法に基づく.この境界メッシュ生成器は,いくぶんゆっくりではあるが,高い確度が達成できる.
Ω = ImplicitRegion[-1 + 2 x^2 ≤ y ≤ x^2, {{x, -1, 1}, {y, -1, 1}}]m = ToElementMesh[Ω, "BoundaryMeshGenerator" -> {"Continuation"}];
m["Wireframe"]4 / 3 - Total[m["MeshElementMeasure"], 2]m = ToElementMesh[Ω, "BoundaryMeshGenerator" -> {"RegionPlot"}];
m["Wireframe"]視覚的に調べると,尖点がうまく解像されていないことが分かる."RegionPlot"の近似は,サンプル点の数を増やす(これはPlotPoints)に似ている)ことによって改善されることがある.
m = ToElementMesh[Ω, "BoundaryMeshGenerator" -> {"RegionPlot", "SamplePoints" -> {51, 21}}];
m["Wireframe"]4 / 3 - Total[m["MeshElementMeasure"], 2]要素メッシュの質
領域がどれほど要素メッシュによってうまく近似されているかにかかわらず,メッシュの要素自体も,数値タスクの解に影響を与える.さまざまな数値アプリケーションについて,要素メッシュに対するさまざまな制約条件が影響してくる.有限要素法では,以下が大きな離散化エラーを引き起す[4].
メッシュの全体的な質は,
から
までの数として表現することができる.
は,与えられた質の推定で最高質,
は最低質を示す.負の質は,要素のインシデントが正しい次数ではないことを示すか,インシデントが自己交差する要素を作成することを示すか,そのどちらもを示すかする.
ElementMeshは,その質の概念を持つ.いったんElementMeshが作成されると,メッシュの質のクエリを行うことができる.
nx = 26;ny = 9;
coordinates = Flatten[ Table[{r Cos[θ], r Sin[θ]}, {r, 1., 2., 1 / (ny - 1)}, {θ, 0., 2 Pi / 3., (2 Pi / 3.) / (nx - 1)}], 1];
incidents = Flatten[Table[{j * nx + i, j * nx + i + 1, (j - 1) * nx + i + 1, (j - 1) * nx + i}, {i, 1, nx - 1}, {j, 1, ny - 1}], 1];mesh = ToElementMesh["Coordinates" -> coordinates, "MeshElements" -> {QuadElement[incidents]}]mesh["Wireframe"]Short[q = mesh["Quality"]]"Quality"計算の結果は,各メッシュ要素の品質推定をメッシュ要素でグループ化したものである.
Through[{Mean, StandardDeviation}[#]]& /@ qMinMax[q]Histogram[q, {0, 1, 0.1}]線分要素メッシュ
三角要素メッシュ
三角要素メッシュの質は,次の公式に従って計算される [1, 2].
ここでは,
は三角形の範囲であり,
は
番目の辺の長さである.
mesh = ToElementMesh["Coordinates" -> {{-1., -1.}, {1., -1.}, {0., Sqrt[3] - 1.}}, "MeshElements" -> {TriangleElement[{{1, 2, 3}}]}];
mesh["Quality"]クワッド要素メッシュ
ここでは,
はクワッドのiの範囲であり,
は
番目の辺の長さである.
mesh = ToElementMesh["Coordinates" -> {{0., 0.}, {1., 0.}, {1., 1.}, {0., 1.}}, "MeshElements" -> {QuadElement[{{1, 2, 3, 4}}]}];
mesh["Quality"]四面体要素メッシュ
四面体要素メッシュの質は,次の公式に従って計算される [1, 3].
ここでは,
は四面体の体積であり,
は
として実装される.
は
番目の辺の長さである.
mesh = ToElementMesh["Coordinates" -> {{1., 0., 0.}, {1., 1., 1.}, {0., 1., 0.}, {0., 0., 1.}}, "MeshElements" -> {TetrahedronElement[{{1, 2, 3, 4}}]}];
mesh["Quality"]六面体要素メッシュ
Hereここでは,
は六面体の体積であり,
は
番目の辺の長さである.
mesh = ToElementMesh["Coordinates" -> MeshElementBaseCoordinates[HexahedronElement, 1], "MeshElements" -> {HexahedronElement[{{1, 2, 3, 4, 5, 6, 7, 8}}]}];
mesh["Quality"]参照
[1] J. R. Shewchuk, "What Is a Good Linear Finite Element? Interpolation, Conditioning, Anisotropy, and Quality Measures (Preprint)," 2002, unpublished.
[2] R. P. Bhatia and K. L. Lawrence, "Two-Dimensional Finite Element Mesh Generation Based on Stripwise Automatic Triangulation," Computers and Structures, 36, 1990 pp. 309–319.
[3] V. N. Parthasarathy, C. M. Graichen, and A. Hathaway, "Fast Evaluation & Improvement of Tetrahedral 3-D Grid Quality," 1991 (manuscript).
[4] S.-W. Cheng, T. K. Dey, and J. R. Shewchuk, Delaunay Mesh Generation, Boca Raton: CRC Press, 2013.
低品質要素を可視化する
低品質の要素を含むメッシュが生成される場合には,これらの要素を可視化してから,その特定の範囲における質に対処することが必要になることがある.
gc = DiscretizeGraphics[ExampleData[{"Geometry3D", "SpaceShuttle"}, "GraphicsComplex"]];
mesh = ToElementMesh[gc]mesh["Wireframe"]q = mesh["Quality"];
{Min /@ q, Mean /@ q, StandardDeviation /@ q}Histogram[q, {0, 1, 0.1}]メッシュの平均的な質を改善することは可能であるが,悪い質の推定を持つ要素が,形状のコーナーにある可能性もある.形状が形に影響するので,必ずしもこれらの要素の質を改善できるとは限らない.しかしながら,質についてわきまえておくことは大切である.
mesh = ToElementMesh[gc, MeshQualityGoal -> 1]q = mesh["Quality"];
{Min /@ q, Mean /@ q, StandardDeviation /@ q}全体的な平均の質が改善されたことが分かる.しかし,メッシュにはより多くの要素が含まれることになり,数値解析に少し余計に計算時間がかかるようになった.
Histogram[q, {0, 1, 0.1}]mesh["Wireframe"]一定の閾値より低い要素を可視化するためには,低品質の要素がメッシュ要素の質のリストから選ばれ,可視化される.
pos = Position[q, _ ? (# ≤ 0.1&)];
mesh["Wireframe"[pos, "MeshElement" -> "MeshElements", Boxed -> False]]pos = Position[q, _ ? (# ≤ 0.1&)];
Show[
mesh["Wireframe"["ElementMeshDirective" -> Directive[EdgeForm[GrayLevel[.6]], FaceForm[]]]],
mesh["Wireframe"[pos, "MeshElement" -> "MeshElements", "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]
, Boxed -> False]部分領域を持つ要素メッシュ
偏微分方程式が,複数の物質からなる領域とインタラクトするということがよくある.偏微分方程式の解は,メッシュ要素が内部材料の境界と交差しない場合の方が高品質のものになる.このことを示すために,変数拡散係数を持つ偏微分方程式をもう一度考慮し,2つの領域について解く(「有限要素法で偏微分方程式を解く」)を参照のこと.).1つの領域は内部境界を尊重するが,もう一つの領域はしない.
sh = 0.2;sh2 = 0.02;sw = 0.3;
coordinates = {{0., 0.}, {1., 0.}, {1., sh}, {1., 1.}, {0., 1.}, {0., sh + sh2}, {sw, sh + sh2}, {sw, sh}, {0., sh}};
el1 = LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 5}, {5, 6}, {6, 7}, {7, 8}, {8, 9}, {9, 1}}];
bMesh1 = ToBoundaryMesh["Coordinates" -> coordinates, "BoundaryElements" -> {el1, LineElement[{{3, 8}}]}];
bMesh2 = ToBoundaryMesh["Coordinates" -> coordinates, "BoundaryElements" -> {el1}];
GraphicsRow[{bMesh1["Wireframe"], bMesh2["Wireframe"]}]ϵr = If[y ≤ sh, {{11.7, 0.}, {0., 11.7}}, {{1., 0}, {0., 1.}}]mesh1 = ToElementMesh[bMesh1];
mesh2 = ToElementMesh[bMesh2];
{mesh1["Wireframe"], mesh2["Wireframe"]}op = Inactive[Div][-ϵr.Inactive[Grad][u[x, y], {x, y}], {x, y}] - 10 ^ -8. / 8.86*^-12;
Subscript[Γ, D] = {DirichletCondition[u[x, y] == 0, x == 1 || y == 1 || y == 0],
DirichletCondition[u[x, y] == 10 ^ 3, 0 < x ≤ sw && sh ≤ y ≤ sh + sh2]};ufun1 = NDSolveValue[{op == 0, Subscript[Γ, D]}, u, {x, y}∈mesh1];ufun2 = NDSolveValue[{op == 0, Subscript[Γ, D]}, u, {x, y}∈mesh2];{
Show[ContourPlot[ufun1[x, y], {x, y}∈mesh1, ColorFunction -> "TemperatureMap", AspectRatio -> Automatic],
bMesh1["Wireframe"]],
Show[ContourPlot[ufun2[x, y], {x, y}∈mesh2, ColorFunction -> "TemperatureMap", AspectRatio -> Automatic],
bMesh2["Wireframe"]]
}Plot3D[ufun1[x, y] - ufun2[x, y], {x, y}∈mesh1, ColorFunction -> "TemperatureMap", PlotRange -> All, Boxed -> False]次の例は,部分領域を持つ円上の領域を示す.部分領域の中には,穴がある.
annulus[x_, y_] := (9 / 10) ^ 2 ≤ x ^ 2 + y ^ 2 ≤ 1 ^ 2
holes[{x0_, y0_}, r_] := ((x + x0) ^ 2 + (y + y0) ^ 2 ≤ (r) ^ 2)
crds = {{-1 / 2, 0}, {1 / 2, 0}, {0, -1 / 2}, {0, 1 / 2}, {2 / 5, 2 / 5}, {-2 / 5, -2 / 5}, {2 / 5, -2 / 5}, {-2 / 5, 2 / 5}};
sd = Or@@(holes[#, 1 / 8]& /@ crds);Show[
RegionPlot[annulus[x, y], {x, -1, 1}, {y, -1, 1}, PlotStyle -> LightOrange],
RegionPlot[x ^ 2 + y ^ 2 < (9 / 10) ^ 2 && !sd, {x, -1, 1}, {y, -1, 1}]
, Frame -> False]内側の境界を守る領域からメッシュを簡単に作るために,領域の設定はすべての内側の境界が作成されるような方法で行われる.
Ω2 = ImplicitRegion[Or[annulus[x, y], sd], {x, y}];
ToElementMesh[Ω2]["Wireframe"]次に,部分領域内の穴の部分と穴ではない部分が,領域の穴を明示的に指定することによって,反転される.
ToElementMesh[Ω2, "RegionHoles" -> crds]["Wireframe"]その後,例えば部分領域の1つをさらに細かくするということも可能である.
ToElementMesh[Ω2, "RegionHoles" -> crds, "RegionMarker" -> {{{0, 0}, 1, 0.01}, {{19 / 20, 0}, 2, 0.001}}]["Wireframe"]一点気を付けなければならないのは,DirichletConditionのような境界条件をTrueの領域述部に適用することである. Trueを述部として指定することによって,内部境界を含めたすべての境界にDirichletConditionが適用される.これは必ずしも意図される結果ではないかもしれない.その場合には,ディリクレ条件が適用されるべき境界部分を明示的に指定した方がよい.DirichletConditionの「考えられる問題」セクションも参照されたい.
部分領域を持つメッシュを作成するためには,穴のある領域を作成することを完全に避けるとよい場合もある.これは,"RegionHoles"オプションをNoneに設定することで指定できる.
mesh = ToElementMesh[Ω2, "RegionHoles" -> None, "RegionMarker" -> Join[
MapThread[{#1, #2, 0.001}&, {crds, Range[Length[crds]]}], {{{0, 0}, Length[crds] + 1, 0.01}, {{19 / 20, 0}, Length[crds] + 2, 0.002}}]];
temp = Most[Range[0, 1, 1 / (Length[crds] + 2)]];
colors = ColorData["BrightBands"][#]& /@ temp;
mesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]三次元において材料領域を作成する方法についての追加情報は,OpenCascadeLinkに記載されている.
複数の材料のシミュレーションについての実際の例は,例えば有限要素法で偏微分方程式を解くやPDEModelsの概要ページに記載されている.分野特定のモノグラフの多くにも複数の材料のシミュレーションについてのセクションがある.あるいは,物理の特定分野,例えば,物質輸送に進むと,その分野の特定のガイドページに応用例や複数の材料の領域等の特定の特徴をリストするモデルの概要セクションがある.
マーカー
マーカーは,メッシュの要素に関連付けられている正の整数番号である.メッシュにおけるマーカーの主な目的は,偏微分方程式指定において,十歳の座標から(境界)述部を切り離すことにある.つまり,偏微分方程式が指定されたときに,偏微分方程式が座標から独立するように指定することができる.この機能は,シミュレーションの領域が変化する可能性があるが,偏微分方程式は変化しないといういう場合に便利である.
次のセクションでは,複数の材料のモデリングにメッシュ要素マーカーを使うことについて説明し,例を挙げる.
- 境界マーカーは,領域境界上のNeumannValueを指定するのに便利であり,"BoundaryElements"の一種である.
- 点マーカーは,境域境界上のDirichletConditionを指定するのに便利であり,"PointElements"の一種である.
mesh = ToElementMesh["Coordinates" -> {{0, 0}, {1, 0}, {1, 1}, {0, 1}}, "MeshElements" -> {TriangleElement[{{1, 2, 3}, {3, 4, 1}}, {10, 20}]}]この場合は,三要素には各要素について1つずつ,2つのマーカーが含まれる.
mesh["Wireframe"["MeshElementMarkerStyle" -> Blue]]物質マーカーの利用を示すために,変数拡散係数を持つ偏微分方程式をもう一度考慮し,解く.(「有限要素法で偏微分方程式を解く」)を参照のこと).
sh = 0.2;sh2 = 0.02;sw = 0.3;
bmesh = ToBoundaryMesh["Coordinates" -> {{0., 0.}, {1., 0.}, {1., sh}, {1., 1.}, {0., 1.}, {0., sh + sh2}, {sw, sh + sh2}, {sw, sh}, {0., sh}}, "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 5}, {5, 6}, {6, 7}, {7, 8}, {8, 9}, {9, 1}, {3, 8}}]}];
bmesh["Wireframe"]ToElementMeshの"RegionMarker"オプションで領域マーカーを指定することが可能である.このためには,領域内の座標と整数マーカーが与えられなければならない.オプションとして,追加の最大セル測定値を指定し,部分領域を精緻化することもできる.
mesh = ToElementMesh[bmesh, "RegionMarker" -> {{{0.1, sh / 2}, 10, 0.001}, {{0.1, sh * 2}, 20}}];
mesh["Wireframe"["MeshElementStyle" -> {Directive[FaceForm[Green]], Directive[FaceForm[Red]]}]]mesh["MeshElementMarkerUnion"]偏微分方程式の係数は,述部のElementMarkerにアクセスすることによって指定できるようになった.
ϵr = If[ElementMarker == 10, {{11.7, 0.}, {0., 11.7}}, {{1., 0}, {0., 1.}}]op = Inactive[Div][-ϵr.Inactive[Grad][u[x, y], {x, y}], {x, y}] - 10 ^ -8. / 8.86*^-12;
Subscript[Γ, D] = {DirichletCondition[u[x, y] == 0, x == 1 || y == 1 || y == 0],
DirichletCondition[u[x, y] == 10 ^ 3, 0 < x ≤ sw && sh ≤ y ≤ sh + sh2]};この場合は非アクティブな偏微分方程式が使われているが,
のような係数を指定する過程は,HeatTransferPDEComponent等の偏微分方程式の演算子が使われる場合と同じであることに注意する.後者の場合,"MassDensity"等のパラメータには以下のような条件文も使える.
"MassDensity"->If[ElementMarker10,{{11.7,0.},{0.,11.7}},{{1.,0},{0.,1.}}].
これについての詳細は,HeatTransferのモノグラフ等,物理の関連分野のモノグラフに記載されている.
ufun = NDSolveValue[{op == 0, Subscript[Γ, D]}, u, {x, y}∈mesh];ContourPlot[ufun[x, y], {x, y}∈mesh, ColorFunction -> "TemperatureMap", AspectRatio -> Automatic]通常,要素メッシュ内の境界と部分領域は,述部のある形式によって与えられる.代りに,境界メッシュ要素は,境界条件指定で使用されるマーカーを持つこともできる.
sh = 0.2;sh2 = 0.02;sw = 0.3;
bmesh = ToBoundaryMesh["Coordinates" -> {{0., 0.}, {1., 0.}, {1., sh}, {1., 1.}, {0., 1.}, {0., sh + sh2}, {sw, sh + sh2}, {sw, sh}, {0., sh}}, "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 5}, {5, 6}, {6, 7}, {7, 8}, {8, 9}, {9, 1}, {3, 8}}, {1, 1, 1, 1, 2, 3, 3, 3, 2, 4}]}];
bmesh["Wireframe"["MeshElementMarkerStyle" -> Blue, "MeshElementStyle" -> {Blue, Green, Red, Yellow}]]mesh = ToElementMesh[bmesh, "RegionMarker" -> {{{0.1, sh / 2}, 10, 0.001}, {{0.1, sh * 2}, 20}}];
mesh["Wireframe"["MeshElementStyle" -> {FaceForm[Green], FaceForm[Red]}]]mesh["Wireframe"["MeshElement" -> "BoundaryElements", "MeshElementMarkerStyle" -> Blue, "MeshElementStyle" -> {Blue, Green, Red, Yellow}]]mesh["BoundaryElementMarkerUnion"]mesh["PointElementMarkerUnion"]境界と点要素の重要性は,偏微分方程式を解いた場合に,境界条件が適用される場所であるということである.一般化されたノイマン(Neumann)境界条件は,境界要素上で積分することによって適用される.ディリクレ境界条件は,点要素で適用される.
新たに挿入された点とセグメントは,正しい境界マーカーを持つことに注意する.
ϵr = If[ElementMarker == 10, {{11.7, 0.}, {0., 11.7}}, {{1., 0}, {0., 1.}}]op = Inactive[Div][-ϵr.Inactive[Grad][u[x, y], {x, y}], {x, y}] - 10 ^ -8. / 8.86*^-12Subscript[Γ, D] = {DirichletCondition[u[x, y] == 0, ElementMarker == 1],
DirichletCondition[u[x, y] == 10 ^ 3, ElementMarker == 3]};ufun = NDSolveValue[{op == 0, Subscript[Γ, D]}, u, {x, y}∈mesh]ContourPlot[ufun[x, y], {x, y}∈mesh, ColorFunction -> "TemperatureMap", AspectRatio -> Automatic]これで,偏微分方程式の係数と境界条件は,座標から独立したものになり,新しい形状は,偏微分方程式や境界条件を変更しなくてもすぐに配置できるようになったことに注意する.
まとめると,要素の集合は,形式 eltype[incidents,markers] で与えられる.incidents は,各要素の座標のリストのリストである.markers は,incidents と同じ長さのリストであり,物質特性のいくつかが異なっている,あるいは別の境界条件が適用される可能性のある領域や境界のさまざまな部分を識別するのに使える.指定の要素について評価する際には,ElementMarkerは事実上その要素のマーカーで置き換えられる.
"MeshElements","BoundaryElements","PointElements"の要素グループにはマーカーが含まれることがある.それぞれのグループでマーカーは異なる目的に使われる."MeshElements"の要素グループのマーカーは,偏微分方程式の係数に使える."BoundaryElements"の要素グループのマーカーは,NeumannValueまたはPeriodicBoundaryConditionで使え,"PointElements"の要素グループのマーカーは,DirichletConditionで使える.
方程式あるいは境界条件でマーカーを指定する際には,マーカーは
等の方程式として指定しなければならない.Several ElementMarker式のいくつかは,AndやOr等のブール述語と組み合せることができる.特定のマーカー以外のすべての境界において境界条件が有効であると表現するためには,系統的論述
が使える.ElementMarkerを指定するその他の方法は,現行バージョンのNDSolveではサポートされていない.
シミュレーションでElementMarkerを使用する過程を表すために,"MeshElements","BoundaryElements","PointElements"の中にElementMarkerを持つ簡単なメッシュを作成する.
mesh = ToElementMesh["Coordinates" -> {{1.293, 0.228}, {1., 0.}, {0.94, 0.342}, {1.293, 0.}, {1.215, 0.442}, {2., 0.}, {1.879, 0.684}}, "MeshElements" -> {TriangleElement[{{1, 3, 2}, {1, 2, 4}, {1, 4, 6}, {1, 6, 7}, {1, 7, 5}, {1, 5, 3}}, {66, 66, 66, 44, 44, 44}]},
"BoundaryElements" -> {LineElement[{{3, 2}, {1, 3}, {2, 4}, {4, 6}, {6, 1}, {6, 7}, {7, 5}, {5, 3}}, {11, 22, 33, 33, 22, 44, 55, 55}]},
"PointElements" -> {PointElement[{{1}, {2}, {3}, {4}, {5}, {6}, {7}}, {1, 2, 2, 3, 3, 4, 4}]}]Show[mesh["Wireframe"], mesh["Wireframe"["MeshElement" -> "PointElements", "MeshElementIDStyle" -> Brown]]]Show[
mesh["Wireframe"],
mesh["Wireframe"["MeshElementMarkerStyle" -> Blue]],
mesh["Wireframe"["MeshElement" -> "BoundaryElements", "MeshElementMarkerStyle" -> Red]],
mesh["Wireframe"["MeshElement" -> "PointElements", "MeshElementMarkerStyle" -> Brown]]
]ノード6からノード7までの境界辺は has a marker of 44のマーカーを持ち,ノード{1,6,7}からなるメッシュ要素も44のマーカーの集合を持つことに注意されたい.
例として,"MeshElements"のマーカーに依存する右辺を持つポワソンタイプの方程式を使う.メッシュ要素が44のマーカーを持つ場合には,右辺の値として10が使われる.その他の場合には1が使われる.また
のNeumannValueは44に設定されたマーカーを持つすべての辺で設定される.DirichletConditionは,1のElementMarkerを持つ点に設定される.
ufun = NDSolveValue[{Laplacian[u[x, y], {x, y}] == If[ElementMarker == 44, 10, 1] + NeumannValue[-1, ElementMarker == 44], DirichletCondition[u[x, y] == 0, ElementMarker == 1]}, u, {x, y}∈mesh]Plot3D[ufun[x, y], {x, y}∈mesh]メッシュ要素と境界要素がどちらも数字44のマーカーを使っても問題ない.If文のような係数で使われるElementMarkerは, "MeshElements"にあるマーカーに動作する.NeumannValueのElementMarkerへの参照は,"BoundaryElements"のマーカーを,DirichletConditionのElementMarkerは"PointElements"のマーカーを参照する.
形状が複雑であることもある.この場合には,述部で要素メッシュの部分を識別することは,必ずしも容易なことではない.
b = {-y, 1 / 25 - (-(3 / 2) + x)^2 - y^2, 1 - x^2 - y^2, -4 + x^2 + y^2, -(x - 2 y) / (√5)};
ContourPlot[Evaluate[(# == 0& /@ b)], {x, 0.8, 2.2}, {y, -0.2, 1.}, Contours -> {0}]bmesh = ToBoundaryMesh[ImplicitRegion[And@@(# <= 0& /@ b), {{x, 0.8, 2.2}, {y, -0.2, 1.}}], "MaxBoundaryCellMeasure" -> 0.2]bmesh["Wireframe"]デフォルトで,境界メッシュ生成の際には,隣接する境界の自動のグループ分けが計算される.自動のグループ分けは,境界の法線を調べ,それらを互いに比べることによって行われる.2つの法線ベクトルのドット積の絶対値が閾値を超える(つまり,ベクトルは同じ方向であることを意味する)場合には,隣接するセグメントが連続していると見なされる.
groups = bmesh["BoundaryElementMarkerUnion"]temp = Most[Range[0, 1, 1 / (Length[groups])]];
colors = ColorData["BrightBands"][#]& /@ tempbmesh["Wireframe"["MeshElementStyle" -> (Directive[#]& /@ colors)]]bn = bmesh["BoundaryNormals"];
mean = Mean /@ GetElementCoordinates[bmesh["Coordinates"], #]& /@ ElementIncidents[bmesh["BoundaryElements"]];
Show[
bmesh["Wireframe"],
Graphics[{DarkGreen, MapThread[Arrow[{#1, #2}]&, {Join@@mean, Join@@(bn / 14 + mean)}]}]]一旦自動の境界要素マーカーが導入されると,点の要素マーカーがそこから作られる.そのために,ノードは接続された境界要素マーカーを調べる.これらのマーカーが同じである場合には,そのマーカーが点の要素マーカーとして使われる.コーナーでは,ノードは異なるマーカーを持つ2つの境界部分に属する.その場合には,新しいマーカー番号が使われる.前のバージョンでは,接続されたノードとそのマーカーの間の選択はランダムであった.
bmesh["Wireframe"["MeshElement" -> "PointElements", "MeshElementMarkerStyle" -> Red]]これは,左の円状の部分におけるDirichletConditionが以下のようになることを意味する.
DirichletCondition[u[x, y] == value, ElementMarker == 12 | ElementMarker == 6 | ElementMarker == 8]ディリクレ条件は,コーナーのノードが境界条件の一部として考慮されるべきかによる.
コーナーノードにおけるマーカー(どのマーカーでも同じことが言えるが)の値を変更する方法の1つとして,ノードの座標とその座標で使われるマーカーを指定することによって,ユーティリティ関数ElementMeshResetPointElementMarkerを利用する方法がある.
bmesh2 = ElementMeshResetPointElementMarker[bmesh, {{2 / Sqrt[5], 1 / Sqrt[5]}, {1, 0}, {13 / 10, 0}, {2, 0}}, {6, 6, 2, 4}]bmesh2["Wireframe"["MeshElement" -> "PointElements", "MeshElementMarkerStyle" -> Red]]The上からの左の円状部分におけるDirichletConditionは以下のようになる.
DirichletCondition[u[x, y] == value, ElementMarker == 6]もう一つの方法として,左側のノードがすべてマーカー6に属し,右側がマーカー4に属する場合には,"PointMakerFunction"を使ってメッシュを生成することもできる.これは
座標に基づいてマーカーを変更する.その他の場合はすべて,自動生成された点マーカーが使われる.
pointMarkerFunction = Compile[{{coords, _Real, 2}, {pMarker, _Integer, 1}},
MapThread[
Block[{x = #1[[1]], y = #1[[2]], autoMarker = #2},
Which[
x ≤ 1, 6,
x ≥ 2, 4,
True, autoMarker]
]&, {coords, pMarker}]];bmesh = ToBoundaryMesh[ImplicitRegion[And@@(# <= 0& /@ b), {{x, 0.8, 2.2}, {y, -0.2, 1.}}], "MaxBoundaryCellMeasure" -> 0.2, "PointMarkerFunction" -> pointMarkerFunction]bmesh["Wireframe"["MeshElement" -> "PointElements", "MeshElementMarkerStyle" -> Red]]完全なメッシュの生成中に,これらの境界マーカーは伝播され,上で説明したように,後程NDSolve内から使うことができる.
境界メッシュ内にマーカーを任意に置くことを可能にする関数を書くこともできる.境界マーカーの配置には2つの関数が使用できる.1つは境界点に作用するもので,もう一つは,2Dでは境界線に,そして3Dでは境界面に作用するものである.
pointMarkerFunction = Compile[{{coords, _Real, 2}, {pMarker, _Integer, 1}}, Block[{x = #[[1]], y = #[[2]], epsilon},
epsilon = 10 ^ -5.;
Which[
Abs[1 / 25 - (-3 / 2 + x) ^ 2 - y ^ 2] ≤ epsilon, 2,
Abs[1 - x ^ 2 - y ^ 2] ≤ epsilon, 3,
Abs[-4 + x ^ 2 + y ^ 2] ≤ epsilon, 4,
Abs[(-x + 2 * y) / Sqrt[5]] ≤ epsilon, 1,
Abs[-y] ≤ epsilon, 1,
True, 0
]]& /@ coords];pm = pointMarkerFunction[bmesh["Coordinates"], Flatten[ElementMarkers[bmesh["PointElements"]]]]Graphics[{PointSize[0.002], Point[{4 / (√5), 2 / (√5)}], MapThread[ Text[ToString[#1], #2]&, {pm, bmesh["Coordinates"]}]}]上のWhich文の順序が重要であることに注意する.座標が境界のいくつかの上にある場合には,マッチする最初の述部が返される.
2つ目の関数は,境界の辺に働くように書くことができる.この関数は,境界要素の座標と,点関数のすでに計算されたマーカーを得る.
boundaryMarkerFunction = Compile[{{boundaryElementCoords, _Real, 3}, {boundaryElementPointMarkres, _Integer, 2}}, Which[
Union[#] == {2}, 2,
Union[#] == {3}, 3,
Union[#] == {4}, 4,
True, 1 ]& /@ boundaryElementPointMarkres];bmesh = ToBoundaryMesh[ImplicitRegion[And@@(# <= 0& /@ b), {{x, 0.8, 2.2}, {y, -0.2, 1.}}], "PointMarkerFunction" -> pointMarkerFunction, "BoundaryMarkerFunction" -> boundaryMarkerFunction,
"MaxBoundaryCellMeasure" -> 0.1]bmesh["BoundaryElements"]bmesh["Wireframe"["MeshElementMarkerStyle" -> Red, "MeshElementStyle" -> {Blue, Yellow, Green, Red}]]ElementIncidents[bmesh["BoundaryElements"]]ElementMarkers[bmesh["BoundaryElements"]]bmesh["PointElements"]bmesh["Wireframe"["MeshElement" -> "PointElements", "MeshElementMarkerStyle" -> Red, "MeshElementStyle" -> {Blue, Yellow, Green, Red}]]ElementIncidents[bmesh["PointElements"]]ElementMarkers[bmesh["PointElements"]]mesh = ToElementMesh[bmesh, "MaxCellMeasure" -> 0.1]mesh["BoundaryElementMarkerUnion"]Show[mesh["Wireframe"], mesh["Wireframe"["MeshElement" -> "BoundaryElements", "MeshElementStyle" -> {Blue, Yellow, Green, Red}]]]mesh["PointElementMarkerUnion"]点マーカーの和集合は,メッシュ要素のスタイル指示子をいくつ与えなければならないかを示す.
Show[mesh["Wireframe"], mesh["Wireframe"["MeshElement" -> "PointElements", "MeshElementStyle" -> {Blue, Yellow, Green, Red}]]]辺と点におけるマーカーの値は,完全なメッシュが生成された後,伝播したことに注意する.
ElementIncidents[mesh["MeshElements"]]//ShortElementMarkers[mesh["MeshElements"]]//Short"RegionMarker"オプションを定義することによって,メッシュ要素マーカーを加えることができる.
別の方法として,完全なメッシュを直接生成することも可能である.
mesh = ToElementMesh[ImplicitRegion[And@@(# <= 0& /@ b), {{x, 0.8, 2.2}, {y, -0.2, 1.}}], "PointMarkerFunction" -> pointMarkerFunction, "BoundaryMarkerFunction" -> boundaryMarkerFunction, "RegionHoles" -> None]3Dにおける平面も部分分割してその部分分割の辺を維持することが可能である.しかし,そうするためには,小さい面に区別できるマーカーを指定し,小さい面の部分分割が完全なメッシュに存在することを示す必要がある.
bmesh = ToBoundaryMesh["Coordinates" -> {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 0, 1 / 2}, {1, 0, 1 / 2}, {1, 1, 1 / 2}, {0, 1, 1 / 2}, {1 / 2, 1 / 2, 1 / 2}}, "BoundaryElements" -> {QuadElement[{{1, 2, 3, 4}, {1, 2, 6, 5}, {2, 3, 7, 6}, {3, 4, 8, 7}, {4, 1, 5, 8}}], TriangleElement[{{5, 6, 9}, {6, 7, 9}, {7, 8, 9}, {8, 5, 9}}]}];
bmesh["Wireframe"]境界要素にマーカーが指定されていない場合には,平面の小さい面は合成される.
ToElementMesh[bmesh]["Wireframe"[PlotRange -> {All, All, {1 / 4, 3 / 3}}]]完全なメッシュを生成すると,上の曲面の分割が失われてしまう.
曲面の部分分割を維持するためには,要素にマーカーを加える必要がある.マーカーが他から区別できるものである場合には,そのマーカーを分割する辺は維持される.
bmesh = ToBoundaryMesh["Coordinates" -> {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 0, 1 / 2}, {1, 0, 1 / 2}, {1, 1, 1 / 2}, {0, 1, 1 / 2}, {1 / 2, 1 / 2, 1 / 2}}, "BoundaryElements" -> {QuadElement[{{1, 2, 3, 4}, {1, 2, 6, 5}, {2, 3, 7, 6}, {3, 4, 8, 7}, {4, 1, 5, 8}}], TriangleElement[{{5, 6, 9}, {6, 7, 9}, {7, 8, 9}, {8, 5, 9}}, {1, 2, 3, 3}]}];
bmesh["Wireframe"]ToElementMesh[bmesh]["Wireframe"[PlotRange -> {All, All, {1 / 4, 3 / 3}}, "MeshElement" -> "BoundaryElements", "MeshElementStyle" -> (FaceForm /@ ColorData["SunsetColors"] /@ {1 / 5, 2 / 5, 3 / 5, 4 / 5, 1})]]この場合,上の曲面の4つの分割面のうち,2つは同じマーカー3を使っているため合成される.マーカー1とマーカー2の曲面要素は完全なメッシュでも維持されている.
ときには,メッシュ要素をそのマーカーで分割すると便利なことがある.
splitByMarker = Flatten[MeshElementSplitByMarker[bmesh["BoundaryElements"]]]Manipulate[Show[bmesh["Wireframe"], Graphics3D[Polygon[GetElementCoordinates[bmesh["Coordinates"], ElementIncidents[group]]]]], {group, splitByMarker, Slider}, SaveDefinitions -> True]点マーカーの計算は,オプション"PointMarkers""BoundaryDeduced"をToElementMeshに指定することによって,よりよく制御することができる.このメソッドを指定すると,新しい点マーカーのIDを作成する際により厳密なアプローチを使うことができる.
reg = RegionUnion[Cuboid[{1, 2, 0}, {3, 5 / 2, 1}], Cylinder[{{3, 1 / 2, 1 / 2}, {3, 3, 1 / 2}}, 1 / 4]];
(mesh = ToElementMesh[reg, "PointMarkers" -> "BoundaryDeduced"])["Wireframe"]bIDs = mesh["BoundaryElementMarkerUnion"];
Manipulate[Show[
mesh["Wireframe"["MeshElementStyle" ->
Directive[Opacity[0.2], FaceForm[LightBlue], EdgeForm[]]]],
mesh["Edgeframe"],
mesh["Wireframe"[ElementMarker == bIDs[[id]], "MeshElementStyle" -> Directive[FaceForm[LightGray], EdgeForm[]]]]
], {{id, 1, "ElementMarker ID"}, 1, Length[bIDs], 1, Appearance -> "Open"}, SaveDefinitions -> True]position = SelectPointMarkerFromBoundaryMarker[mesh, {6, 8}]position = SelectPointMarkerFromBoundaryMarker[mesh, 6 | 8]Show[
mesh["Wireframe"["MeshElementStyle" -> Directive[Opacity[0.2], FaceForm[LightBlue], EdgeForm[]]]],
mesh["Edgeframe"],
mesh["Wireframe"[position, "MeshElement" -> "PointElements", "MeshElementStyle" -> Directive[Red]]]
]すべての曲面に繋がっている点のみを考慮するような交点を求める.IDが7の境界マーカーは,大きな円柱の曲面である.
position = SelectPointMarkerFromBoundaryMarker[mesh, {6 | 8, 7}, Complement[#1, #2]&]Show[
mesh["Wireframe"["MeshElementStyle" -> Directive[Opacity[0.2], FaceForm[LightBlue], EdgeForm[]]]],
mesh["Edgeframe"],
mesh["Wireframe"[position, "MeshElement" -> "PointElements", "MeshElementStyle" -> Directive[Red]]]
]上の例では,円柱と選択した曲面を繋ぐ辺上には,点マーカーがないことに注意する.
制約条件が満たされない場合には,PointMarkerFromBoundaryMarkerはFalseを返し,DirichletConditionで使われる場合には,「FalseがTrueである場所は境界上にないので,DirichletCondition[u1,False]は事実上無視されます.」というメッセージが返される.
SelectPointMarkerFromBoundaryMarker[mesh, 1000]3Dで複数の材料領域を作成し表示する例については,OpenCascadeLinkのドキュメントを参照されたい.
マーカー付きの式を評価する
ElementMarkerに依存する式を評価したい場合がある.これは関数EvaluateOnElementMeshを使って行うことができる.
mesh = ToElementMesh[Annulus[], "RegionMarker" -> {{{0, 0}, 1}, {{3 / 4, 0}, 2}}, "RegionHoles" -> None]mesh["Wireframe"["MeshElementStyle" -> {FaceForm[Green], FaceForm[Red]}]]expr = If[ElementMarker == 1, 1, 2]fun = EvaluateOnElementMesh[{x, y}, expr, mesh]メッシュ上で式を評価すると,他の関数と同じように使える補間関数が返される.
NIntegrate[fun[x, y], {x, y}∈mesh]2 * Integrate[1, {x, y}∈Annulus[]] + Integrate[1, {x, y}∈Disk[{0, 0}, 1 / 2]]//Nfun[2, 0]この振舞いは,ElementMeshInterpolationに説明されるように変更することができる.
fun = EvaluateOnElementMesh[{x, y}, expr, mesh, "ExtrapolationHandler" -> {Function[Indeterminate], "WarningMessage" -> False}]fun[2, 0]EvaluateOnElementMeshの使用は,非連続な係数のオーバーシュートまたはアンダーシュートの問題のセクションで説明されるように,オーバーシュートおよびアンダーシュートに繋がることがある.
単位
ときには,メッシュのスケールを再設定しなければならないことがある.例えば,インポートされたメッシュが正しいスケールではない場合や,許容範囲の問題を避けるためにメッシュが別のスケールで生成され,後のステップで調整しなければならない場合がある.
bmesh = ToBoundaryMesh[Annulus[{1, 0}, {1, 2}]];bmesh["LengthUnit"]SetLengthUnit[bmesh, "Milimeters"]bmesh2 = ElementMeshCoordinateRescale[bmesh, "Meters"]しかしながら,後に行う有限要素解析で,結果の数値のよりよい安定性のために,偏微分方程式係数を再設定したほうがよいこともある.
他の関数での要素メッシュ
領域メンバーシップのテスト
ある点について,その点が領域にあるかどうかは,RegionMemberを使ってテストすることができる.
Ω = ImplicitRegion[ 25 x y ≥ 1 && x^2 + y^2 ≤ 1, {{x, 0, 1}, {y, 0, 1}}];
RegionMember[Ω, {.5, .6}]いったん領域が指定されると,それを有限要素法で使うためには,領域をほぼ網羅する要素にそれを部分分割する,つまりメッシュにしなければならない.メッシュがすでに利用可能である場合には,ElementMeshRegionMemberが使える.
mesh = ToElementMesh[Ω];
ElementMeshRegionMember[mesh, {.5, .6}]点の集合がメッシュにされた領域内にあるか,それとも領域外にあるかをテストする.
pts = RandomReal[1, {10000, 2}];
members = ElementMeshRegionMember[mesh, pts];
inside = Pick[pts, members, True];
outside = Pick[pts, members, False];
Graphics[{Point[outside], {Blue, Point[inside]}}]領域の境界に近い点については,誤った結果が出される可能性がある.これは,二次近時が連続境界を正しく表すのに十分ではない場合でさえもあり得ることである.
pts = RandomReal[1, {1000000, 2}];
members = ElementMeshRegionMember[mesh, pts];
exact = RegionMember[Ω][pts];
possibleBoundaryPos = Position[MapThread[Xnor, {exact, members}], False]
f[{x_, y_}] := (25 x y - 1)(x ^ 2 + y ^ 2 - 1)
f /@ Extract[pts, possibleBoundaryPos]Show[
mesh["Wireframe"],
Graphics[{Red, PointSize[0.03], Point[Extract[pts, possibleBoundaryPos]]}]
]これについて一つ可能なことは,領域の境界を精緻化することである.
mesh = ToElementMesh[Ω, "MaxBoundaryCellMeasure" -> 0.01]mesh["Wireframe"]members = ElementMeshRegionMember[mesh, pts];
possibleBoundaryPos = Position[MapThread[Xnor, {exact, members}], False]
f[{x_, y_}] := (25 x y - 1)(x ^ 2 + y ^ 2 - 1)
f /@ Extract[pts, possibleBoundaryPos]補間
NDSolveが有限要素法を通して偏微分方程式の解を計算する場合には,返されたInterpolatingFunctionには要素メッシュが含まれる.
ufun = NDSolveValue[{Laplacian[u[x, y], {x, y}] == 20, DirichletCondition[u[x, y] == 0, True]}, u, {x, y}∈Disk[]];ufun["ElementMesh"]ElementMeshからInterpolatingFunctionを構築することも可能である.
m = ToElementMesh[Disk[]]data = Sqrt[Total[m["Coordinates"] ^ 2, {2}]];efun = ElementMeshInterpolation[{m}, data]Plot3D[efun[x, y], {x, y}∈m]補間関数は,偏微分方程式の係数として使うことができる.補間関数のメッシュが偏微分方程式の解の過程のメッシュと同じである場合には,解の過程はより速くなる.
mesh = ToElementMesh[Rectangle[], MaxCellMeasure -> 0.0005];
fun = ElementMeshInterpolation[mesh, Sqrt[Total[mesh["Coordinates"] ^ 2, {2}]]];
RepeatedTiming[sol = NDSolveValue[{-Laplacian[u[x, y], {x, y}] + fun[x, y] * u[x, y] == 1, DirichletCondition[u[x, y] == 0, x == 0]}, u, Element[{x, y}, mesh]];]mesh = ToElementMesh[Rectangle[], MaxCellMeasure -> 0.00025];
fun = ElementMeshInterpolation[mesh, Sqrt[Total[mesh["Coordinates"] ^ 2, {2}]]];
mesh2 = ToElementMesh[Rectangle[], MaxCellMeasure -> 0.0005]
RepeatedTiming[sol = NDSolveValue[{-Laplacian[u[x, y], {x, y}] + fun[x, y] * u[x, y] == 1, DirichletCondition[u[x, y] == 0, x == 0]}, u, Element[{x, y}, mesh2]];]