異方性メッシュ族をPythonで再構成する：6つのmesh familyの節点座標・要素接点・境界辺接点表からgood/badを読む

Introduction

このページでは，Interpolation error analysis using a new geometric parameter Section 9 で与えられている異方性メッシュ族をコードに沿って解説します．

方針は単純です．まず 1 本の完成版コード (Appendix)で，節点座標対応表，要素接点対応表，境界辺接点対応表，メッシュ図，そして MinAngle / MaxAngle 型の比較表までまとめて出せるようにします．そのうえで，Section 9.1 の good family，Section 9.2 の指標，Section 9.3 の bad family を順に読み解きます．

1．このページで何を解説するのか

ここでは，記事本文は「Section 9 の意味を説明すること」に集中し，詳細な実装は完成版コード (Appendix)にまとめました．これにより，読者はまず「何を再現したいのか」を理解し，そのあとで「どの関数が何を担当しているのか」を追えるようになります．

特に Section 9 では，9.1 に (I) Standard，(II) Shishkin，(III) anisotropic mesh from [25]，(IV) graded mesh が与えられ，9.2 では MinAngle と MaxAngle 型の量が導入され，9.3 では bad family として (V) Shishkin mesh，(VI) Graded mesh が挙げられています．この記事では，この流れそのものを崩さないことを優先します．注意: (V)は bad と分類していますが，古典的意味で収束はします．理論に係数が悪くなる可能性があるという意味で bad に分類しています．

2．完成版コード (Appendix)の役割

完成版コードの役割は次の6点です．

Section 9.1 の (I)–(IV) の節点列を生成すること．
Section 9.3 の (V)，(VI) を，x 方向は staggered row，y 方向は Shishkin/graded で構成すること．
節点座標対応表を出すこと．
要素接点対応表と境界辺接点対応表を出すこと．
メッシュ図を描き，必要なら節点番号・要素番号・境界辺番号を重ねること．
MinAngle と MaxAngle 型の量を計算し，6つの family を一つの表で比較すること．

つまり，このファイル 1 本を見れば，Section 9 の mesh family を「図」と「表」と「指標」の三つの視点からまとめて確認できます．

3．Section 9.1：good family (I)–(IV)

Section 9.1 の (I)–(IV) は，まず節点列 \((x_1^i, x_2^i)\) を与え，その tensor-product grid を三角形に分割することで得られます．完成版コードでは，この部分を section9_lines(kind, N, ...) が担当しています．

(I) Standard mesh：x 方向・y 方向ともに一様分割です．
(II) Shishkin mesh：x 方向は一様，y 方向のみ Shishkin 型に分割します．
(III) Anisotropic mesh from [25]：x 方向・y 方向とも cosine clustering を使います．
(IV) Graded mesh：x 方向は一様，y 方向は \((i/N)^2\) とします．

この4つは good family 側の導入例として使いやすく，節点座標対応表と要素接点対応表を出すだけでも，異方性の入り方の違いがかなりよく見えます．

4．Section 9.2：`MinAngle` と `MaxAngle` をどう読むか

Section 9.2 では，各三角形の 3 辺長を \(|L_1| \le |L_2| \le |L_3|\) の順に並べたうえで，

\[ \mathrm{MinAngle}:=\max_{T\in\mathcal T_h}\frac{|L_3|^2}{|T|_2}, \qquad \mathrm{MaxAngle}:=\max_{T\in\mathcal T_h}\frac{|L_1||L_2|}{|T|_2} \]

を計算しています．完成版コードでは，tri_edge_lengths，tri_area，section9_indicators がこの部分を担当します．

ここで大事なのは，辺長は厳密に L1 < L2 < L3 ではなく，L1 <= L2 <= L3 と並べることです．等辺や二等辺の要素では辺長が一致することがあるからです．完成版コードでも，np.sort(...) によりこの順序を保証しています．

5．Section 9.3：bad family (V), (VI)

Section 9.3 では (V) Shishkin mesh と (VI) Graded mesh が bad family として現れます．ここで注意すべきなのは，(V)，(VI) を「別の y 方向節点列」として考えるだけでは不十分で，x 方向の切り方そのものが 9.1 と異なることです．

