NUMCXX  0.13.20181108
Numerical library for small projects and teaching purposes
46-nonlin-fvm.cxx
Go to the documentation of this file.
1 
7 #include <cstdio>
8 #include <iostream>
9 #include <ctime>
10 #include <cmath>
11 #include <numcxx/numcxx.hxx>
12 #include <numcxx/simplegrid.hxx>
13 #include <numcxx/fvm2d.hxx>
14 #ifdef VTKFIG
16 #endif
17 
18 #include "vtkfigFrame.h"
19 #include "vtkfigDataSet.h"
20 #include "vtkfigGridView.h"
21 #include "vtkfigScalarView.h"
22 
23 
24 
25 int main(void)
26 {
27  numcxx::Geometry Geometry;
28  Geometry.set_points({
29  {-2,0},
30  {0,0},
31  {2,0},
32  {2,2},
33  {0.5,2},
34  {0,1},
35  {-0.5,2},
36  {-2,2}
37  });
38 
39  Geometry.set_bfaces({
40  {0,1},
41  {1,2},
42  {2,3},
43  {3,4},
44  {4,5},
45  {5,6},
46  {6,7},
47  {7,0},
48  {5,1},
49  });
50 
51 
52  Geometry.set_bfaceregions({1,1,2,3,4,4,3,2,5});
53 
54  Geometry.set_regionpoints({
55  {-0.5,1},
56  {0.5,1}
57  });
58  Geometry.set_regionnumbers({1,2});
59  double vol=0.001;
60  Geometry.set_regionvolumes({vol,vol});
61 
62 
63 
64  numcxx::SimpleGrid grid(Geometry,"zpaAqDV");
65 
66 
67 
68  numcxx::DArray1 bcfac(8);
69  numcxx::DArray1 bcval(8);
70  bcfac=0;
71  bcval=0;
72  bcfac(4)=fvm2d::Dirichlet;
73  bcfac(1)=fvm2d::Dirichlet;
74  bcval(4)=1.0;
75  bcval(1)=0.0;
76 
77  auto nnodes=grid.npoints();
78 
79  numcxx::DArray1 source(nnodes);
80  source=0;
81 
82  auto fkappa = [](double u, double &kappa,double &dkappa )
83  {
84  kappa=0.001+1000.0*u*u*u*u*u*u;
85  dkappa=6000.0*u*u*u*u*u;
86  };
87 
88 
89  numcxx::DSparseMatrix SGlobal(nnodes,nnodes);
90  numcxx::DSolverUMFPACK Solver(SGlobal);
91 
92  numcxx::DArray1 Rhs(nnodes);
93  numcxx::DArray1 Sol(nnodes);
94  numcxx::DArray1 Res(nnodes);
95  numcxx::DArray1 Upd(nnodes);
96 
97  Sol=0.0;
98  fvm2d::initialize_bc(grid,bcval,Sol);
99 
100 
101 
102  int iter=0;
103  double norm=1.0;
104  double oldnorm=1.0;
105  double d=0.1;
106  double ddelta=1.2;
107 
108 
109  while (iter<100 && norm >1.0e-13)
110  {
111  fvm2d::assemble_and_apply_nonlinear_heat(grid,bcfac,bcval, source,fkappa, SGlobal, Sol, Res);
112  Solver.update(SGlobal);
113  Solver.solve(Upd,Res);
114  oldnorm=norm;
115  norm=numcxx::norm2(Upd);
116  printf("iter=%d norm=%8.4e contract=%8.5e\n", iter,norm,norm/oldnorm);
117  Sol=Sol-d*Upd;
118  iter++;
119  d= std::min(1.0, d*ddelta);
120  }
121 
122 
123 
124 #ifdef VTKFIG
125  auto griddata=numcxx::vtkfigDataSet(grid);
126  griddata->SetPointScalar(Sol ,"Sol");
127 
128  auto frame=vtkfig::Frame::New();
129  frame->SetSize(800,400);
130  frame->SetLayout(2,1);
131 
132  auto gridview=vtkfig::GridView::New();
133  gridview->SetData(griddata);
134  frame->AddFigure(gridview,0);
135 
136  auto solview=vtkfig::ScalarView::New();
137  solview->SetData(griddata,"Sol");
138  frame->AddFigure(solview,1);
139 
140 
141 
142  frame->Interact();
143 #endif
144 }
145 
146 
void initialize_bc(numcxx::SimpleGrid &grid, numcxx::DArray1 &g, numcxx::DArray1 &Sol)
Definition: fvm2d.cxx:291
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.
A::value_type norm2(const A &a)
Euklidean norm of array or expression.
Definition: util.ixx:47
void solve(TArray< T > &Sol, const TArray< T > &Rhs)
Solve LU factorized system.
Main header of the library.
A::value_type min(const A &a)
Minimum of of array or expression.
Definition: util.ixx:59
void assemble_and_apply_nonlinear_heat(const numcxx::SimpleGrid &grid, const numcxx::DArray1 &bcfac, const numcxx::DArray1 &bcval, const numcxx::DArray1 &source, std::function< void(const double, double &, double &)> fkappa, numcxx::DSparseMatrix &SGlobal, numcxx::DArray1 &Sol, numcxx::DArray1 &Rhs)
Definition: fvm2d.cxx:219
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
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
int main(void)
Class containing data for simple grid data structure.
Definition: simplegrid.hxx:19
Header for adapter beteween vtkfig dataset and simple grid.
Bridge class for using umfpack as solver for vmatrix.
const int npoints() const
Return number of points.
Definition: simplegrid.hxx:66
A::value_type norm(const A &a)
Euklidean norm of array or expression.
Definition: util.hxx:54
const double Dirichlet
BC value marking Dirichlet boundary condition.
Definition: fvm2d.hxx:13
void update()
Perform actual computation of LU factorization.