全体のコード（まとめ）
あとは余計なコメントや出力などを消して完成とします。
※最終的なコードは本編から少し変えました。
「Diff1d.h」に微分の関数定義を書き、「Diff1d.cpp」に微分処理を書いていましたが、Diff1d.h内のクラス定義をやめて直接微分の処理を書くようにしました。
それによってmain.cppで以下のインスタンス化部分がいらなくなります。
1 2 3 4 5 | // delete fd1d::Diff1d diff1d; // modify Fields::vectorField1d Diff = fd1d::gradX1d(U); |
main.cpp
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 | #include <iostream> #include <vector> #include "Mesh.h" //メッシュ設定 #include "Fields.h" //場の定義 #include "Diff1d.h" //空間微分 #include "Solution1D.h" //時間ループ #include "FileWriter.h" //ファイル出力 #include "PrintVector.h" //場のprint文 using namespace std; using std::vector; int main() { int Nx =100; double totallength = 2.0; Mesh mesh_(Nx, totallength);//Meshクラスのインスタンス化 //Fieldsの設定 #include "creatFields.h"; //initial condition cout << "Initial :: Iteration completed, writing the final file" << endl; for (int i = 0; i < U.size(); i++) { U[i].value = 1.0; if ((U[i].x1d > 0.5) && (U[i].x1d < 1.0)) { U[i].value = 2.0; } } Unew = U; string U_name = "Unew"; PrintVector(Unew, U_name); //file writer string name = "1dLinearConvec"; FileWriter filewriterObject_; filewriterObject_.write1dField(Unew, name); // solution int ntimestep = 500; double dt = 0.001; Solution1d sol1d(ntimestep, dt); for (int nt = 0; nt < sol1d.nt; nt++) { cout << "Iteration no : " << nt << endl; U= Unew; double wavespeed = 1.0; double cdt = sol1d.dt * wavespeed; Fields::vectorField1d Diff = fd1d::gradX1d(U); Unew = U - (cdt * Diff); // du/dt = c*du/dx => Unew{i} = U{i} - c*dt(u{i} - u{i-1})/dx // ノイマン条件 Unew[0].value = Unew[1].value; Unew[U.size() - 1].value = Unew[U.size() - 2].value; } //file writer string name2 = "1dLinearConvec_final"; FileWriter filewriterObject_1; filewriterObject_1.write1dField(Unew, name2); return 0; } |
creatFields.hを作ってUの設定は別のヘッダーファイルで定義するようにしました。
creatFields.h
1 2 3 4 5 6 7 | // UFields(typedef vector<Fields> vectorField1d;) Fields FieldsOperations; Fields::vectorField1d U(mesh_.NumberOfNode); Fields::vectorField1d Unew(mesh_.NumberOfNode); U = FieldsOperations.passGridinfomation(U, mesh_.xpts); Unew = FieldsOperations.passGridinfomation(Unew, mesh_.xpts); |
さらにプログラム作成途中でU[i].vauleの値を確認したい場合に出力するプログラムをヘッダーファイルに書きました。
main.cppで
1 2 | //Fieldsの設定 #include "creatFields.h"; |
によりインクルードしています。
PrintVector.h
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 | #include "Fields.h" #include <iostream> #include <string> void PrintVector(Fields::vectorField1d& vec, std::string& name) { std::cout << "====== " << name << ".x1d ====== " << std::endl; for (int i = 0; i < vec.size(); i++) { std::cout << name << "[" << i << "] = " << vec[i].x1d << std::endl; } std::cout << "====== " << name << ".value ====== " << std::endl; for (int i = 0 ; i < vec.size(); i++) { std::cout << name << "[" << i << "] = " << vec[i].value << std::endl; } } |
main.cppで
1 2 | string U_name = "Unew"; PrintVector(Unew, U_name); |
の2行でよりでU[i].vauleの値を確認できます。
Mesh.h
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 | #ifndef MESH_H #define MESH_H #include "Mesh.h" #include <vector> using std::vector; class Mesh { public: Mesh(); //constructor Mesh(int&, double&); //引数ありconstructor virtual ~Mesh();//deconstructor int NumberOfNode;//x分割数 double lengthx;//x方向長さ double dx1d;//分割幅 vector<double> xpts;//x座標 void discritized1dgrid(vector<double>&); //x[i]の作成 }; #endif //MESH_H |
Mesh.cpp
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 | #include "Mesh.h" #include <iostream> Mesh::Mesh() { } Mesh::Mesh(int& n, double& txlengthx): NumberOfNode(n),//x分割数 lengthx(txlengthx),//x方向長さ xpts(NumberOfNode, 0.0)//vector(要素数, 値) { dx1d = lengthx / (NumberOfNode-1); discritized1dgrid(xpts); //add } void Mesh::discritized1dgrid(vector<double>& vec) { for (int i = 1; i < vec.size(); i++) //vec[0]=0なのでi=1からループする { vec[i] = vec[i - 1] + dx1d; } } Mesh::~Mesh() { } |
Fields.h
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 | #ifndef FIELDS_H #define FIELDS_H #include "Mesh.h" #include <vector> using std::vector; class Fields { public: Fields();//constructor virtual ~Fields(); //deconstructor typedef vector<Fields> vectorField1d; double x1d; //x座標 double value; //Uの値 Fields::vectorField1d passGridinfomation(Fields::vectorField1d&, vector<double>&); //演算子のオーバーロード friend Fields::vectorField1d operator+(const Fields::vectorField1d&, const Fields::vectorField1d&); friend Fields::vectorField1d operator-(const Fields::vectorField1d&, const Fields::vectorField1d&); friend Fields::vectorField1d operator&&(const Fields::vectorField1d&, const Fields::vectorField1d&); friend Fields::vectorField1d operator*(const double&, const Fields::vectorField1d&); friend Fields::vectorField1d operator*(const Fields::vectorField1d&, const double&); }; #endif //FIELDS_H |
Fields.cpp
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 | #include "Fields.h" #include <iostream> Fields::Fields(): x1d(0.0), value(2.0) { } Fields::vectorField1d Fields::passGridinfomation(Fields::vectorField1d& vec, vector<double>& xpts_) { /* typedef vector<Fields> vectorField1d; double value; double x1d; */ for (int i = 0; i < vec.size(); i++) { vec[i].x1d = xpts_[i]; } return vec; } Fields::vectorField1d operator+(const Fields::vectorField1d& vecA, const Fields::vectorField1d& vecB) { Fields::vectorField1d temp = vecA; for (int i = 0; i < vecA.size(); i++) { temp[i].value += vecB[i].value; } return temp; } Fields::vectorField1d operator-(const Fields::vectorField1d& vecA, const Fields::vectorField1d& vecB) { Fields::vectorField1d temp = vecA; for (int i = 0; i < vecA.size(); i++) { temp[i].value -= vecB[i].value; } return temp; } Fields::vectorField1d operator&&(const Fields::vectorField1d& vecA, const Fields::vectorField1d& vecB) { Fields::vectorField1d temp = vecA; for (int i = 0; i < vecA.size(); i++) { temp[i].value *= vecB[i].value; } return temp; } Fields::vectorField1d operator*(const double& dbvalue, const Fields::vectorField1d& vecB) { Fields::vectorField1d temp = vecB; for (int i = 0; i < vecB.size(); i++) { temp[i].value = dbvalue * vecB[i].value; } return temp; } Fields::vectorField1d operator*(const Fields::vectorField1d& vecB, const double& dbvalue) { Fields::vectorField1d temp = vecB; for (int i = 0; i < vecB.size(); i++) { temp[i].value = dbvalue * vecB[i].value; } return temp; } Fields::~Fields() { } |
Diff1d.h
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 | #include "Fields.h" #include <vector> using std::vector; namespace fd1d { // U U U U //*-----*----*----* // i-1 i Fields::vectorField1d gradX1d(Fields::vectorField1d& vec) { Fields::vectorField1d temp = vec; for (int i = 1; i < vec.size(); i++) { temp[i].value = (vec[i].value - vec[i - 1].value) / (vec[i].x1d - vec[i - 1].x1d); } return temp; } } |
Solution1D.h
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | #ifndef SOLUTION_H #define SOLUTION_H class Solution1d { public: Solution1d(int&, double&); // constructor virtual ~Solution1d(); //deconstructor int nt; //ステップ数 double dt; //時間刻み幅 }; #endif |
Solution1D.cpp
1 2 3 4 5 6 7 8 9 10 11 12 13 | #include <iostream> #include "Solution1D.h" Solution1d::Solution1d(int& nt_, double& dt_) : nt(nt_), dt(dt_) { std::cout << "Solution1d()" << std::endl; } Solution1d::~Solution1d() { } |
FileWriter.h
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 | #ifndef FILEWRITWE_H #define FILEWRITWE_H #include <stdio.h> #include <iostream> #include <vector> #include <cmath> #include <string> #include <cstring> #include <sstream> #include <fstream> //ファイルの書き込みに必要 #include "Fields.h" using namespace std; using std::vector; using std::string; class FileWriter { public: FileWriter();//constructor virtual ~FileWriter(); //destructor void write1dField(Fields::vectorField1d&, string&); }; #endif |
FileWriter.cpp
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 | #include "FileWriter.h" #include <iostream> FileWriter::FileWriter() { } void FileWriter::write1dField(Fields::vectorField1d& vec, string& name_) { string myfileType = ".dat"; //拡張子 string path = "C:\\work\\C_lang\\20220413_diff_C_puls\\\004_NS_Eq_1D_blog\\1DNS\\1DNS\\"; //絶対パス string pathName = path.append(name_).append(myfileType); //"myOutput;" + ".dat" => "myOutput.dat" std::cout << "pathName = " << pathName << " : pathName.size() = " << pathName.size() << std::endl; FILE* outfile; fopen_s(&outfile, pathName.c_str(), "w+t"); //https://www.paveway.info/entry/2018/12/25/cpp_std_string_char std::cout << "outfile = " << &outfile << std::endl; for (int i = 0; i < vec.size(); i++) { double xpos = vec[i].x1d; double uvel = vec[i].value; if (outfile) { fprintf(outfile, "%5.8lf\t%5.8lf\n", xpos, uvel); } } fclose(outfile); } FileWriter::~FileWriter() { } |
他にパラメータ、初期条件、境界条件などもプログラムの部品としてまとめておくと、後で色々な組み合わせで組みあげることができて便利になります。
今回は1次元移流方程式を対象にしましたが、拡散方程式やナビエストークス、次元数を3次元に増やしたりと条件に応じて実装できるように設計することができます。
Follow @t_kun_kamakiri