OpenFOAM ユーザーのための C++ 入門

2011年6月5日
春日 悠

はじめに

この文書では、C++ に慣れていない OpenFOAM ユーザー向けに、C++ の初学者が特につまづきやすいと思われる概念について解説を試みます。読者は C++ の基本は学習済みと想定します。

プログラミングの技術は、抽象的なものにしろ具体的なものにしろ、できるだけコーディングをサボるためにあります。以下で説明する機能も、そういうものとして見ると理解しやすいと思います。

継承

オブジェクト指向プログラミングでは、あるクラスの機能を引き継いで新しいクラスを作ることができます。この機能のことを「継承」といいます。継承を使うと、新しいクラスと継承元のクラスとで異なっている部分だけをプログラミングすればよいので、これを「差分プログラミング」と言ったりします。

さて、ここで例として、あるファンタジー系のゲームのキャラクタのクラスを実装することを考えましょう。「人間」という種族がいて、その中には「戦士」あるいは「魔法使い」を職業とする人がいるものとします。それぞれの職業にはそれ特有の「技」や「魔法」などの特殊技能があるものとします。

#include <iostream>
#include <string>


using namespace std;


class Human
{
public:
	Human(const char *name) : _name(name) {};
	void sayName() { cout << _name << endl; };
	void special() { cout << "special" << endl; };
private:
	string _name;
};


class Fighter : public Human
{
public:
	Fighter(const char *name) : Human(name) {};
	void special() { skill(); };
private:
	void skill() { cout << "skill" << endl; };
};


class Magician : public Human
{
public:
	Magician(const char *name) : Human(name) {};
	void special() { magic(); };
private:
	void magic() { cout << "magic" << endl; };
};


int main()
{
	Human h("human1");
	Fighter f("fighter1");
	Magician m("magician1");

	h.sayName();
	h.special();

	f.sayName();
	f.special();

	m.sayName();
	m.special();

	return 0;
}

「人間」(Human) としての特性はいろいろありますが、ここでは簡単のために「名前」(name) だけを持つものとしています。「名前」を出力するメンバ関数 sayName() を用意しています。また、「人間」の各職業にはそれぞれ特有の特殊技能があるものとしたので、それを実行するメンバ関数 special() を持たせています。

「人間」(Human) を継承して「戦士」(Fighter)、「魔法使い」(Magician) クラスを作っています。それぞれ「技」(skill) と「魔法」(magic) の技能を持っていて、「人間」(Human) から継承したメンバ関数 special() をオーバーライド (上書き) してそれぞれの技能を実行するようにしています。

main() 関数では、「人間」(Human)、「戦士」(Fighter)、「魔法使い」(Magician) それぞれのオブジェクトを作成して、sayName() と special() を実行させています。

出力は次のようになります。

human1
special
fighter1
skill
magician1
magic

それぞれの名前と技能が出力されています。

ところで、あるクラスのオブジェクトを、その親クラスのポインタで指すことができます。

Human *fp = &f;
Human *mp = &m;

fp->sayName();
fp->special();

mp->sayName();
mp->special();

この場合、Fighter は Fighter として、Magician は Magician として振舞ってほしいところですが、そうはなりません。出力は次のようになります。

fighter1
special
magician1
special

親クラスのポインタで子クラスを指したときに、指しているそれぞれの子クラスとしての振る舞いを期待するならば、その場合はオーバーライドされるメンバ関数を「仮想関数」にします。

