label について

Post date : 2024-08-29 ( update : 2024-12-29 )

label を確認したい場合

buildlayer() 等で label 付けをした場合に、想定通りについているか plot() で確認できます。

mesh Th=square (10,10);

fespace Vh(Th,P1);

int[int] l0=labels(Th); // array of labels

cout << "# of labels=" << l0.n << ", labels:" << l0 << endl; // print information about labels

for(int il=0; il<l0.n; il++){ // l0.n=# of labels

  int l=l0[il]; // label #

  varf vb1(u,v)= on(l,u=1); // set u=1 at the boundary #l

  Vh u;

  u[] = vb1(0,Vh,tgv=1);// u=1 on l-th border ..

  plot(u,wait=1,cmm=" border l="+l); // draw l-th border

}

注意：

• Vh u=label;  では角がわからない

• 複数の領域からなる2次元メッシュを用いて buildlayer() する場合に、境界に label を設定しても region はそのままになっている(領域ごとに番号がふられる)

  • 解決方法：2次元メッシュの region (領域内のラベル)を 0 にリセット　Th=change(Th,fregion=0);

参考：

• https://community.freefem.org/t/get-labels-of-boundary-nodes-for-p2-elements/471/7

mesh データからラベルを取得する方法

mesh Th  のデータへアクセス

• Th(i) : (全体で) i 番目の節点

• Th[k] : k 番目の三角形  <- 全体で Th.nt 個

  • Th[k][i] : k 番目の三角形における i 番目の節点

    • Th[k][i].x, Th[k][i].x : x座標、y座標

    • Th[k][i].label : ラベル番号

例：

for (int i = 0; i < Th.nt; i++)

  for (int j = 0; j < 3; j++)

    cout << i << " " << j << " - Th[i][j] = " << Th[i][j]

         << ", x = " << Th[i][j].x << ", y= " << Th[i][j].y

         << ", label=" << Th[i][j].label << endl;

参考：

• https://doc.freefem.org/documentation/mesh-generation.html#mesh-connectivity-and-data

trunc() で作られた境界のラベリング

作成済みのメッシュから trunc() で領域削除を行うと，新しくできた境界は同じラベルが付けられてしまいます．

面ごとに異なるラベルを設定したい場合は次の方法で実現することができます．

1. trunc(..., label=XXX); // trunc() 実行時に，未使用のラベルを設定

2. func fl=(label == XXX ? (....): lablel); // 各面ごとにラベルづけする関数を作成

3. Th=change(Th, flabel=fl); // 作成した関数を用いてラベル変更

例．立方体領域から小さな直方体を削除し，各面に新たなラベルをつける

load "msh3"

mesh3 Th = cube(20, 20, 20);

// remove sub-region

real xmin=0.1, xmax=0.3, ymin=0.1, ymax=0.2, zmin=0, zmax=0.3;

int ltr=10; // unused label #

Th = trunc(Th, !(xmin < x && x < xmax &&

                ymin < y && y < ymax &&

                zmin < z && z < zmax),label=ltr);

// relabel

func fl=(label==10 ? (x<=xmin ? ltr+1:(x>=xmax ? ltr+2:(y<=ymin ? ltr+3:(y>=ymax ? ltr+4:ltr+5)))):label);

Th=change(Th,flabel=fl);

plot(Th, wait=1);