このコードでは，odd row は通常格子，even row は半ステップずれた格子になっており，さらに odd strip/even strip ごとに三角形の貼り方を変えています．完成版コードでは，この考え方を mesh1_style_rows と build_mesh1_style_mesh に移しました．

したがって，(V)，(VI) の構成は次のように読むのが自然です．

(V)：y 方向は Shishkin，x 方向は staggered row．
(VI)：y 方向は graded，x 方向は staggered row．

これにより，図だけでは分かりにくかった「どこが 9.1 の family と違うのか」が，節点列と要素接続の両方から読めるようになります．

6．A型/B型と good/bad の見方

本ページでは，メッシュ退化を二つに分けて読みます．

A型：形状パラメータが一様有界な相似縮小型．
B型：精密化と同時に形状パラメータが発散する graded 退化型．

A型では古典理論によって収束次数自体は既知であり，本質は誤差係数の退化依存性にあります．一方，B型では shape-regularity が壊れていても new geometric condition 側では still good になり得るかどうかが中心問題になります．Section 9 を family ごとに表と図で並べる意義は，まさにここにあります．

Section 9.3 では Mesh (V) を Bad にしています．(V) は　shape-regularity　を満たしますが，誤差係数が退化依存します．これにより，Bad に分類しています．しかしながら，収束次数自体は古典理論によって保証されます．

完成版コードは，この good/bad の違いを indicator_table_all という比較表にまとめます．したがって，読者は6つの family を一つの表で見たあとに，その背後にある mesh data を節点表・要素表・境界辺表へ戻って確認できます．

7．完成版コードの読み方

完成版コードは，次の順番で読むと分かりやすいです．

section9_lines：9.1 の (I)–(IV) を作る．
build_tri_mesh_from_tensor_grid：tensor-product grid を三角形化する．
tables_from_mesh と show_tables：三つの表を出す．
section9_indicators：Section 9.2 の指標を計算する．
section93_lines，mesh1_style_rows，build_mesh1_style_mesh：9.3 の (V)，(VI) を作る．
最後の demo 部分：family ごとの図や比較表を出す．

この構成にしておくと，記事本文では個々の関数を必要なときだけ参照すればよく，HTML の中に長大なコードを何度も埋め込まずに済みます．

8．実行のしかた

この整理版では，まず完成版コードを Colab または Jupyter へ貼り付けて実行し，その出力を見ながら本文を読むことを想定しています．最初の example では N=4 の小さい例で節点表・要素表・境界辺表を確認し，その後で6つの family の比較表へ進みます．

記事を読む順番としては，

まず節点座標対応表・要素接点対応表・境界辺接点対応表を見る．
次にメッシュ図で節点番号・要素番号・境界辺番号の対応を確認する．
最後に MinAngle / MaxAngle 型の量で family 全体を比較する．

と進むのが一番分かりやすいです．

9．まとめ

完成版コードを1本に固定し，本文は Section 9 の意味を説明することに集中させました．その結果，

Section 9.1 の good family，
Section 9.2 の指標，
Section 9.3 の bad family，
A型 / B型と good / bad の見方，
staggered rowの役割，

がそれぞれ整理しやすくなりました．今後は，この整理版を土台にして，FEniCSx や Gmsh 側の mesh data とどう対応づけるかを次段階として書いていく方針です．

Appendix 1．完成版コード全文（HTML 版）

以下に，コードの全文を掲載します．この版では，(I) から (VI) までの 6 つの family を 2×3 の図としてまとめて表示する plot_all_six_families も含めています．DELTA をコード冒頭で定義し，Shishkin 型 family では，delta を省略した呼び出しでも自動的に DELTA が使われるようにしてあります．

コード全文を開く/閉じる

# Colab / Jupyter 用の全コード（I)-(VI) 比較つき）
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display

pd.set_option("display.max_rows", None)

# -------------------------
# Global configuration
# -------------------------
# Shishkin-type families (II) and (V) depend on delta.
# Change this value freely when you want to compare different transition layers.
DELTA = 1/128


def x_standard(N):
    return np.linspace(0.0, 1.0, N + 1)