class Human
{
public:
	...
	virtual void special() { cout << "special" << endl; };
	...

メンバ関数の定義の先頭に "virtual" を付けると、それは仮想関数になります。仮想関数が子クラスによってオーバーライドされた場合、子クラスのインスタンスを指す親クラスのポインタからその関数が呼び出されると、子クラスの関数が実行されます。つまり、先ほどのコードの出力が次のようになります。

fighter1
skill
magician1
magic

このように、操作は同じであるがインスタンスに合わせて挙動が変わる性質を「多態性 (ポリモーフィズム)」といいます。

「人間」というのは総称であって、実在ではありません。「人間」に属するものはすべて特殊技能を持ちますが、どういう特殊技能を持つかはそれぞれの職業で異なります。したがって、「人間」は「特殊技能を持つ」ということだけを知っていればよく、それが実際にどのようなものかは「人間」自体は知る必要がありません。このような状況は、「純粋仮想関数」を使うことで実現できます。

class Human
{
public:
	...
	virtual void special() = 0;
	...

仮想関数の宣言の末尾に "= 0" を付けると純粋仮想関数になります。これは、メンバ関数としては存在するが実装を持たず、実装は子クラスに任せるという宣言です。具体的な実装を持たないため、Human クラスのオブジェクトを作ることはできません。このように純粋仮想関数を持つクラスのことを抽象クラスと呼びます。すべてのメンバ関数が純粋仮想関数であるクラスのことを、クラスのインターフェイスだけを定義しているという意味で「インターフェイスクラス」と呼ぶことがあります。

抽象クラスは次のような使い方ができます。

#include <iostream>
#include <string>
#include <cassert>


using namespace std;


class Human
{
public:
	Human(const char *name) : _name(name) {};
	virtual ~Human() {};
	void sayName() { cout << _name << endl; };
	virtual void special() = 0;
private:
	string _name;
};


class Fighter : public Human
{
public:
	Fighter(const char *name) : Human(name) {};
	void special() { skill(); };
private:
	void skill() { cout << "skill" << endl; };
};


class Magician : public Human
{
public:
	Magician(const char *name) : Human(name) {};
	void special() { magic(); };
private:
	void magic() { cout << "magic" << endl; };
};


enum CharacterType {
	fighter,
	magician,
};


Human* makeCharacter(const char *name, CharacterType type)
{
	switch(type)
	{
	case fighter:
		return new Fighter(name);
		break;
	case magician:
		return new Magician(name);
		break;
	default:
		assert(0);
		break;
	}
}


int main()
{
	Human *h1;
	Human *h2;

	h1 = makeCharacter("fighter1", fighter);
	h2 = makeCharacter("magician1", magician);

	h1->sayName();
	h1->special();

	h2->sayName();
	h2->special();

	delete h1;
	delete h2;

	return 0;
}

名前と職業を指定するとそのインスタンスを作成してくれる makeCharacter() という関数を定義しています。インスタンスは「人間」(Human) のポインタで受けます。中身は「戦士」(Fighter) だったり「魔法使い」(Magician) だったりしますが、その詳細は気にせずにあくまで「人間」として扱うことができます。こうすることで、共通の処理を統一して記述することができます。

OpenFOAM の例

OpenFOAM のライブラリは、だいたい 1 つの抽象クラスとそれを継承して作られた複数の具体的実装クラスのセットの形になっています。たとえば、非圧縮性流体用のレイノルズ平均応力モデルを見てみましょう。これは RASModel が抽象クラスで、その実装は RASModel を継承した kEpsilon (標準 k-εモデル) や kOmegaSST (SST k-ω モデル), LRR (Launder-Reece-Rodi レイノルズ応力モデル) といったクラスで行われています。たとえば RASModel では k(), epsilon() などが純粋仮想関数として定義されており、その実装は各乱流モデルに任されています。kEpsilon ではこれはそれぞれそのまま k とεを返せばよいだけですが、kOmega (標準 k-ω モデル) ではεを k とωから計算して返しています。

ソルバーは各乱流モデルのライブラリではなく、非圧縮性流体用のレイノルズ平均応力モデル一般である RASModel を使います。具体的にどの乱流モデルのクラスが機能するかは、ケースの設定によって決まります。こうしてインターフェイスと実装を分離することで、ソルバー実装者は各乱流モデルの実装に煩わされずに済みます。乱流モデルごとにソルバーを作り替えたり、あるいは乱流モデルごとに処理を場合分けしたりする必要がなく、柔軟に乱流モデルを切り替えることができます。また、新しい乱流モデルを追加することも簡単に行えます。

テンプレート

二つの変数の値を入れ替える関数 swap() を考えます。整数でも倍精度浮動小数点数でも使いたい場合、難しいことを考えなければ、それぞれ専用の関数を定義すればよいでしょう。

void swapInt(int *i, int *j)
{
	int t;

	t = *j;
	*j = *i;
	*i = t;
}


void swapDouble(double *a, double *b)
{
	double t;

	t = *b;
	*b = *a;
	*a = t;
}

swapInt() と swapDouble() は変数の型が違うだけで、処理の記述は同じです。これは労力の無駄に感じます。関数の多重定義 (オーバーロード) を使うと、関数名は統一できます。

void swap(int *i, int *j)
{
	int t;

	t = *j;
	*j = *i;
	*i = t;
}


void swap(double *a, double *b)
{
	double t;

	t = *b;
	*b = *a;
	*a = t;
}

変数の型によって関数を使い分ける必要はなくなりましたが、同じ処理を 2 度書いていることには変わりありません。

プログラミングの技術として、ジェネリック (総称) プログラミングという考え方があります。これは、型を引数にしてしまうというものです。これであれば、型の引数を使って代表的な処理を 1 つ記述しておき、その引数に個々の型を当てはめることでそれぞれの型用の関数ができあがります。

C++ はジェネリックプログラミング用に「テンプレート」という仕組みを用意しています。テンプレートで swap() を定義すると、つぎのようになります。

template<class T> void swap(T *a, T *b)
{
	T t;

	t = *b;
	*b = *a;
	*a = t;
}

この swap() を呼び出す場合、引数に int 型変数のポインタが渡されると T = int、double 型変数のポインタが渡されると T = double となり、swap() はそれぞれの型用の関数として柔軟に機能します。

テンプレートに適用される型は、上記のように類推される場合もあれば、明示的に指定する必要がある場合もあります。たとえば STL (標準テンプレートライブラリ) の map (連想配列) を使う場合、キーの型とデータの型をそれぞれ指定する必要があります。

map<string, int> m;

OpenFOAM の例

OpenFOAM では基本的なクラスはテンプレートで実装され、より具体的な型はそれらに型を当てはめる形で定義されています。たとえば fvPatchField や GeometricField はテンプレートで記述され、fvPatchScalarField や volVectorField はそれらの特殊化として定義されています。

fvPatchField.H

template<class Type>
class fvPatchField
:
    public Field<Type>
{
    ...

fvPatchFieldsFwd.H

typedef fvPatchField<scalar> fvPatchScalarField;
typedef fvPatchField<vector> fvPatchVectorField;

GeometricField.H

template<class Type, template<class> class PatchField, class GeoMesh>
class GeometricField
:
    public DimensionedField<Type, GeoMesh>
{
    ...

volFieldsFwd.H

typedef GeometricField<scalar, fvPatchField, volMesh> volScalarField;
typedef GeometricField<vector, fvPatchField, volMesh> volVectorField;

演算子の多重定義

C++ では演算子 ("+" や "-" など) の挙動をユーザーが定義 (多重定義) することができます。

たとえば、ベクトルを扱うモジュールを考えましょう。

#include <iostream>


using namespace std;


struct Vector {
	double value[3];
};


void Vector_show(Vector *v)
{
	for(int i = 0; i < 3; i++){
		cout << v->value[i] << endl;
	}
}


void Vector_add(Vector *a, Vector *b, Vector *c)
{
	for(int i = 0; i < 3; i++){
		c->value[i] = a->value[i] + b->value[i];
	}
}


int main()
{
	Vector a, b, c;

	a.value[0] = 1.; a.value[1] = 2.; a.value[2] = 3.;
	b.value[0] = 4.; b.value[1] = 5.; b.value[2] = 6.;

	cout << "a =" << endl;
	Vector_show(&a);

	cout << "b =" << endl;
	Vector_show(&b);

	Vector_add(&a, &b, &c);

	cout << "c =" << endl;
	Vector_show(&c);

	return 0;
}

できればもっとベクトルらしく "c = a + b" などと記述できればうれしいところです。演算子の多重定義を使うとそれを実現できます。

#include <iostream>


using namespace std;


class Vector {
public:
	void setValue(double x, double y, double z){
		_value[0] = x;
		_value[1] = y;
		_value[2] = z;
	};

	void show(){
		for(int i = 0; i < 3; i++){
			cout << _value[i] << endl;
		}
	};

	double& operator[](size_t i){
		return _value[i];
	};

	friend const Vector operator+(const Vector &a, const Vector &b);
private:
	double _value[3];
};


const Vector operator+(const Vector &a, const Vector &b)
{
	Vector c;

	for(int i = 0; i < 3; i++){
		c._value[i] = a._value[i] + b._value[i];
	}

	return c;
}


int main()
{
	Vector a, b, c;

	a.setValue(1., 2., 3.);
	b.setValue(4., 5., 6.);

	cout << "a =" << endl;
	a.show();

	cout << "b =" << endl;
	b.show();

	c = a + b;

	cout << "c =" << endl;
	c.show();

	a[0] = 7.; a[1] = 8.; a[2] = 9.;

	cout << "a =" << endl;
	for(int i = 0; i < 3; i++){
		cout << a[i] << endl;
	}

	return 0;
}

ベクトル用のクラス Vector を定義しています。operator+() が演算子 "+" の処理の定義です。二項演算子として定義するために、クラスに属さない関数として定義しています。しかし、Vector の成分をプライベートなメンバ変数としているため、それにアクセスするために Vector の宣言の中で operator+() をフレンド関数として指定しています。

ついでに operator[]() も定義しています。これで Vector のオブジェクトを "a[0]" のようにあたかも配列であるかのように扱えます。operator[]() は Vector の成分を入れる配列の要素の参照を返しているため、"a[0] = 1." のように代入することができます。これは関数に代入しているように見える記述 ("a[0] = 1." は "a.operator[](0) = 1." とも書ける) になっているため、よく考えると奇妙ですが、こういう柔軟さが C++ の強力さであり、また、わかりにくいところでもあります。

ついでにもう一つ、演算子の多重定義でわかりにくい例を示しましょう。

b.getA()();

"()" が続いています。これはいったいどういうことでしょう？ 多重定義できる演算子には "()" も含まれます。たとえば次のようなクラスを定義できます。

class A{
public:
	void operator()(){
		cout << "a()" << endl;
	};
};

このクラスのオブジェクトは、次のように使うことができます。

A a;
a();

このように関数のように振る舞うオブジェクトのことを「関数オブジェクト」と呼びます。先ほどのわかりにくい例は、この A のオブジェクトの参照を返す関数を定義すれば実現できます。

class B{
public:
	A& getA(){
		return _a;
	}
private:
	A _a;
};

つまり、こういうことです。

B b;
b.getA()();

オブジェクトの参照を返す関数の返り値にただちに演算子 "()" を適用したため、"()" が続いたのでした。これもまた C++ の柔軟さを示していますが、わかりにくいので自分で書くときは避けた方が賢いでしょう。

OpenFOAM の例

演算子の多重定義は、これがなければ OpenFOAM は存在し得ないと言えるほど、OpenFOAM にとって重要な機能です。OpenFOAM の最大の特徴は、解くべき方程式を離散化方程式の形ではなく、支配方程式に近い形で書けるところにあります。これを実現しているのが演算子の多重定義です。

solve
(
    fvm::ddt(T) - fvm::laplacian(DT, T)
);

上記のコードは laplacianFoam のものですが、演算子の多重定義と多態性を上手く使って支配方程式から離散化方程式を構成する手順を実装することで、抽象レベルの高い記述によるソルバーの実装を可能にしています。