【OpenFOAM用のC++勉強(1)】ベクトル、テンソル

こんにちは。
OpenFOAM用のC++の勉強のメモを自分のために残しています。

OpenFOAM v2112

１．基本的なC++

まずはC++の基本的な構文を書いてみます。
後で、OpenFOAM用のコードに書き換えます。

#include <iostream>
#include <math.h>

using namespace std;

int main()
{
    cout << "===== integer operations =====" << endl;
    int a = 5, b = 10;
    int c = 15;

    c = a + b;

    cout << "a = " << a << endl;
    cout << "b = " << b << endl;
    cout << "a + b = " << c << endl;

    cout << "===== double precisions operations =====" << endl;
    double a_d = 12.5;
    double b_d = sqrt(a_d);
    double c_d;
    c_d = pow(a_d,2) + pow(b_d,3);
    
    cout << "a_d = " << a_d << endl;
    cout << "b_d = " << b_d << endl;
    cout << "a_d^2 + b_d^3 = " << c_d << endl;
}

g++を使ってコンパイルすると、

g++ main.cpp 

a.outという実行ファイルができます。
これを実行することで計算がされます。

 ./a.out 

===== integer operations =====
a = 5
b = 5
a + b = 15
===== double precisions operations =====
a_d = 5
b_d = 5
a_d^2 + b_d^3 = 200.444

２．OpenFOAM用に書き換え

同じことをOpenFOAMの環境に合わせて実行してみましょう。
同じ階層にMakeというフォルダを作り、その中に
・files
・options
というファイルを作成します。

files

main.cpp

EXE = $(FOAM_USER_APPBIN)/basic001

echo $FOAM_USER_APPBIN 

で確認すると以下のフォルダにbasic001というコマンドを作成します。
/home/kamakiri/OpenFOAM/kamakiri-v2012/platforms/linux64Gcc63DPInt32Opt/bin

options
何も書かない
ファイル構成は以下のようになっています。

.
├── Make
│ ├── files
│ └── options
├── a.out
└── mainOF.C

#include "Field.H"

using namespace Foam;

int main()
{
    Info << "===== integer operations =====" << endl;
    label a = 5, b = 10;
    label c = 15;

    c = a + b;

    Info << "a = " << a << endl;
    Info << "b = " << b << endl;
    Info << "a + b = " << c << endl;

    Info << "===== double precisions operations =====" << endl;
    scalar a_d = 12.5;
    scalar b_d = Foam::sqrt(a_d);
    scalar c_d;
    c_d = pow(a_d,2) + pow(b_d,3);
    
    Info << "a_d = " << a_d << endl;
    Info << "b_d = " << b_d << endl;
    Info << "a_d^2 + b_d^3 = " << c_d << endl;

    return 0;
}

以下のように置き換えました。
int ➡ label
double ➡ scalar
cout ➡ Info
その他にも後ほどでてくるvectorやtensorや文字列を定義するwordなどがあります。

label a;
label a(10);
label a = 10;
scalar b;
scalar b(1.5);
scalar b = 1.5;
vector v1(1, 2, 3);
tensor t1(1, 2, 3, 4, 5, 6, 7, 8, 9);
word text = "OpenFOAM"; //文字列を定義
List<vector> a(10); // 配列を定義
a[0] = vector(1, 2, 3);

では、先ほどのコードをwmakeを実行してコンパイルします。

wmake

コンパイルすると/home/kamakiri/OpenFOAM/kamakiri-v2012/platforms/linux64Gcc63DPInt32Opt/bin
にbasic001という実行ファイルができるのでコマンドで実行しします。

basic001

結果さ先ほどと同じになります。

３．配列を定義

配列もやってみましょう。「１．」のコードの続きに書いていきます。

#include <iostream>
#include <math.h>

using namespace std;

