NUMCXX  0.13.20181108
Numerical library for small projects and teaching purposes
44-transient-heat-fe.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/fem2d.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  // Parameters to test
28  double h=0.05; // Space step size
29  double tau=0.01; // Time step size
30  double theta=1.0; // Implicit Euler vs CN vs explicit Euler
31  double T=100.0; // Length of time interval
32  bool lump=false; // Mass lumping
33 
34 
35 
36 
38  // Describe geometry
39  numcxx::Geometry Geometry;
40  Geometry.set_points({
41  {0,0},
42  {1,0},
43  {1,1},
44  {0,1},
45  {0.7,0.7},
46  });
47 
48  Geometry.set_bfaces({
49  {0,1},
50  {1,2},
51  {2,3},
52  {3,0}
53  });
54 
55 
56  Geometry.set_bfaceregions({1,2,3,4});
57 
58  Geometry.set_regionpoints({
59  {0.5,0.5}
60  });
61 
62  Geometry.set_regionnumbers({1});
63 
64  double vol= h*h*0.25;
65  Geometry.set_regionvolumes({vol});
66 
67 
69  // Create grid
70  numcxx::SimpleGrid grid(Geometry,"zpaAqDV");
71 
72 
74  // Prepare visualization
75 
76  auto griddata=numcxx::vtkfigDataSet(grid);
77 
78  auto frame=vtkfig::Frame::New();
79  frame->SetSize(1200,600);
80  frame->SetLayout(2,1);
81 
82  auto gridview=vtkfig::GridView::New();
83  gridview->SetData(griddata);
84  frame->AddFigure(gridview,0);
85 
86  auto solview=vtkfig::ScalarView::New();
87  solview->SetData(griddata,"Sol");
88  solview->SetValueRange(0,1);
89  frame->AddFigure(solview,1);
90 
91 
93  // Set up problem data
94  numcxx::DArray1 bcfac(6);
95  numcxx::DArray1 bcval(6);
96  bcfac=0;
97  bcval=0;
98 
99  bcfac(3)=fem2d::Dirichlet;
100  bcfac(1)=fem2d::Dirichlet;
101  bcval(3)=1.0;
102  bcval(1)=0.0;
103 
104  auto nnodes=grid.npoints();
105 
106  numcxx::DArray1 source(nnodes);
107  numcxx::DArray1 kappa(nnodes);
108  kappa=1.0e-2;
109  source=0;
110 
111  int N=T/tau;
112 
113  numcxx::DSparseMatrix SGlobal(nnodes,nnodes);
114  numcxx::DArray1 Rhs(nnodes);
115  numcxx::DArray1 Sol(nnodes);
116  numcxx::DArray1 OldSol(nnodes);
117  Sol=0.0;
118 
119  numcxx::DSolverUMFPACK Solver(SGlobal);
120 
121  for (int n=0;n<N;n++)
122  {
123  OldSol=Sol;
124  fem2d::assemble_transient_heat_problem(grid,bcfac,bcval,source,kappa,tau,theta, lump, OldSol, SGlobal, Rhs);
125 
126  if (theta>0.0 || !lump)
127  {
128  Solver.update(SGlobal);
129  Solver.solve(Sol,Rhs);
130  }
131  else
132  {
133  for (int i=0;i<Sol.size();i++)
134  Sol(i)=Rhs(i)/SGlobal(i,i);
135  }
136 
137 #ifdef VTKFIG
138 
139  frame->SetFrameTitle("Time="+std::to_string(tau*n));
140  griddata->SetPointScalar(Sol ,"Sol");
141  frame->Show();
142 
143 #endif
144  }
145  frame->Interact();
146 
147 }
148 
149 
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.
size_t size() const
Obtain size of array.
Definition: tarray.hxx:58
void solve(TArray< T > &Sol, const TArray< T > &Rhs)
Solve LU factorized system.
Main header of the library.
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
int main(void)
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 assemble_transient_heat_problem(const numcxx::SimpleGrid &grid, const numcxx::DArray1 &bcfac, const numcxx::DArray1 &bcval, const numcxx::DArray1 &source, const numcxx::DArray1 &kappa, double tau, double theta, bool lump, numcxx::DArray1 &OldSol, numcxx::DSparseMatrix &SGlobal, numcxx::DArray1 &Rhs)
Definition: fem2d.cxx:165
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.
Bridge class for using umfpack as solver for vmatrix.
const int npoints() const
Return number of points.
Definition: simplegrid.hxx:66
void update()
Perform actual computation of LU factorization.