MEMORAND

粒子法でもやるか

  • math
  • CFD
  • wip

C++で粒子法CFDをつくる

Scalar処理

後で拡張しやすいように、1変数はScalar型として定義しておく

using Scalar = double;

別に↓でもいいが、ここではusingを使用した

typedef double Scalar;

Vector処理

ベクトルを作る

classでデータ構造を定義して、ベクトル計算の処理もクラス内関数として作成する
#include <array>を使用した
理由はベクトルごと比較演算子が使えるから

class Vector {
    std::array<Scalar,3> v;
public:
    Vector() : v{0,0,0} {}
    Vector(Scalar x, Scalar y, Scalar z) : v{x,y,z} {}

Vector(Scalar x...はコンストラクタ定義で、このクラスのインスタンスを生成するときに初期値を設定できるようになる。設定しないと(0, 0, 0)になる
コンストラクタ定義なしでScalar x=0, Scalar y=0, Scalar z=0;と書くと、
例えばVector v(12, 6, 9);みたいな変数定義ができない

成分呼び出し

x成分だけ取り出したいこともあるので、以下のようにした

Scalar& x() { return v[0]; }
Scalar& y() { return v[1]; }
Scalar& z() { return v[2]; }
//const用
const Scalar& x() const { return v[0]; }
const Scalar& y() const { return v[1]; }
const Scalar& z() const { return v[2]; }

Vector vec;という変数を作ったなら、このx成分はvec.x()で呼び出せる
const Vector vec;という定数定義でも飛び出せるようにconst用も作成

演算子オーバーロード

解説はここが分かりやすい
operator+()こういうやつを書いておくと、あとでベクトルを使った式が書けるらしい
a = b + c;という具合に

加算は

$$\mathbf{a}+\mathbf{b}=\begin{bmatrix}a_x\\a_y\\a_z\end{bmatrix}+\begin{bmatrix}b_x\\b_y\\b_z\end{bmatrix}=\begin{bmatrix}a_x+b_x\\a_y+b_y\\a_z+b_z\end{bmatrix}$$
Vector operator+(const Vector& vctr) const {
	return Vector(v[0]+vctr.v[0], v[1]+vctr.v[1], v[2]+vctr.v[2]);
}

Vector vctrはVector型の変数のこと。returnの部分は式と同じ
代入演算子も作る
*thisは演算子前の変数のポインタ
返り値がポインタなのでVector型変数のアドレス（←???）

Vector& operator+=(const Vector& vctr) {
	v[0] += vctr.v[0];
	v[1] += vctr.v[1];
	v[2] += vctr.v[2];
	return *this;
}

a += b;だったら、aを指す
足した数値を返すというより、元々あった変数に足す演算子なので、わざわざ*thisを使っている
*this = *this + vctrでも動くが、コピーの手間で遅いらしい（ChatGPT調べ）

減算も同じようにやる

ベクトルの乗算はスカラとしかできないので、

$$n\mathbf{a}=n\begin{bmatrix}a_x\\a_y\\a_z\end{bmatrix}=\begin{bmatrix}na_x\\na_y\\na_z\end{bmatrix}$$
Vector operator*(Scalar sclr) const {
	return Vector(v[0]*sclr, v[1]*sclr, v[2]*sclr);
}

演算子は左辺の変数の型を見て振る舞いを決めるので、このままだとa * nはできてもn * aができない。（ここでいう左辺右辺は演算子から見て左か右かで、等号は無関係）この場合は非メンバ関数というものを使う、と言っても引数が2個あるだけ
operator*(左辺, 右辺)と定義することで左辺がスカラでも計算できるようになる
非メンバ関数にはprivate変数にアクセスする権限がないため、friendとすることでアクセス可能にする。また、関数の処理前にconstは付けない

friend Vector operator*(Scalar sclr, const Vector& vctr) {
	return Vector(sclr*vctr.v[0], sclr*vctr.v[1], sclr*vctr.v[2]);
}

あと代入演算子

Vector& operator*=(Scalar sclr) {
	v[0] *= sclr;
	v[1] *= sclr;
	v[2] *= sclr;
	return *this;
}

除算は交換しないので非メンバ関数は要らない

OpenFOAMでは、内積は&、外積は^だが、実はC++では&も^も加減より優先度が低い（数学的には乗算なので優先度が高くなければならない）
演算子オーバーロードは演算子の処理を変えられるだけなので、優先順位までは変えられない。
なので、OpenFOAMは(a & b) + cという具合に括弧が必須


参照：wikipedia

というわけなので、ここはOpenFOAMに則らず、関数として処理していく

Scalar dot(const Vector& vctr) const{
	return (v[0] * vctr.v[0]) + (v[1] * vctr.v[1]) + (v[2] * vctr.v[2]);
}

クラス外で以下のように書く

inline Scalar dot(const Vector& vctr1, const Vector& vctr2) { return vctr1.dot(vctr2); }

a.dot(b)とも書けるし、dot(a, b)とも書ける
なんかどっちの書き方も好きだから両方残した

外積

Vector cross(const Vector& vctr) const {
	return Vector(
		v[1]*vctr.v[2]-v[2]*vctr.v[1],
		v[2]*vctr.v[0]-v[0]*vctr.v[2],
		v[0]*vctr.v[1]-v[1]*vctr.v[0]
	);
}

内積と同じように

inline Vector cross(const Vector& vctr1, const Vector& vctr2) { return vctr1.cross(vctr2); }

比較演算子
わざわざstd::array<Scalar,3> v;にしたのはこのため

bool operator==(const Vector& vctr) const { return v == vctr.v; }
bool operator!=(const Vector& vctr) const { return v != vctr.v; }

不等号は数学的にスカラじゃないと比較できない

あと標準出力用の表示も定義する

friend std::ostream& operator<<(std::ostream& os, const Vector& vctr) {
	os << "(" << vctr.v[0] << ", " << vctr.v[1] << "," << vctr.v[2] << ")";
	return os;
}

非メンバ関数を使ってシフト演算子<<が来たときに(x, y)を返すようにする
非メンバ関数なので例に漏れずfriendを付ける
また、一番上に#include <iostream>を入れる

Vector.h

#include <iostream>
#include <array>

//Scalar
using Scalar = double;
//Vector
class Vector {
    std::array<Scalar,3> v
public:
    Vector() : v{0,0,0} {}
    Vector(Scalar x, Scalar y, Scalar z) : v{x,y,z} {}
    
    Scalar& x() { return v[0]; }
    Scalar& y() { return v[1]; }
    Scalar& z() { return v[2]; }
    //const用
    const Scalar& x() const { return v[0]; }
    const Scalar& y() const { return v[1]; }
    const Scalar& z() const { return v[2]; }
    //加算
    Vector operator+(const Vector& vctr) const {
        return Vector(v[0]+vctr.v[0], v[1]+vctr.v[1], v[2]+vctr.v[2]);
    }
    Vector& operator+=(const Vector& vctr) {
        v[0] += vctr.v[0];
        v[1] += vctr.v[1];
        v[2] += vctr.v[2];
        return *this;
    }
    //減算
    Vector operator-(const Vector& vctr) const {
        return Vector(v[0]-vctr.v[0], v[1]-vctr.v[1], v[2]-vctr.v[2]);
    }
    Vector& operator-=(const Vector& vctr) {
        v[0] -= vctr.v[0];
        v[1] -= vctr.v[1];
        v[2] -= vctr.v[2];
        return *this;
    }
    //乗算
    Vector operator*(Scalar sclr) const {
        return Vector(v[0]*sclr, v[1]*sclr, v[2]*sclr);
    }
    friend Vector operator*(Scalar sclr, const Vector& vctr) {
        return Vector(sclr*vctr.v[0], sclr*vctr.v[1], sclr*vctr.v[2]);
    }
    Vector& operator*=(Scalar sclr) {
        v[0] *= sclr;
        v[1] *= sclr;
        v[2] *= sclr;
        return *this;
    }
    //除算
    Vector operator/(Scalar sclr) const {
        return Vector(v[0]/sclr, v[1]/sclr, v[2]/sclr);
    }
    Vector& operator/=(Scalar sclr) {
        v[0] /= sclr;
        v[1] /= sclr;
        v[2] /= sclr;
        return *this;
    }
    //内積
    Scalar dot(const Vector& vctr) const{
        return (v[0] * vctr.v[0]) + (v[1] * vctr.v[1]) + (v[2] * vctr.v[2]);
    }
    //外積
    Vector cross(const Vector& vctr) const {
        return Vector(
            v[1]*vctr.v[2]-v[2]*vctr.v[1],
            v[2]*vctr.v[0]-v[0]*vctr.v[2],
            v[0]*vctr.v[1]-v[1]*vctr.v[0]
        );
    }
    //比較
    bool operator==(const Vector& vctr) const { return v == vctr.v; }
    bool operator!=(const Vector& vctr) const { return v != vctr.v; }
    //標準出力用
    friend std::ostream& operator<<(std::ostream& os, const Vector& vctr) {
        os << "(" << vctr.v[0] << ", " << vctr.v[1] << ", " << vctr.v[2] << ")";
        return os;
    }
};

inline Scalar dot(const Vector& vctr1, const Vector& vctr2) { return vctr1.dot(vctr2); }
inline Vector cross(const Vector& vctr1, const Vector& vctr2) { return vctr1.cross(vctr2); }

test.cpp

#include "Vector.h"

int main(){
    Vector a(15, 5, 20);
    Vector b(2, 18, 6);
    double n=5;
    std::cout << a+b << std::endl;
    std::cout << a-b << std::endl;
    std::cout << a*n << std::endl;
    std::cout << n*a << std::endl;
    std::cout << a/n << std::endl;
    a+=b;
    std::cout << a << std::endl;
    a-=b;
    std::cout << a << std::endl;
    a*=n;
    std::cout << a << std::endl;
    a/=n;
    std::cout << a << std::endl;
    std::cout << dot(a, b) << std::endl;
    std::cout << cross(a, b) << std::endl;
    std::cout << a.norm() << std::endl;
    if (a == b) std::cout << "a=b" << std::endl;
    else std::cout << "a!=b" << std::endl;
    return 0;
}

結果

(17, 23, 26)
(13, -13, 14)
(75, 25, 100)
(75, 25, 100)
(3, 1, 4)
(17, 23, 26)
(15, 5, 20)
(75, 25, 100)
(15, 5, 20)
240
(-330, -50, 260)
25.4951
a!=b