def x_shishkin_y(N, delta=None):
    if delta is None:
        delta = DELTA
    tau = 2 * delta * abs(np.log(N))
    tau = min(max(tau, 0.0), 1.0)
    y = np.zeros(N + 1)
    half = N // 2
    for i in range(N + 1):
        if i <= half:
            y[i] = tau * i / half if half > 0 else 0.0
        else:
            y[i] = tau + (1 - tau) * (i - half) / (N - half)
    return y


def x_cosine_cluster(N):
    i = np.arange(N + 1, dtype=float)
    return 0.5 * (1.0 - np.cos(i * np.pi / N))


def x_graded_square(N):
    i = np.arange(N + 1, dtype=float)
    return (i / N) ** 2


def section9_lines(kind, N, delta=None):
    if delta is None:
        delta = DELTA
    if kind == "I":
        x1 = x_standard(N)
        x2 = x_standard(N)
    elif kind == "II":
        x1 = x_standard(N)
        x2 = x_shishkin_y(N, delta=delta)
    elif kind == "III":
        x1 = x_cosine_cluster(N)
        x2 = x_cosine_cluster(N)
    elif kind == "IV":
        x1 = x_standard(N)
        x2 = x_graded_square(N)
    else:
        raise ValueError("kind must be one of 'I', 'II', 'III', 'IV'.")
    return x1, x2


def build_tri_mesh_from_tensor_grid(x1, x2, diagonal="right", flip_corner_cells=False):
    x1 = np.asarray(x1, dtype=float)
    x2 = np.asarray(x2, dtype=float)
    nx = len(x1) - 1
    ny = len(x2) - 1

    nodes = []
    for j in range(ny + 1):
        for i in range(nx + 1):
            nodes.append([x1[i], x2[j]])
    nodes = np.array(nodes, dtype=float)

    def nid(i, j):
        return j * (nx + 1) + i

    elements = []
    for j in range(ny):
        for i in range(nx):
            n00 = nid(i, j)
            n10 = nid(i + 1, j)
            n01 = nid(i, j + 1)
            n11 = nid(i + 1, j + 1)
            use_diagonal = diagonal
            if flip_corner_cells:
                is_top_left = (i == 0 and j == ny - 1)
                is_bottom_right = (i == nx - 1 and j == 0)
                if is_top_left or is_bottom_right:
                    use_diagonal = "left" if diagonal == "right" else "right"
            if use_diagonal == "right":
                elements.append([n00, n10, n11])
                elements.append([n00, n11, n01])
            elif use_diagonal == "left":
                elements.append([n00, n10, n01])
                elements.append([n10, n11, n01])
            else:
                raise ValueError("diagonal must be 'right' or 'left'.")
    return nodes, np.array(elements, dtype=int)


def boundary_edges_from_elements(elements):
    edge_count = {}
    for tri in elements:
        a, b, c = map(int, tri)
        edges = [tuple(sorted((a, b))), tuple(sorted((b, c))), tuple(sorted((c, a)))]
        for e in edges:
            edge_count[e] = edge_count.get(e, 0) + 1
    boundary_edges = np.array([e for e, c in edge_count.items() if c == 1], dtype=int)
    idx = np.lexsort((boundary_edges[:, 1], boundary_edges[:, 0]))
    return boundary_edges[idx]


def tables_from_mesh(nodes, elements, boundary_edges):
    node_table = pd.DataFrame({
        "node": np.arange(len(nodes), dtype=int),
        "x": nodes[:, 0],
        "y": nodes[:, 1],
    })
    element_table = pd.DataFrame({
        "element": np.arange(len(elements), dtype=int),
        "node0": elements[:, 0],
        "node1": elements[:, 1],
        "node2": elements[:, 2],
    })
    boundary_table = pd.DataFrame({
        "boundary_edge": np.arange(len(boundary_edges), dtype=int),
        "node0": boundary_edges[:, 0],
        "node1": boundary_edges[:, 1],
    })
    return node_table, element_table, boundary_table


def show_tables(node_table, element_table, boundary_table):
    print("=== node table ===")
    display(node_table.style.format({"x": "{:16.12e}", "y": "{:16.12e}"}))
    print("=== element table ===")
    display(element_table)
    print("=== boundary-edge table ===")
    display(boundary_table)