int main()
{
    cout << "===== integer operations =====" << endl;
    int a = 5, b = 10;
    int c = 15;

    c = a + b;

    cout << "a = " << a << endl;
    cout << "b = " << b << endl;
    cout << "a + b = " << c << endl;

    cout << "===== double precisions operations =====" << endl;
    double a_d = 12.5;
    double b_d = sqrt(a_d);
    double c_d;
    c_d = pow(a_d,2) + pow(b_d,3);
    
    cout << "a_d = " << a_d << endl;
    cout << "b_d = " << b_d << endl;
    cout << "a_d^2 + b_d^3 = " << c_d << endl;

    cout << "===== vector (array) operatios =====" << endl;
    double v1[3] = {1, 2, 3};
    
    cout << "total v1 size is " << sizeof(v1) << endl;
    cout << "double size is " << sizeof(double) << endl;
    cout << "v1 element size is " << sizeof(v1)/sizeof(double) << endl;

    double v2[3];
    v2[0] = 10;
    v2[1] = 20;
    v2[2] = 30;

    double v_pro;

    for (int i=0; i< sizeof(v1)/sizeof(double); i++)
    {
        cout << "v1[" << i << "] = " << v1[i] << endl; 
    }
    cout << endl;

    for (int i=0; i<sizeof(v2)/sizeof(double); i++)
    {
        cout << "v2[" << i << "] = " << v2[i] << endl; 
    }
    cout << endl;

    for (int i=0; i<sizeof(v2)/sizeof(double); i++)
    {
        cout << "v1[" << i << "] * " << "v2[" << i << "] = "
        << v1[i] * v2[i] << endl; 
        v_pro += v1[i] * v2[i];
    }
    cout << "v1 * v2 = " << v_pro << endl; 

    return 0;
}

wmakeコンパイルします。
※OpenFOAM用のコードは何も使っていないのでg++で良いですが・・・
optionsでのbasic003として実行ファイルを作成します。
basic003を実行。

===== integer operations =====
a = 5
b = 10
a + b = 15
===== double precisions operations =====
a_d = 12.5
b_d = 3.53553
a_d^2 + b_d^3 = 200.444
===== vector (array) operatios =====
total v1 size is 24
double size is 8
v1 element size is 3
v1[0] = 1
v1[1] = 2
v1[2] = 3

v2[0] = 10
v2[1] = 20
v2[2] = 30

v1[0] * v2[0] = 10
v1[1] * v2[1] = 40
v1[2] * v2[2] = 90
v1 * v2 = 140

物理現象を解くにあたってベクトルやテンソルなど配列要素を使うことは必須ですが、自分で配列を定義したら負けというのがあります。
C++のvectorコンテナを使うのが効率が良いのですが、OpenFOAMではさらにベクトルとテンソルを定義しなおして使えるようになっています。

４．ベクトル定義をOpenFOAM用に書き換える

以上のC++のコードをOpenFOAM用のコードに置き換えていきます。

OpenFOAM: API Guide: src/finiteVolume/cfdTools/general/include/fvCFD.H Source File

#include "fvCFD.H"

using namespace Foam;

int main()
{
    Info << "===== integer operations =====" << endl;
    label a = 5, b = 10;
    label c = 15;

    c = a + b;

    Info << "a = " << a << endl;
    Info << "b = " << b << endl;
    Info << "a + b = " << c << endl;

    Info << "===== double precisions operations =====" << endl;
    scalar a_d = 12.5;
    scalar b_d = Foam::sqrt(a_d);
    scalar c_d;
    c_d = pow(a_d,2) + pow(b_d,3);
    
    Info << "a_d = " << a_d << endl;
    Info << "b_d = " << b_d << endl;
    Info << "a_d^2 + b_d^3 = " << c_d << endl;

    Info << "===== vector (array) operatios =====" << endl;
    vector v1(1, 2, 3);
    vector v2(10, 20, 30);

    scalar v_pro;

    forAll(v1, i)
    {
        Info << "v1[" << i << "] = " << v1[i] << endl; 
    }
    Info << endl;

    forAll(v2, i)
    {
        Info << "v2[" << i << "] = " << v2[i] << endl; 
    }
    Info << endl;

    v_pro = v1 & v2;
    
    Info << "v1 * v2 = " << v_pro << endl; 

    return 0;
}

files

main.C

EXE = $(FOAM_USER_APPBIN)/basic004

fvCFD.Hを読み込むためのパスを指定します。

options

EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude

EXE_LIBS = \
    -lfiniteVolume \
    -lmeshTools

wmakeを実行してコンパイルしbasic004を実行します。

===== integer operations =====
a = 5
b = 10
a + b = 15
===== double precisions operations =====
a_d = 12.5
b_d = 3.53553
a_d^2 + b_d^3 = 200.444
===== vector (array) operatios =====
v0[0] = 1
v0[1] = 2
v0[2] = 3
v0[3] = 4

v1[0] = 1
v1[1] = 2
v1[2] = 3

v2[0] = 10
v2[1] = 20
v2[2] = 30

v1 * v2 = 140

５．テンソルをOpenFOAMのクラスで定義

