NUMCXX  0.13.20181108
Numerical library for small projects and teaching purposes
42-convtest-fem-sin.cxx
Go to the documentation of this file.
1 #include <cstdio>
7 #include <iostream>
8 #include <ctime>
9 #include <cmath>
10 #include <numcxx/numcxx.hxx>
11 #include <numcxx/simplegrid.hxx>
12 #include <numcxx/fem2d.hxx>
13 #ifdef VTKFIG
15 
16 #include "vtkfigFrame.h"
17 #include "vtkfigDataSet.h"
18 #include "vtkfigGridView.h"
19 #include "vtkfigScalarView.h"
20 #include "vtkfigXYPlot.h"
21 #endif
22 
23 int main(void)
24 {
25  double maxref=9;
26  numcxx::Geometry Geometry;
27  Geometry.set_points({
28  {0,0},
29  {1,0},
30  {1,1},
31  {0,1},
32  {0.7,0.7},
33  });
34 
35  Geometry.set_bfaces({
36  {0,1},
37  {1,2},
38  {2,3},
39  {3,0}
40  });
41 
42 
43  Geometry.set_bfaceregions({1,2,3,4});
44 
45  Geometry.set_regionpoints({
46  {0.5,0.5}
47  });
48  Geometry.set_regionnumbers({1});
49 
50 
51  auto frame=vtkfig::Frame::New();
52  frame->SetSize(800,800);
53  frame->SetLayout(2,2);
54 
55  auto gridview=vtkfig::GridView::New();
56  frame->AddFigure(gridview,0);
57 
58  auto solview=vtkfig::ScalarView::New();
59  frame->AddFigure(solview,1);
60 
61  auto errplot=vtkfig::XYPlot::New();
62  frame->AddFigure(errplot,2);
63  errplot->SetXTitle("log10(h)");
64  errplot->SetYTitle("log10(error)");
65 
66  auto timeplot=vtkfig::XYPlot::New();
67  frame->AddFigure(timeplot,3);
68  timeplot->SetXTitle("log10(N)");
69  timeplot->SetYTitle("log10(time)");
70 
71 
72 
73 
74  std::vector<double> H;
75  std::vector<double> N;
76  std::vector<double> TAsm;
77  std::vector<double> TLuf;
78  std::vector<double> TLus;
79  std::vector<double> TGen;;
80 
81 
82  std::vector<double> L2Error;
83  std::vector<double> H1Error;
84 
85 
86  for (int iref=0;iref<maxref;iref++)
87  {
88  double t0;
89  double h_intended=0.5*pow(2.0,-iref);
90  double vol= h_intended*h_intended*0.25;
91 
92  Geometry.set_regionvolumes({vol});
93  t0=numcxx::cpu_clock();
94  numcxx::SimpleGrid grid(Geometry,"zpaAqDV");
95  double t_gen=numcxx::cpu_clock()-t0;
96 
97  double hmin,hmax;
98  grid.calc_hminmax(hmin,hmax);
99  printf("h_intended: %8.3g h_max: %8.3g\n",h_intended,hmax);
100 
101 
102  numcxx::DArray1 bcfac(6);
103  numcxx::DArray1 bcval(6);
104  bcfac=0;
105  bcval=0;
106 
107  bcfac(1)=fem2d::Dirichlet;
108  bcfac(2)=fem2d::Dirichlet;
109  bcfac(3)=fem2d::Dirichlet;
110  bcfac(4)=fem2d::Dirichlet;
111  bcval(1)=0.0;
112  bcval(2)=0.0;
113  bcval(3)=0.0;
114  bcval(4)=0.0;
115 
116  auto nnodes=grid.npoints();
117  auto nodes=grid.get_points();
118 
119  numcxx::DArray1 source(nnodes);
120  numcxx::DArray1 exact(nnodes);
121  numcxx::DArray1 kappa(nnodes);
122  kappa=1;
123 
124  for (int i=0; i<nnodes;i++)
125  {
126  double x=nodes(i,0);
127  double y=nodes(i,1);
128  exact[i]=sin(M_PI*x)*sin(M_PI*y);
129  source[i]=2.0*M_PI*M_PI*exact[i];
130  }
131 
132  numcxx::DSparseMatrix SGlobal(nnodes,nnodes);
133  numcxx::DArray1 Rhs(nnodes);
134  numcxx::DArray1 Sol(nnodes);
135 
136  t0=numcxx::cpu_clock();
137  fem2d::assemble_heat_problem(grid,bcfac,bcval,source,kappa,SGlobal, Rhs);
138  double t_asm=numcxx::cpu_clock()-t0;
139 
140  numcxx::DSolverUMFPACK Solver(SGlobal);
141 
142  t0=numcxx::cpu_clock();
143  Solver.update(SGlobal);
144  double t_luf=numcxx::cpu_clock()-t0;
145 
146  t0=numcxx::cpu_clock();
147  Solver.solve(Sol,Rhs);
148  double t_lus=numcxx::cpu_clock()-t0;
149 
150  double l2error=fem2d::l2norm(grid,Sol-exact);
151  double h1error=fem2d::h1norm(grid,Sol-exact);
152 
153  printf("h: %8.3e l2: %8.3e h1: %8.3e\n",hmax,l2error,h1error);
154  H.push_back(log10(hmax));
155  H1Error.push_back(log10(h1error));
156  L2Error.push_back(log10(l2error));
157 
158  printf("time/ms gen: %8.3f asm: %8.3f luf: %8.3f lus: %8.3f\n",
159  t_gen*1000.0, t_asm*1000.0, t_luf*1000.0, t_lus*1000.0);
160 
161  N.push_back(log10(grid.npoints()));;
162  TAsm.push_back(log10(t_asm));
163  TLuf.push_back(log10(t_luf));
164  TLus.push_back(log10(t_lus));
165  TGen.push_back(log10(t_gen));
166 
167 
168 #ifdef VTKFIG
169  auto griddata=numcxx::vtkfigDataSet(grid);
170  griddata->SetPointScalar(Sol ,"Sol");
171 
172  gridview->SetData(griddata);
173  solview->SetData(griddata,"Sol");
174 
175  frame->SetFrameTitle("iref="+std::to_string(iref));
176  if (iref>0)
177  {
178  double x0=-3;
179  double x1=0;
180  auto X=numcxx::DArray1{x0,x1};
181  auto OH1=numcxx::DArray1{x0,x1};
182  auto OH2=numcxx::DArray1{2*x0,2*x1};
183 
184  errplot->Clear();
185 
186 
187  errplot->SetLegendPosition(0.7,0.4);
188  errplot->SetXRange(x0,x1);
189  errplot->SetYRange(-6,0);
190 
191 
192  errplot->SetPlotLineType("-");
193  errplot->SetPlotMarkerType("o");
194  errplot->SetMarkerSize(2);
195 
196 
197  errplot->SetPlotColor(0.2,0.5,0.2);
198  errplot->SetPlotLegend("O(h**2)");
199  errplot->AddPlot(X,OH2);
200  errplot->SetPlotColor(0.5,0.2,0.2);
201  errplot->SetPlotLegend("O(h)");
202  errplot->AddPlot(X,OH1);
203 
204  errplot->SetPlotColor(0.0,0.5,0);
205  errplot->SetPlotLegend("L2");
206 
207  errplot->AddPlot(H,L2Error);
208  errplot->SetPlotLegend("H1");
209  errplot->SetPlotColor(0.5,0,0);
210  errplot->AddPlot(H,H1Error);
211 
212 
213  double n0=0;
214  double n1=7;
215  double n11=6;
216  double t0=-5;
217  double t1=2;
218 
219  auto XN=numcxx::DArray1{n0,n1};
220  auto XN1=numcxx::DArray1{n0,n11};
221  auto ON1=numcxx::DArray1{n0-6.0,n1-6.0};
222  auto ON32=numcxx::DArray1{1.5*n0-8.0, 1.5*n11-8.0};
223 
224 
225  timeplot->SetLegendPosition(0.3,0.6);
226  timeplot->SetLegendSize(0.15,0.3);
227  timeplot->SetXRange(n0,n1);
228  timeplot->SetYRange(t0,t1);
229 
230 
231  timeplot->Clear();
232  timeplot->SetPlotLineType("-");
233  timeplot->SetPlotMarkerType("o");
234  timeplot->SetMarkerSize(2);
235 
236 
237  timeplot->SetPlotLegend("Gen");
238  timeplot->SetPlotColor(0.5,0,0);
239  timeplot->AddPlot(N,TGen);
240 
241  timeplot->SetPlotLegend("Asm");
242  timeplot->SetPlotColor(0.0,0.5,0);
243  timeplot->AddPlot(N,TAsm);
244 
245  timeplot->SetPlotLegend("LU fact");
246  timeplot->SetPlotColor(0.0,0.0,0.5);
247  timeplot->AddPlot(N,TLuf);
248 
249  timeplot->SetPlotLegend("LU solve");
250  timeplot->SetPlotColor(0.25,0.25,0);
251  timeplot->AddPlot(N,TLus);
252 
253  timeplot->SetPlotLegend("O(N**(3/2))");
254  timeplot->SetPlotColor(0.0,0.0,1.0);
255  timeplot->AddPlot(XN1,ON32);
256 
257  timeplot->SetPlotLegend("O(N)");
258  timeplot->SetPlotColor(1.0,0,0);
259  timeplot->AddPlot(XN,ON1);;
260 
261 
262  }
263  frame->Interact();
264 #endif
265  }
266 
267 }
268 
269 
Class collecting data for the description of piecewise linear geometries.
Definition: geometry.hxx:17
Sparse matrix class using CRS storage scheme.
void set_regionpoints(const std::initializer_list< std::initializer_list< double >> &il)
Set member via intializer list.
Definition: geometry.hxx:32
void set_regionvolumes(const std::initializer_list< double > &il)
Set member via intializer list.
Definition: geometry.hxx:38
Header for simple grid data class.
int main(void)
void solve(TArray< T > &Sol, const TArray< T > &Rhs)
Solve LU factorized system.
const TArray2< double > & get_points() const
Get array of point coordinates.
Definition: simplegrid.hxx:42
void assemble_heat_problem(const numcxx::SimpleGrid &Grid, const numcxx::DArray1 &BCfac, const numcxx::DArray1 &BCval, const numcxx::DArray1 &Source, const numcxx::DArray1 &Kappa, numcxx::DSparseMatrix &SGlobal, numcxx::DArray1 &Rhs)
Definition: fem2d.cxx:74
Main header of the library.
double h1norm(const numcxx::SimpleGrid &grid, const numcxx::DArray1 &u)
Definition: fem2d.cxx:302
void set_points(const std::initializer_list< std::initializer_list< double >> &il)
Set member via intializer list.
Definition: geometry.hxx:23
void set_regionnumbers(const std::initializer_list< int > &il)
Set member via intializer list.
Definition: geometry.hxx:35
const double Dirichlet
BC value marking Dirichlet boundary condition.
Definition: fem2d.hxx:11
void set_bfaces(const std::initializer_list< std::initializer_list< int >> &il)
Set member via intializer list.
Definition: geometry.hxx:26
void set_bfaceregions(const std::initializer_list< int > &il)
Set member via intializer list.
Definition: geometry.hxx:29
One dimensional array class.
Definition: tarray1.hxx:31
Class containing data for simple grid data structure.
Definition: simplegrid.hxx:19
Header for adapter beteween vtkfig dataset and simple grid.
double cpu_clock()
cpu time in seconds
Definition: util.ixx:101
Bridge class for using umfpack as solver for vmatrix.
void calc_hminmax(double &hmin, double &hmax) const
Calculate some grid data.
Definition: simplegrid.cxx:55
const int npoints() const
Return number of points.
Definition: simplegrid.hxx:66
void update()
Perform actual computation of LU factorization.
double l2norm(const numcxx::SimpleGrid &grid, const numcxx::DArray1 &u)
Definition: fem2d.cxx:275