InstallingDealii: diffusion-equation.cc

File diffusion-equation.cc, 8.3 KB (added by leemo@…, 14 years ago)

diffusion equation solver -- backup

Line 
1#include <grid/tria.h>
2#include <dofs/dof_handler.h>
3#include <grid/grid_generator.h>
4
5#include <grid/tria_accessor.h>
6#include <grid/tria_iterator.h>
7#include <dofs/dof_accessor.h>
8
9#include <fe/fe_q.h>
10
11#include <dofs/dof_tools.h>
12
13#include <fe/fe_values.h>
14#include <base/quadrature_lib.h>
15
16#include <base/function.h>
17#include <numerics/vectors.h>
18#include <numerics/matrices.h>
19
20#include <lac/vector.h>
21#include <lac/full_matrix.h>
22#include <lac/sparse_matrix.h>
23#include <lac/solver_cg.h>
24#include <lac/vector_memory.h>
25#include <lac/precondition.h>
26
27#include <numerics/data_out.h>
28#include <fstream>
29#include <iostream>
30
31#include <sstream>
32
33
34class DiffusionProblem 
35{
36  public:
37    DiffusionProblem ();
38
39    void run ();
40   
41  private:
42    void make_grid_and_dofs ();
43    void set_initial_condition ();
44    void assemble_system ();
45    void solve ();
46    void output_results (int time_step_counter) const;
47
48    Triangulation<2>     triangulation;
49    FE_Q<2>              fe;
50    DoFHandler<2>        dof_handler;
51
52    SparsityPattern      sparsity_pattern;
53    SparseMatrix<double> system_matrix;
54
55    Vector<double>       solution;
56    Vector<double>       system_rhs;
57   
58    const double         dt;
59   
60    bool                 matrix_is_assembled;
61   
62    Tensor<2,2>          conductivity_tensor;          // first two=rank (ie matrix), second two=SPACE_DIM
63    double               surface_area_to_volume_ratio;
64    double               capacitance;
65};
66
67
68DiffusionProblem::DiffusionProblem () :
69                fe (1),
70                dof_handler (triangulation),
71                dt (0.01),
72                matrix_is_assembled(false)
73{
74    surface_area_to_volume_ratio = 1400; // 1/cm
75    capacitance = 1.0; // uF/cm^2
76   
77    conductivity_tensor[0][0] = 1.75;
78    conductivity_tensor[1][0] = 0.0;
79    conductivity_tensor[0][1] = 0.0;
80    conductivity_tensor[1][1] = 1.75;   
81}
82
83void DiffusionProblem::set_initial_condition ()
84{
85//  Triangulation<2>::active_cell_iterator cell, endc;
86//  cell = triangulation.begin_active();
87//  endc = triangulation.end();
88 
89    DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(),
90                                      endc = dof_handler.end();
91 
92  const unsigned int   dofs_per_cell = fe.dofs_per_cell;
93  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
94 
95  for (; cell!=endc; ++cell)
96  {
97      cell->get_dof_indices (local_dof_indices);
98      for (unsigned int vertex=0;
99           vertex < GeometryInfo<2>::vertices_per_cell;
100           ++vertex)
101        {
102          const Point<2> position
103            = (cell->vertex(vertex));
104           
105
106         
107          if (position(0)==0 && position(1)==0)
108            {
109               std::cout << "Found the origin" <<std::endl;
110               solution(local_dof_indices[vertex]) = 1.0;
111            }
112          else
113          {
114            solution(local_dof_indices[vertex]) = 0.0;
115          }
116        };
117  }
118}
119
120
121void DiffusionProblem::make_grid_and_dofs ()
122{
123  GridGenerator::hyper_cube (triangulation, -1, 1);
124  triangulation.refine_global (5);
125  std::cout << "Number of active cells: "
126        << triangulation.n_active_cells()
127        << std::endl;
128  std::cout << "Total number of cells: "
129            << triangulation.n_cells()
130            << std::endl;
131 
132  dof_handler.distribute_dofs (fe);
133
134  std::cout << "Number of degrees of freedom: "
135            << dof_handler.n_dofs()
136            << std::endl;
137
138  sparsity_pattern.reinit (dof_handler.n_dofs(),
139                           dof_handler.n_dofs(),
140                           dof_handler.max_couplings_between_dofs());
141  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
142  sparsity_pattern.compress();
143
144  system_matrix.reinit (sparsity_pattern);
145
146  solution.reinit (dof_handler.n_dofs());
147  system_rhs.reinit (dof_handler.n_dofs());
148}
149
150
151void DiffusionProblem::assemble_system () 
152{
153  QGauss<2>  quadrature_formula(2);
154  FEValues<2> fe_values (fe, quadrature_formula, 
155                         UpdateFlags(update_values    |
156                                     update_gradients |
157                                     update_JxW_values));
158
159  const unsigned int   dofs_per_cell = fe.dofs_per_cell;
160  const unsigned int   n_q_points    = quadrature_formula.n_quadrature_points;
161
162  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
163  Vector<double>       cell_rhs (dofs_per_cell);
164
165  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
166
167  DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(),
168                                      endc = dof_handler.end();
169                                     
170  double dudt_coeff_term = 1; //capacitance*surface_area_to_volume_ratio;
171                                       
172  for (; cell!=endc; ++cell)
173  {
174      fe_values.reinit (cell);
175
176      cell_matrix = 0;
177      cell_rhs = 0;
178     
179      cell->get_dof_indices (local_dof_indices);
180
181      if(!matrix_is_assembled)
182      {
183        for (unsigned int i=0; i<dofs_per_cell; ++i)
184        {
185          for (unsigned int j=0; j<dofs_per_cell; ++j)
186          {
187            for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
188            {
189              cell_matrix(i,j) += (fe_values.shape_grad (i, q_point) *
190                                   (conductivity_tensor *
191                                   fe_values.shape_grad (j, q_point)) *
192                                   fe_values.JxW (q_point));
193                                 
194              cell_matrix(i,j) += fe_values.shape_value(i, q_point) *
195                                  fe_values.shape_value(j, q_point) *
196                                  fe_values.JxW (q_point) *
197                                  dudt_coeff_term / dt;
198            }
199          }
200        }
201      }
202     
203      for (unsigned int i=0; i<dofs_per_cell; ++i)
204      {
205        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
206        {
207          cell_rhs(i) += (fe_values.shape_value (i, q_point) *
208                          1 *
209                          fe_values.JxW (q_point));
210                         
211          cell_rhs(i) += solution(local_dof_indices[i]) *
212                         fe_values.shape_value(i, q_point) *
213                         fe_values.JxW (q_point) *
214                         dudt_coeff_term / dt;                 
215        }
216      }
217
218      if(!matrix_is_assembled) 
219      {
220        for (unsigned int i=0; i<dofs_per_cell; ++i)
221        {
222          for (unsigned int j=0; j<dofs_per_cell; ++j)
223          {
224            system_matrix.add (local_dof_indices[i],
225                               local_dof_indices[j], 
226                               cell_matrix(i,j));
227          }
228        }
229      }
230
231      for (unsigned int i=0; i<dofs_per_cell; ++i)
232      {
233        system_rhs(local_dof_indices[i]) += cell_rhs(i);
234      }
235   };
236   
237   matrix_is_assembled = true;
238
239
240  std::map<unsigned int,double> boundary_values;
241  VectorTools::interpolate_boundary_values (dof_handler,
242                                            0,
243                                            ZeroFunction<2>(),
244                                            boundary_values);
245  MatrixTools::apply_boundary_values (boundary_values,
246                                      system_matrix,
247                                      solution,
248                                      system_rhs);
249}
250
251
252void DiffusionProblem::solve () 
253{
254  SolverControl           solver_control (1000, 1e-12);
255  PrimitiveVectorMemory<> vector_memory;
256  SolverCG<>              cg (solver_control, vector_memory);
257
258  cg.solve (system_matrix, solution, system_rhs,
259            PreconditionIdentity());
260}
261
262
263void DiffusionProblem::output_results (int time_step_counter) const
264{
265  DataOut<2> data_out;
266  data_out.attach_dof_handler (dof_handler);
267  std::stringstream filename;
268  filename << "solution_" << time_step_counter << ".gpl";
269  data_out.add_data_vector (solution, "solution");
270  data_out.build_patches ();
271
272  std::ofstream output (filename.str().c_str());
273  data_out.write_gnuplot (output);
274}
275
276
277void DiffusionProblem::run () 
278{
279  make_grid_and_dofs ();
280  set_initial_condition ();
281  double t=0;
282  int time_step_counter=0;
283  while (t<0.1)
284  {
285    assemble_system ();
286    solve ();
287    t += dt;
288    time_step_counter++;
289    std::cout << "Time = " << t << std::endl;
290    output_results (time_step_counter);
291  }
292}
293
294   
295
296int main () 
297{
298  DiffusionProblem diffusion_problem;
299  diffusion_problem.run ();
300  return 0;
301}