def plot_mesh_with_labels(nodes, elements, boundary_edges=None,
                          show_node_ids=True, show_elem_ids=False,
                          show_boundary_ids=False,
                          title="mesh"):
    plt.figure(figsize=(6, 6))
    for tri in elements:
        tri = list(map(int, tri))
        poly = np.vstack([nodes[tri], nodes[tri[0]]])
        plt.plot(poly[:, 0], poly[:, 1], color="0.45", lw=0.9)
    if boundary_edges is not None:
        for k, (a, b) in enumerate(boundary_edges):
            seg = nodes[[a, b]]
            plt.plot(seg[:, 0], seg[:, 1], color="black", lw=1.5)
            if show_boundary_ids:
                m = seg.mean(axis=0)
                plt.text(m[0], m[1], str(k), fontsize=8, color="tab:purple")
    if show_node_ids:
        for k, (x, y) in enumerate(nodes):
            plt.text(x, y, str(k), fontsize=8, color="tab:blue")
    if show_elem_ids:
        for k, tri in enumerate(elements):
            pts = nodes[np.array(tri, dtype=int)]
            c = pts.mean(axis=0)
            plt.text(c[0], c[1], str(k), fontsize=8, color="tab:red")
    plt.gca().set_aspect("equal")
    plt.xlabel("x")
    plt.ylabel("y")
    plt.title(title)
    plt.show()


def tri_edge_lengths(pts):
    L = [
        np.linalg.norm(pts[1] - pts[0]),
        np.linalg.norm(pts[2] - pts[1]),
        np.linalg.norm(pts[0] - pts[2]),
    ]
    return np.sort(np.array(L, dtype=float))  # L1 <= L2 <= L3


def tri_area(pts):
    a = pts[1] - pts[0]
    b = pts[2] - pts[0]
    return 0.5 * abs(a[0] * b[1] - a[1] * b[0])


def section9_indicators(nodes, elements):
    min_like = 0.0
    max_like = 0.0
    for tri in elements:
        pts = nodes[tri]
        L1, L2, L3 = tri_edge_lengths(pts)
        area = tri_area(pts)
        min_like = max(min_like, (L3 ** 2) / area)
        max_like = max(max_like, (L1 * L2) / area)
    return float(min_like), float(max_like)


# --- bad family based on mesh1.h ---
def section93_lines(kind, N, delta=None):
    if delta is None:
        delta = DELTA
    """
    Section 9.3 の (V), (VI) 用の y 方向節点列．
    (V)  : Shishkin mesh なので 9.1 の (II) と同じ縦方向節点列
    (VI) : Graded mesh   なので 9.1 の (IV) と同じ縦方向節点列
    """
    if kind == "V":
        return x_shishkin_y(N, delta=delta)
    elif kind == "VI":
        return x_graded_square(N)
    else:
        raise ValueError("kind must be 'V' or 'VI'.")


def mesh1_style_rows(N):
    """
    mesh1.h 型の x 方向節点列を行ごとに作る．
    偶数番目の行(0-based)は通常格子，奇数番目の行は半ステップずれた格子にする．
    """
    x_base = np.linspace(0.0, 1.0, N + 1)
    mids = 0.5 * (x_base[:-1] + x_base[1:])
    rows = []
    for j in range(N + 1):
        if j % 2 == 0:
            rows.append(x_base.copy())
        else:
            rows.append(np.concatenate(([x_base[0]], mids, [x_base[-1]])))
    return rows


