まずmain()ですね．
#include <iostream>
#include "heat.hpp"
int main(int argc, const char * argv[]) {
int N=50,nt=0;
double dt=0.1;
Heat Solver(N,dt);
Solver.vtk(nt);
double errmax=1e-5,err;
while (errmax<(err=Solver.Step(nt++)))
{
Solver.vtk(nt);
std::cout << "nt= " << nt << " err= " << err << std::endl;
}
return 0;
}Heatクラスの仕様は上で完全に決定しているので，あとは実装するだけですわ.
#include <sstream>
#include <boost/filesystem.hpp>
#include "libSPARSE.hpp"
#include "myArray.hpp"
class Heat
{
double xmax=1.;
double u_x0=0.,u_x1=0.,u_y0=1.,u_y1=0.; //境界条件
double dx,dt,cfl;
myArray<double> mySOL;
SPARSE::Coef myMAT;
public:
Heat(int N,double dt):mySOL(N-1,N-1)
{
dx=xmax/N;
cfl=dt/dx/dx;
for(int i=0; i<mySOL.size(0);i++)
for(int j=0; j<mySOL.size(1);j++)
{
double x=(i+1)*dx,y=(j+1)*dx;
mySOL(i,j)=(1+y-2*y*y)*(1-4*(x-0.5)*(x-0.5));
}
//行列サイズは mySOL.size()xmySOL.size()になる
for(int k=0;k<mySOL.size();k++) myMAT.New(k,1.0+4*cfl);//対角成分は 1+4*cfl
for(int k=0;k<mySOL.size();k++) myMAT[k][myMAT.Dimension()]=0.;
for(int k=0;k<mySOL.size();k++)//他の要素を入れる
{
int i=mySOL.index(0,k); //k番目の未知数のi
int j=mySOL.index(1,k); //k番目の未知数のj
if (i==0) myMAT[k][myMAT.Dimension()]+=cfl*u_x0;
else myMAT[k][mySOL.serial(i-1,j)]=-cfl;
if (i==mySOL.size(0)-1) myMAT[k][myMAT.Dimension()]+=cfl*u_x1;
else myMAT[k][mySOL.serial(i+1,j)]=-cfl;
if (j==0) myMAT[k][myMAT.Dimension()]+=cfl*u_y0;
else myMAT[k][mySOL.serial(i,j-1)]=-cfl;
if (j==mySOL.size(1)-1) myMAT[k][myMAT.Dimension()]+=cfl*u_y1;
else myMAT[k][mySOL.serial(i,j+1)]=-cfl;
}
boost::filesystem::path folder=getenv("HOME");
boost::filesystem::current_path(folder/"tes_test");
}
double Step(int nt)
{
auto thisMAT=myMAT;
for(int k=0;k<mySOL.size();k++)
thisMAT[k][thisMAT.Dimension()]+=mySOL[k];
thisMAT.MakeIndex();
thisMAT.EraseUp();
thisMAT.EraseDown();
thisMAT.ClearIndex();
double err=0.0;
for(int k=0;k<mySOL.size();k++)
{
err=std::max(std::abs(thisMAT[k][thisMAT.Dimension()]-mySOL[k]),err);
mySOL[k]=thisMAT[k][thisMAT.Dimension()];
}
return err;
}
void vtk(int nt)
{
std::ostringstream oss;
oss << "vtk_" << std::setfill('0') << std::right << std::setw(6) << nt << ".vtk";
std::ofstream vout(oss.str());
vout << "# vtk DataFile Version 2.0" << std::endl << "ひゃっほう" << std::endl << "ASCII" << std::endl;
//Geometry指定
vout << "DATASET STRUCTURED_POINTS" << std::endl;
vout << "DIMENSIONS " << mySOL.size(0)+2 << " " << mySOL.size(1)+2 << " " << 1 << std::endl;
vout << "ORIGIN 0 0 0" << std::endl;
vout << "SPACING " << 1.0/(mySOL.size(0)+1) << " " << 1.0/(mySOL.size(1)+1) << " " << 1.0 << std::endl;
//DataSet指定
vout << "POINT_DATA " << (mySOL.size(0)+2)*(mySOL.size(1)+2) << std::endl;
vout << "SCALARS U double" << std::endl;
vout << "LOOKUP_TABLE default" << std::endl;
for(int j=0;j<mySOL.size(1)+2;j++)
for(int i=0;i<mySOL.size(0)+2;i++)
{
if (i==0) vout << u_x0 <<std::endl;
else if (i==mySOL.size(0)+1) vout << u_x1 <<std::endl;
else
{
if (j==0) vout << u_y0 <<std::endl;
else if (j==mySOL.size(1)+1) vout << u_y1 <<std::endl;
else vout << mySOL(i-1,j-1) <<std::endl;
}
}
}
};
#endif /* heat_hpp */新しい文法要素はありませんね．念のためmyArray.hppも示しておきます:
#ifndef myArray_hpp
#define myArray_hpp
#include <iostream>
#include <vector>
template <class T> class myArray
{
std::vector<int> _sizes;
std::vector<int> _upto;
std::vector<T> _data;
public:
T* data() {return _data.data();}
const int size(){return _upto.back();}
const int size(const int index){return _sizes.at(index);}
T& operator[](const int k){return _data[k];}
typename std::vector<T>::iterator begin(){return _data.begin();}
typename std::vector<T>::iterator end(){return _data.end();}
const int index(const int dim,const int k)
{
int _index=k;
for(int i=0;i<dim;i++) _index/=_sizes[i];
if (dim==_sizes.size()-1) return _index;
return _index%_sizes[dim];
}
myArray(const int n0,const int n1)
{
_sizes.resize(2); //2次元配列である
_upto.resize(2);
_sizes[0]=n0; _sizes[1]=n1;_upto[0]=n0; _upto[1]=_upto[0]*n1;
_data.resize(size());
}
const int serial(const int i0,const int i1){return i0+_upto[0]*i1;}
T& operator()(const int i0,const int i1){return _data[serial(i0,i1)];}
myArray(const int n0,const int n1,const int n2)
{
_sizes.resize(3); //3次元配列である
_upto.resize(3);
_sizes[0]=n0; _sizes[1]=n1; _sizes[2]=n2;_upto[0]=n0; _upto[1]=_upto[0]*n1; _upto[2]=_upto[1]*n2;
_data.resize(size());
}
const int serial(const int i0,const int i1,const int i2){return i0+_upto[0]*i1+_upto[1]*i2;}
T& operator()(const int i0,const int i1,const int i2){return _data[serial(i0,i1,i2)];}
};
#endif /* myArray_hpp */以上