LoginSignup
7
7

More than 5 years have passed since last update.

@shohirose(広瀬 翔)

メッシュとイテレータとサーキュレータ

Last updated at Posted at 2019-03-16

はじめに

計算力学・数値流体力学の分野において領域の差分化はほぼ必須であり、したがって、差分化から得られるメッシュのデータ構造およびそれに対する操作をどのように定義し実装するかというのは非常に重要な問題です。

最も単純なアイデアは、頂点・辺・面・セルのそれぞれに固有の番号を振り、この番号を基に各メッシュ要素およびそれに付随するデータにアクセスするというものです。各要素が（std::vectorなどの）配列を使って保存されている場合は「固有の番号＝配列の要素番号」であるため、特にFortranやCでプログラムする場合に有効な手法です。

C++で実装するとしたら、下記のように各要素をstd::vectorで保持し、アクセス関数を提供する単純なクラスになります。

volume_mesh
#include <vector>

// ... point/edge/face/cellのヘッダーファイルを読み込み ... //

class volume_mesh {
 public:
  // ... ctor/dtorを定義 ... //

  // 各配列へのアクセス関数を提供
  const auto& points() const noexcept { return points_; }
  const auto& edges() const noexcept { return edges_; }
  const auto& faces() const noexcept { return faces_; }
  const auto& cells() const noexcept { return cells_; }

  // ... その他アクセス関数を定義 ... //

 private:
  std::vector<point> points_;
  std::vector<edge> edges_;
  std::vector<face> faces_;
  std::vector<cell> cells_;

  // ... その他のデータ ... //
};

ただ、メッシュ要素に固有の番号としてint/unsigned intなどを使うため、

  • 異なるメッシュ要素の番号を間違って使ってしまう可能性がある。
  • メッシュに対する操作が内部のデータ構造に依存してしまうため、アルゴリズム×メッシュの数だけ実装が必要となる。

などの問題が生じます。

この記事では、C++のメッシュライブラリがこのような問題をどのように解決しているか紹介します。

Strongly-typed index

異なるメッシュ要素の番号を区別するため、Boost Graph Library(BGL)、OpenMesh、CGAL - Surface Meshではstrongly-typed indexを使用しています。これは、インデックスの型にどのメッシュ要素に属するかを型情報として含めるものです。

strongly-typed-index
template <typename Derived>
class index {
 public:
  using size_type = std::size_t;
  explicit index(size_type idx) : idx_{idx} {}
  // size_typeにstatic_castできるよう演算子を定義
  explicit operator size_type() const noexcept { return idx_; }
  // ... その他コンストラクタ、演算子、アクセス関数等を定義 ... //
 private:
  size_type idx_{std::numeric_limits<size_type>::max()};
};

class vertex_index : public index<vertex_index> {
  using base_type = index<vertex_index>;
 public:
  explicit vertex_index(std::size_t idx) : base_type(idx) {}
  // ... その他コンストラクタを定義 ... //
};

class edge_index : public index<edge_index> { /* ... */ };
class face_index : public index<face_index> { /* ... */ };
class cell_index : public index<cell_index> { /* ... */ };

より詳しくは過去の記事

もどうぞ。

Strongly-typed indexはメッシュ要素に固有の番号であり、より一般的にはメッシュ要素へのハンドルとみなすことができます。したがって、配列の要素番号とは全く関係ない値を使っても、対応するメッシュ要素へのアクセス方法さえ提供しているのならば問題ありません。

PropertyMap

メッシュ要素固有の番号を元にして、そのメッシュ要素に関連するプロパティにアクセスする方法を提供するのがBGLのPropertyMapです。PropertyMapについては、SlideShare: Boost.勉強会 #3 関西 プレゼン資料がよくまとめられていてわかりやすいです。