def build_mesh1_style_mesh(N, y_lines):
    """
    mesh1.h の odd/even row と要素接続を Python に移した 2D 三角形メッシュ生成．
    x 方向は staggered row，y 方向は外から与える．
    """
    y_lines = np.asarray(y_lines, dtype=float)
    if len(y_lines) != N + 1:
        raise ValueError("y_lines must have length N+1")

    rows_x = mesh1_style_rows(N)

    nodes = []
    row_starts = []
    count = 0
    for j in range(N + 1):
        row_starts.append(count)
        for x in rows_x[j]:
            nodes.append([x, y_lines[j]])
            count += 1
    nodes = np.array(nodes, dtype=float)

    def nid(row, k):
        return row_starts[row] + k

    elements = []
    for j in range(N):
        if j % 2 == 0:
            # mesh1.h の odd strip に対応（0-based では j 偶数）
            for i in range(N):
                p1 = nid(j, i)
                p2 = nid(j, i + 1)
                p3 = nid(j + 1, i + 1)
                elements.append([p1, p2, p3])
                p1 = nid(j, i)
                p2 = nid(j + 1, i + 1)
                p3 = p2 - 1
                elements.append([p1, p2, p3])
            p1 = nid(j, N)
            p2 = nid(j + 1, N + 1)
            p3 = p2 - 1
            elements.append([p1, p2, p3])
        else:
            # mesh1.h の even strip に対応（0-based では j 奇数）
            p1 = nid(j, 0)
            p2 = p1 + 1
            p3 = nid(j + 1, 0)
            elements.append([p1, p2, p3])
            for i in range(N):
                p1 = nid(j, i + 1)
                p2 = p1 + 1
                p3 = nid(j + 1, i + 1)
                elements.append([p1, p2, p3])
                p1 = nid(j, i + 1)
                p2 = nid(j + 1, i + 1)
                p3 = p2 - 1
                elements.append([p1, p2, p3])
    return nodes, np.array(elements, dtype=int)


def section93_reconstruction(kind, N, delta=None):
    if delta is None:
        delta = DELTA
    """
    Section 9.3 の (V), (VI) を作る．
    x 方向は mesh1.h 型，y 方向は (V) で Shishkin, (VI) で graded を使う．
    """
    y_lines = section93_lines(kind, N, delta=delta)
    return build_mesh1_style_mesh(N, y_lines)


def build_family_mesh(kind, N, delta=None):
    if delta is None:
        delta = DELTA
    """
    kind = "I", ..., "VI" に対して nodes, elements, boundary_edges を返す．
    """
    if kind in ["I", "II", "III", "IV"]:
        x1, x2 = section9_lines(kind, N, delta=delta)
        nodes, elements = build_tri_mesh_from_tensor_grid(
            x1, x2, diagonal="right", flip_corner_cells=True
        )
    elif kind in ["V", "VI"]:
        nodes, elements = section93_reconstruction(kind, N, delta=delta)
    else:
        raise ValueError("kind must be one of 'I', 'II', 'III', 'IV', 'V', 'VI'.")

    boundary_edges = boundary_edges_from_elements(elements)
    return nodes, elements, boundary_edges


def plot_family_on_ax(ax, kind, N=16, delta=None,
                      
                      show_node_ids=False,
                      show_elem_ids=False,
                      show_boundary_ids=False):
    """
    1 つの family を与えられた Axes に描く．
    """
    if delta is None:
        delta = DELTA
    nodes, elements, boundary_edges = build_family_mesh(kind, N, delta=delta)

    for tri in elements:
        tri = list(map(int, tri))
        poly = np.vstack([nodes[tri], nodes[tri[0]]])
        ax.plot(poly[:, 0], poly[:, 1], color="0.45", lw=0.9)

    for k, (a, b) in enumerate(boundary_edges):
        seg = nodes[[a, b]]
        ax.plot(seg[:, 0], seg[:, 1], color="black", lw=1.4)
        if show_boundary_ids:
            m = seg.mean(axis=0)
            ax.text(m[0], m[1], str(k), fontsize=7, color="tab:purple")

    if show_node_ids:
        for k, (x, y) in enumerate(nodes):
            ax.text(x, y, str(k), fontsize=7, color="tab:blue")

    if show_elem_ids:
        for k, tri in enumerate(elements):
            pts = nodes[np.array(tri, dtype=int)]
            c = pts.mean(axis=0)
            ax.text(c[0], c[1], str(k), fontsize=7, color="tab:red")

    name_map = {
        "I": "I : Standard",
        "II": "II : Shishkin",
        "III": "III : Anisotropic",
        "IV": "IV : Graded",
        "V": "V : Shishkin (bad family)",
        "VI": "VI : Graded (bad family)",
    }

    ax.set_aspect("equal")
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_title(f"{name_map[kind]}  (N={N})")


