NUMCXX  0.13.20181108
Numerical library for small projects and teaching purposes
42-convtest-fem-sin.cxx

Stationary convergence test for fem method

#include <cstdio>
#include <iostream>
#include <ctime>
#include <cmath>
#include <numcxx/fem2d.hxx>
#ifdef VTKFIG
#include "vtkfigFrame.h"
#include "vtkfigDataSet.h"
#include "vtkfigGridView.h"
#include "vtkfigScalarView.h"
#include "vtkfigXYPlot.h"
#endif
int main(void)
{
double maxref=9;
numcxx::Geometry Geometry;
Geometry.set_points({
{0,0},
{1,0},
{1,1},
{0,1},
{0.7,0.7},
});
Geometry.set_bfaces({
{0,1},
{1,2},
{2,3},
{3,0}
});
Geometry.set_bfaceregions({1,2,3,4});
Geometry.set_regionpoints({
{0.5,0.5}
});
Geometry.set_regionnumbers({1});
auto frame=vtkfig::Frame::New();
frame->SetSize(800,800);
frame->SetLayout(2,2);
auto gridview=vtkfig::GridView::New();
frame->AddFigure(gridview,0);
auto solview=vtkfig::ScalarView::New();
frame->AddFigure(solview,1);
auto errplot=vtkfig::XYPlot::New();
frame->AddFigure(errplot,2);
errplot->SetXTitle("log10(h)");
errplot->SetYTitle("log10(error)");
auto timeplot=vtkfig::XYPlot::New();
frame->AddFigure(timeplot,3);
timeplot->SetXTitle("log10(N)");
timeplot->SetYTitle("log10(time)");
std::vector<double> H;
std::vector<double> N;
std::vector<double> TAsm;
std::vector<double> TLuf;
std::vector<double> TLus;
std::vector<double> TGen;;
std::vector<double> L2Error;
std::vector<double> H1Error;
for (int iref=0;iref<maxref;iref++)
{
double t0;
double h_intended=0.5*pow(2.0,-iref);
double vol= h_intended*h_intended*0.25;
Geometry.set_regionvolumes({vol});
numcxx::SimpleGrid grid(Geometry,"zpaAqDV");
double t_gen=numcxx::cpu_clock()-t0;
double hmin,hmax;
grid.calc_hminmax(hmin,hmax);
printf("h_intended: %8.3g h_max: %8.3g\n",h_intended,hmax);
numcxx::DArray1 bcfac(6);
numcxx::DArray1 bcval(6);
bcfac=0;
bcval=0;
bcfac(1)=fem2d::Dirichlet;
bcfac(2)=fem2d::Dirichlet;
bcfac(3)=fem2d::Dirichlet;
bcfac(4)=fem2d::Dirichlet;
bcval(1)=0.0;
bcval(2)=0.0;
bcval(3)=0.0;
bcval(4)=0.0;
auto nnodes=grid.npoints();
auto nodes=grid.get_points();
numcxx::DArray1 source(nnodes);
numcxx::DArray1 exact(nnodes);
numcxx::DArray1 kappa(nnodes);
kappa=1;
for (int i=0; i<nnodes;i++)
{
double x=nodes(i,0);
double y=nodes(i,1);
exact[i]=sin(M_PI*x)*sin(M_PI*y);
source[i]=2.0*M_PI*M_PI*exact[i];
}
numcxx::DSparseMatrix SGlobal(nnodes,nnodes);
numcxx::DArray1 Rhs(nnodes);
numcxx::DArray1 Sol(nnodes);
fem2d::assemble_heat_problem(grid,bcfac,bcval,source,kappa,SGlobal, Rhs);
double t_asm=numcxx::cpu_clock()-t0;
numcxx::DSolverUMFPACK Solver(SGlobal);
Solver.update(SGlobal);
double t_luf=numcxx::cpu_clock()-t0;
Solver.solve(Sol,Rhs);
double t_lus=numcxx::cpu_clock()-t0;
double l2error=fem2d::l2norm(grid,Sol-exact);
double h1error=fem2d::h1norm(grid,Sol-exact);
printf("h: %8.3e l2: %8.3e h1: %8.3e\n",hmax,l2error,h1error);
H.push_back(log10(hmax));
H1Error.push_back(log10(h1error));
L2Error.push_back(log10(l2error));
printf("time/ms gen: %8.3f asm: %8.3f luf: %8.3f lus: %8.3f\n",
t_gen*1000.0, t_asm*1000.0, t_luf*1000.0, t_lus*1000.0);
N.push_back(log10(grid.npoints()));;
TAsm.push_back(log10(t_asm));
TLuf.push_back(log10(t_luf));
TLus.push_back(log10(t_lus));
TGen.push_back(log10(t_gen));
#ifdef VTKFIG
auto griddata=numcxx::vtkfigDataSet(grid);
griddata->SetPointScalar(Sol ,"Sol");
gridview->SetData(griddata);
solview->SetData(griddata,"Sol");
frame->SetFrameTitle("iref="+std::to_string(iref));
if (iref>0)
{
double x0=-3;
double x1=0;
auto X=numcxx::DArray1{x0,x1};
auto OH1=numcxx::DArray1{x0,x1};
auto OH2=numcxx::DArray1{2*x0,2*x1};
errplot->Clear();
errplot->SetLegendPosition(0.7,0.4);
errplot->SetXRange(x0,x1);
errplot->SetYRange(-6,0);
errplot->SetPlotLineType("-");
errplot->SetPlotMarkerType("o");
errplot->SetMarkerSize(2);
errplot->SetPlotColor(0.2,0.5,0.2);
errplot->SetPlotLegend("O(h**2)");
errplot->AddPlot(X,OH2);
errplot->SetPlotColor(0.5,0.2,0.2);
errplot->SetPlotLegend("O(h)");
errplot->AddPlot(X,OH1);
errplot->SetPlotColor(0.0,0.5,0);
errplot->SetPlotLegend("L2");
errplot->AddPlot(H,L2Error);
errplot->SetPlotLegend("H1");
errplot->SetPlotColor(0.5,0,0);
errplot->AddPlot(H,H1Error);
double n0=0;
double n1=7;
double n11=6;
double t0=-5;
double t1=2;
auto XN=numcxx::DArray1{n0,n1};
auto XN1=numcxx::DArray1{n0,n11};
auto ON1=numcxx::DArray1{n0-6.0,n1-6.0};
auto ON32=numcxx::DArray1{1.5*n0-8.0, 1.5*n11-8.0};
timeplot->SetLegendPosition(0.3,0.6);
timeplot->SetLegendSize(0.15,0.3);
timeplot->SetXRange(n0,n1);
timeplot->SetYRange(t0,t1);
timeplot->Clear();
timeplot->SetPlotLineType("-");
timeplot->SetPlotMarkerType("o");
timeplot->SetMarkerSize(2);
timeplot->SetPlotLegend("Gen");
timeplot->SetPlotColor(0.5,0,0);
timeplot->AddPlot(N,TGen);
timeplot->SetPlotLegend("Asm");
timeplot->SetPlotColor(0.0,0.5,0);
timeplot->AddPlot(N,TAsm);
timeplot->SetPlotLegend("LU fact");
timeplot->SetPlotColor(0.0,0.0,0.5);
timeplot->AddPlot(N,TLuf);
timeplot->SetPlotLegend("LU solve");
timeplot->SetPlotColor(0.25,0.25,0);
timeplot->AddPlot(N,TLus);
timeplot->SetPlotLegend("O(N**(3/2))");
timeplot->SetPlotColor(0.0,0.0,1.0);
timeplot->AddPlot(XN1,ON32);
timeplot->SetPlotLegend("O(N)");
timeplot->SetPlotColor(1.0,0,0);
timeplot->AddPlot(XN,ON1);;
}
frame->Interact();
#endif
}
}