ベクトルに加えてテンソルも定義します。

#include "Field.H"

using namespace Foam;

int main()
{
    Info << "===== integer operations =====" << endl;
    label a = 5, b = 10;
    label c = 15;

    c = a + b;

    Info << "a = " << a << endl;
    Info << "b = " << b << endl;
    Info << "a + b = " << c << endl;

    Info << "===== double precisions operations =====" << endl;
    scalar a_d = 12.5;
    scalar b_d = Foam::sqrt(a_d);
    scalar c_d;
    c_d = pow(a_d,2) + pow(b_d,3);
    
    Info << "a_d = " << a_d << endl;
    Info << "b_d = " << b_d << endl;
    Info << "a_d^2 + b_d^3 = " << c_d << endl;

    Info << "===== vector (array) operatios =====" << endl;
    vector v1(1, 2, 3);
    vector v2(10, 20, 30);

    forAll(v1, i)
    {
        Info << "v1[" << i << "] = " << v1[i] << endl; 
    }
    Info << endl;

    forAll(v2, i)
    {
        Info << "v2[" << i << "] = " << v2[i] << endl; 
    }
    Info << endl;

    scalar v_pro;
    v_pro = v1 & v2;
    Info << "scalar v_pro = " << v_pro << endl; 
    Info << "v1 * v2 = " << v1 * v2 << endl; 

    Info << "===== Fields operatios =====" << endl;
    Field<label> fvector(3,0);
    forAll(fvector, i)
    {
        Info << "fvector[" << i << "] = " << fvector[i] << endl; 
    }

    Info << "===== tensor operatios =====" << endl;
    tensor t1(1, 2, 3, 4, 5, 6, 7 ,8, 9);
    tensor t2(10, 20, 30, 40, 50, 60, 70 , 80, 90);
    tensor t3(tensor::zero);
    Info << "tensor t3 = " << t1 << endl;
    Info << "Transpose t2 = " << t2.T() << endl;
    
    tensor t4;
    t4 = t1 & t2;
    Info << "t1 & t2 = " << t4 << endl;

    return 0;

}

wmakeコンパイルします。
optionsでのbasic005として実行ファイルを作成します。
basic005を実行。

===== integer operations =====
a = 5
b = 10
a + b = 15
===== double precisions operations =====
a_d = 12.5
b_d = 3.53553
a_d^2 + b_d^3 = 200.444
===== vector (array) operatios =====
v1[0] = 1
v1[1] = 2
v1[2] = 3

v2[0] = 10
v2[1] = 20
v2[2] = 30

scalar v_pro = 140
v1 * v2 = (10 20 30 20 40 60 30 60 90)
===== Fields operatios =====
fvector[0] = 0
fvector[1] = 0
fvector[2] = 0
===== tensor operatios =====
tensor t3 = (1 2 3 4 5 6 7 8 9)
Transpose t2 = (10 40 70 20 50 80 30 60 90)
t1 & t2 = (300 360 420 660 810 960 1020 1260 1500)

OpenFOAM: API Guide: src/OpenFOAM/fields/Fields/fieldTypes.H File Reference

６．テンソルから固有値と固有ベクトルを計算

テンソルの固有値や固有値ベクトルを求めるには以下のtensor.Hで定義された関数を使います。

OpenFOAM: API Guide: src/OpenFOAM/primitives/Tensor/floats/tensor.H Source File

#include "Field.H"

using namespace Foam;