def plot_all_six_families(N=16, delta=None,
                          
                          show_node_ids=False,
                          show_elem_ids=False,
                          show_boundary_ids=False):
    """
    I から VI までの 6 つの family を 2x3 に並べて描く．
    """
    if delta is None:
        delta = DELTA
    fig, axes = plt.subplots(2, 3, figsize=(15, 9))
    for ax, kind in zip(axes.ravel(), ["I", "II", "III", "IV", "V", "VI"]):
        plot_family_on_ax(
            ax, kind, N=N, delta=delta,
            show_node_ids=show_node_ids,
            show_elem_ids=show_elem_ids,
            show_boundary_ids=show_boundary_ids,
        )
    plt.tight_layout()
    plt.show()


# --- examples ---
# You can freely change DELTA above.
# Families (II) and (V) use this value automatically when delta is omitted.
kind = "IV"
N = 4
x1, x2 = section9_lines(kind, N)
nodes, elements = build_tri_mesh_from_tensor_grid(x1, x2, diagonal="right", flip_corner_cells=True)
boundary_edges = boundary_edges_from_elements(elements)
node_table, element_table, boundary_table = tables_from_mesh(nodes, elements, boundary_edges)
show_tables(node_table, element_table, boundary_table)
plot_mesh_with_labels(nodes, elements, boundary_edges, show_node_ids=True, show_elem_ids=True,
                      show_boundary_ids=False, title=f"family {kind}: nodes + elements")
plot_mesh_with_labels(nodes, elements, boundary_edges, show_node_ids=False, show_elem_ids=False,
                      show_boundary_ids=True, title=f"family {kind}: boundary-edge ids")

# full comparison (I)-(VI)
rows_all = []
for kind in ["I", "II", "III", "IV"]:
    for N in [16, 32, 64, 128]:
        x1, x2 = section9_lines(kind, N)
        nodes_g, elements_g = build_tri_mesh_from_tensor_grid(x1, x2, diagonal="right", flip_corner_cells=True)
        min_like, max_like = section9_indicators(nodes_g, elements_g)
        rows_all.append({"Mesh": kind, "Type": "good family", "N": N,
                         "MinAngle": min_like, "MaxAngle": max_like})
for kind in ["V", "VI"]:
    for N in [16, 32, 64, 128]:
        nodes_b, elements_b = section93_reconstruction(kind, N)
        min_like, max_like = section9_indicators(nodes_b, elements_b)
        rows_all.append({"Mesh": kind, "Type": "bad-family reconstruction", "N": N,
                         "MinAngle": min_like, "MaxAngle": max_like})

indicator_table_all = pd.DataFrame(rows_all)
display(indicator_table_all.style.format({"MinAngle": "{:16.12e}", "MaxAngle": "{:16.12e}"}))



# -------------------------
# Demo: plot bad family reconstruction (V), (VI)
# -------------------------
fig, axes = plt.subplots(1, 2, figsize=(10, 4.5))
for ax, kind in zip(axes, ["V", "VI"]):
    nodes_b, elements_b = section93_reconstruction(kind, 16)
    boundary_b = boundary_edges_from_elements(elements_b)
    for tri in elements_b:
        tri = list(map(int, tri))
        poly = np.vstack([nodes_b[tri], nodes_b[tri[0]]])
        ax.plot(poly[:, 0], poly[:, 1], color="0.4", lw=0.8)
    for e in boundary_b:
        pts = nodes_b[list(e)]
        ax.plot(pts[:, 0], pts[:, 1], color="tab:red", lw=1.2)
    ax.set_aspect("equal")
    ax.set_title(f"family {kind}")
    ax.set_xlabel("x")
    ax.set_ylabel("y")
plt.tight_layout()
plt.show()


# -------------------------
# Demo: plot all six families
# -------------------------
plot_all_six_families(N=16)

Appendix 2．完成版コード全文（HTML 版）

この版では，N，DELTA，EPSILON を global configuration としてまとめ，それ以外の箇所では数値を直書きしないように整理しています．CASE = make_case(...) を1回与えれば，1 family の結果，6 family の比較表，6 family の図を同じ設定で再利用できます．

以下に，現時点の完成版コード全文を HTML 表示用に埋め込みます．この版では，DELTA をコード冒頭で定義し，Shishkin 型 family では，delta を省略した呼び出しでも自動的に DELTA が使われるようにしてあります．