メッシュ要素固有の番号が配列の要素番号と同一の場合は、[]演算子の引数としてstrongly-typed indexを取る配列クラスを定義するのが最も簡単な方法です。例えばstd::vectorを継承してPropertyMap`を作ってみましょう。

property_array
template <typename T, typename Key, typename Allocator = std::allocator<T>>
class property_array : public std::vector<T, Allocator> {
  using base_type = std::vector<T, Allocator>;
 public:
  using typename base_type::size_type;
  /* ... その他value_type等を定義 ... */

  T& operator[](const Key& key) {
    return this->base_type::operator[](static_cast<size_type>(key));
  }
  const T& operator[](const Key& key) const {
    return this->base_type::operator[](static_cast<size_type>(key));
  }
};

template <typename T>
using vertex_property = property_array<T, vertex_index>;
template <typename T>
using edge_property = property_array<T, edge_index>;
template <typename T>
using face_property = property_array<T, face_index>;
template <typename T>
using cell_property = property_array<T, cell_index>;

このproperty_arrayを使うと、基本的に配列要素へKeyを使ってアクセスするように強制されます。例えばvertex_propertyならばvertex_indexを使って要素を指定することになり、配列がメッシュのどの要素に紐付けられているのかというのが型情報に含まれることになります。

property_array
face_property<point> face_centers;
/* ... 初期化 ... */
for (std::size_t i = 0; i < mesh.nfaces(); ++i) {
  const auto f = face_index(i);
  const auto& fc = face_centers[f]; // face_indexを通して要素へアクセス
  // ...
}

Index iterator

前節のようにPropertyMapを定義すると、要素間でインデックスを取り違える心配がありません。一方、前節の例のように全てのループで

  1. std::size_tを使ってforループを回す
  2. ループ内で各要素のインデックスに変換
  3. コンテナにアクセス

するという風に手間がかかります。

CGALのSurface Meshでは、Index iteratorを使ってこの手間を省いています。Index iteratorは、strongly-typed indexを内部に保持し、イテレータと同様の操作を提供します。またdereferenceするとstrongly-typed index自体を返します。

index_iterator
template <typename Index>
class index_iterator {
 public:
  index_iterator(Index idx) : idx_(idx) {}
  // ... その他ctor/dtorを定義 ... //

  index_iterator& operator++() noexcept {
    ++idx_;
    return *this;
  }
  // ... その他演算子を定義 ... //

  const Index& operator*() const noexcept { return idx_; }
  Index& operator*() noexcept { return idx_; }

 private:
  Index idx_;
};

vertex_iterator = index_iterator<vertex_index>;
edge_iterator = index_iterator<edge_index>;
face_iterator = index_iterator<face_index>;
cell_iterator = index_iterator<cell_index>;

またrange-based for用にiterator_rangeを定義します。

iterator_range
template <typename Iterator>
class iterator_range {
 public:
  iterator_range(const Iterator& b, const Iterator& e) : first_(b), second_(e) {}
  // ... ctor/dtorを定義 ... //

  Iterator begin() const noexcept { return first_; }
  Iterator end() const noexcept { return second_; }
  // ... その他関数 ... //
 private:
  Iterator first_;
  Iterator second_;
};

using vertex_range = iterator_range<vertex_iterator>;
using edge_range = iterator_range<edge_iterator>;
using face_range = iterator_range<face_iterator>;
using cell_range = iterator_range<cell_iterator>;

template <typename Iterator>
iterator_range<Iterator> make_iterator_range(const Iterator& b, const Iterator& e) {
  return iterator_range<Iterator>(b, e);
}

メッシュクラスからiterator_rangeオブジェクトを返すようにするとrange-based forを使って書けるようになります。例えば

volume_mesh
class volume_mesh {
 public:
  // ...
  auto nfaces() const noexcept { return face_.size(); }
  face_iterator faces_begin() const noexcept { return face_iterator(0); }
  face_iterator faces_end() const noexcept { return face_iterator(this->nfaces()); }
  // [0,nfaces)の範囲を返す
  face_range faces() const noexcept {
    return make_iterator_range(this->faces_begin(), this->faces_end());
  }
  // ...
 private:
  // ...
  std::vector<face> faces_;
  // ...
};

とvolume_meshを書き換えると、前節の例は次のようにシンプルになります。

face_property<point> face_centers;
/* ... 初期化 ... */
for (auto&& f : mesh.faces()) {
  const auto& fc = face_centers[f];
  // ...
}

とりあえずここまで。後でCirculatorについて追記します。

Circulator

メッシュに関連する操作には、例えば

  • ある頂点に接続する辺の全てに対する操作
  • ある面を構成する全ての頂点に対する操作

など、基準となるメッシュ要素の周りを循環して行う操作が多いです。
サーキュレータ（circulator）とは、この循環操作を一般化したイテレータの一種です。より具体的には

  • 基準となるメッシュ要素＝anchorをもとに作られる。
  • サーキュレータの表す範囲は[c,c)、すなわち範囲の先頭と末尾が同じ。

という点が通常のイテレータと異なります。

最も簡単なサーキュレータの定義は、以下のように先頭および末尾のイテレータを保持し、

  • 範囲の先頭でdecrementされたら末尾を指し示す。
  • 範囲の末尾でincrementされたら先頭を指し示す。

というものです。

circulator
template <typename Iterator>
class circulator {
 public:
  circulator(const Iterator& b, const Iterator& e)
    : begin_(b), end_(e), current_(b) {}
  circulator(const Iterator& b, const Iterator& e, const Iterator& c)
    : begin_(b), end_(e), current_(c) {}
  // ... その他ctor/dtor、アクセス関数などを定義 ... //

  circulator& operator++() {
    Expects(current_ != end_);
    ++current_;
    if (current_ == end_) current_ = begin_;
    return *this;
  }
  circulator& operator--() {
    if (current_ == begin_) current_ = end_;
    --current_;
    return *this;
  }
  // ... その他演算子を定義 ... //

 private:
  Iterator begin_;
  Iterator end_;
  Iterator current_;
};

サーキュレータの範囲は[c,c)であり、これはrange-based forでは空の範囲と判断されてしまうため、range-based forが使えません。そのため、通常はdo-while文を使ってループします。

face_property<point> face_centers;
// 面の中心座標を計算
for (auto&& f : mesh.faces()) {
  point c(0,0,0);
  double n = 0.0;
  auto vc = vc_end = mesh.vertices(f);  // faceを基準にvertexを循環するサーキュレータ
  do {
    c += mesh.point(*vc);
    ++n;
  } while (++vc != vc_end);
  face_centers[f] = c / n;
}

とりあえずここまで。mesh.verticesなどの実装はまた後ほど。

7
7
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
7
7