int main()
{
    Info << "===== integer operations =====" << endl;
    label a = 5, b = 10;
    label c = 15;

    c = a + b;

    Info << "a = " << a << endl;
    Info << "b = " << b << endl;
    Info << "a + b = " << c << endl;

    Info << "===== double precisions operations =====" << endl;
    scalar a_d = 12.5;
    scalar b_d = Foam::sqrt(a_d);
    scalar c_d;
    c_d = pow(a_d,2) + pow(b_d,3);
    
    Info << "a_d = " << a_d << endl;
    Info << "b_d = " << b_d << endl;
    Info << "a_d^2 + b_d^3 = " << c_d << endl;

    Info << "===== vector (array) operatios =====" << endl;
    vector v1(1, 2, 3);
    vector v2(10, 20, 30);

    forAll(v1, i)
    {
        Info << "v1[" << i << "] = " << v1[i] << endl; 
    }
    Info << endl;

    forAll(v2, i)
    {
        Info << "v2[" << i << "] = " << v2[i] << endl; 
    }
    Info << endl;

    scalar v_pro;
    v_pro = v1 & v2;
    Info << "scalar v_pro = " << v_pro << endl; 
    Info << "v1 * v2 = " << v1 * v2 << endl; 

    Info << "===== Fields operatios =====" << endl;
    Field<label> fvector(3,0);
    forAll(fvector, i)
    {
        Info << "fvector[" << i << "] = " << fvector[i] << endl; 
    }

    Info << "===== tensor operatios =====" << endl;
    tensor t1(1, 2, 3, 4, 5, 6, 7 ,8, 9);
    tensor t2(10, 20, 30, 40, 50, 60, 70 , 80, 90);
    tensor t3(tensor::zero);
    Info << "tensor t3 = " << t1 << endl;
    Info << "Transpose t2 = " << t2.T() << endl;
    Info << "t2 xz = " << t2.xz() << endl;

    // tensor.H
    // https://www.openfoam.com/documentation/guides/latest/api/tensor_8H_source.html
    Info << "t2 engin value = " << eigenValues(t2) << endl;
    Info << "t2 engin vectors = " << eigenVectors(t2) << endl;

    tensor t4;
    t4 = t1 & t2;
    Info << "t1 & t2 = " << t4 << endl;

    return 0;
}

wmakeを実行してコンパイルします。
コンパイルするとbasic006という実行ファイルができますので実行すると以下のようになります。

===== integer operations =====
a = 5
b = 10
a + b = 15
===== double precisions operations =====
a_d = 12.5
b_d = 3.53553
a_d^2 + b_d^3 = 200.444
===== vector (array) operatios =====
v1[0] = 1
v1[1] = 2
v1[2] = 3

v2[0] = 10
v2[1] = 20
v2[2] = 30

scalar v_pro = 140
v1 * v2 = (10 20 30 20 40 60 30 60 90)
===== Fields operatios =====
fvector[0] = 0
fvector[1] = 0
fvector[2] = 0
===== tensor operatios =====
tensor t3 = (1 2 3 4 5 6 7 8 9)
Transpose t2 = (10 40 70 20 50 80 30 60 90)
t2 xz = 30
t2 engin value = ((161.168 0) (-11.1684 0) (0 0))
t2 engin vectors = ((0.231971 0) (0.525322 0) (0.818673 0) (0.78583 0) (0.0867513 0) (-0.612328 0) (-0.408248 -0) (0.816497 0) (-0.408248 -0))
t1 & t2 = (300 360 420 660 810 960 1020 1260 1500)

Field.Hしかインクルードしていないのにどうしてテンソルのクラスなどが使えるのかと不思議に思いましたが、Field.Hを探っているとTensorクラスを継承していることがわかりました。

OpenFOAM: API Guide: src/OpenFOAM/fields/Fields/Field/Field.H Source File

継承関係を追っていきます。

継承関係を追っていくとtensor.Hを継承していることがわかります。
さらにtensor.HはTensor.Hを継承しています。

行列の行列式と逆行列

行列には行列式や逆行列などがあります。

//#include "fvCFD.H"
//#include "fieldTypes.H"	
#include "Field.H"
// #include "tensor.H"

using namespace Foam;

int main()
{
    Info << "===== integer operations =====" << endl;
    label a = 5, b = 10;
    label c = 15;

    c = a + b;

    Info << "a = " << a << endl;
    Info << "b = " << b << endl;
    Info << "a + b = " << c << endl;

    Info << "===== double precisions operations =====" << endl;
    scalar a_d = 12.5;
    scalar b_d = Foam::sqrt(a_d);
    scalar c_d;
    c_d = pow(a_d,2) + pow(b_d,3);
    
    Info << "a_d = " << a_d << endl;
    Info << "b_d = " << b_d << endl;
    Info << "a_d^2 + b_d^3 = " << c_d << endl;

    Info << "===== vector (array) operatios =====" << endl;
    vector v1(1, 2, 3);
    vector v2(10, 20, 30);

    forAll(v1, i)
    {
        Info << "v1[" << i << "] = " << v1[i] << endl; 
    }
    Info << endl;

    forAll(v2, i)
    {
        Info << "v2[" << i << "] = " << v2[i] << endl; 
    }
    Info << endl;

    scalar v_pro;
    v_pro = v1 & v2;
    Info << "scalar v_pro = " << v_pro << endl; 
    Info << "v1 * v2 = " << v1 * v2 << endl; 

    Info << "===== Fields operatios =====" << endl;
    Field<label> fvector(3,0);
    forAll(fvector, i)
    {
        Info << "fvector[" << i << "] = " << fvector[i] << endl; 
    }

    Info << "===== tensor operatios =====" << endl;
    tensor t1(1, 2, 3, 4, 5, 6, 7 ,8, 9);
    tensor t2(1, 1, -1, -2, -1, 1, -1, -2, 1);
    tensor t3(tensor::zero);
    Info << "tensor t3 = " << t1 << endl;
    Info << "Transpose t2 = " << t2.T() << endl;
    Info << "t2 xz = " << t2.xz() << endl;

    // tensor.H
    // https://www.openfoam.com/documentation/guides/latest/api/tensor_8H_source.html
    Vector< complex >  e(eigenValues(t2));
    Tensor< complex > ev = eigenVectors(t2);
    //vector  ed(eigenValues(t2)); // error
    //tensor edv = eigenVectors(t2); // error

    Info << "t2 engin value = " << e << endl;
    Info << "t2 engin vectors = " << ev << endl;
    //Info << "t2 engin value = " << ed << endl;
    //Info << "t2 engin vectors = " << edv << endl;

    tensor t4;
    t4 = t1 & t2;
    Info << "t1 & t2 = " << t4 << endl;

    Info << "e.x = " << eigenValues(t2).x() << endl;
    Info << "e.y = " << eigenValues(t2).y() << endl;
    Info << "e.z = " << eigenValues(t2).z() << endl;
    // TentorI.H
    //https://www.openfoam.com/documentation/guides/latest/api/TensorI_8H_source.html
    Info << "determinant of tensor det(t2) = " << det(t2) << endl;
    tensor inv_t2( inv(t2)); // Tensor< Cmpt > 	inv (const Tensor< Cmpt > &t)
    Info << "inverce of tensor t2 = " << inv_t2 << endl;

    return 0;
}

eigenValues(t2)で返ってくる型は3次元ベクトルで要素が複素数である「Vector< complex >」なので、

Vector< complex >  e(eigenValues(t2));
Vector< complex >  e(eigenValues(t2));

などと定義します。
以下で定義してみるとエラーがでます。

vector e(eigenValues(t2));

実は、vectorは「vector.H」で以下のように定義されています。

OpenFOAM: API Guide: src/OpenFOAM/primitives/Vector/floats/vector.H File Reference

vector.H

vectorはVector<scalar>と定義されているため、scalarを持つベクトルなのというわけです。eigenValues(t2)で返ってくる型は3次元ベクトルで要素が複素数なのでVector< complex >としなければいけないというわけです。

wmakeコンパイルします。
optionsでのbasic005として実行ファイルを作成します。
basic005を実行。

===== integer operations =====
a = 5
b = 10
a + b = 15
===== double precisions operations =====
a_d = 12.5
b_d = 3.53553
a_d^2 + b_d^3 = 200.444
===== vector (array) operatios =====
v1[0] = 1
v1[1] = 2
v1[2] = 3

v2[0] = 10
v2[1] = 20
v2[2] = 30

scalar v_pro = 140
v1 * v2 = (10 20 30 20 40 60 30 60 90)
===== Fields operatios =====
fvector[0] = 0
fvector[1] = 0
fvector[2] = 0
===== tensor operatios =====
tensor t3 = (1 2 3 4 5 6 7 8 9)
Transpose t2 = (1 -2 -1 1 -1 -2 -1 1 1)
t2 xz = -1
t2 engin value = ((-0.392647 0) (0.696323 -1.43595) (0.696323 1.43595))
t2 engin vectors = ((0.288278 0) (0.445913 0) (0.847383 0) (-0.42447 0.18794) (0.646488 0) (0.247714 -0.552444) (-0.42447 -0.18794) (0.646488 0) (0.247714 0.552444))
t1 & t2 = (-6 -7 4 -12 -13 7 -18 -19 10)
e.x = (-0.392647 0)
e.y = (0.696323 -1.43595)
e.z = (0.696323 1.43595)
determinant of tensor det(t2) = -1
inverce of tensor t2 = (-1 -1 -0 -1 -0 -1 -3 -1 -1)

理論と合っているか確かめてみます。

引き続きC++の勉強を続けます。

プログラミング言語C++ 第4版

Effective Modern C++ ―C++11/14プログラムを進化させる42項目

Twitter➡@t_kun_kamakiri
ブログ➡宇宙に入ったカマキリ（物理ブログ）