Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
matrix_out.h
1 // ---------------------------------------------------------------------
2 // @f$Id: matrix_out.h 31932 2013-12-08 02:15:54Z heister @f$
3 //
4 // Copyright (C) 2001 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__matrix_out_h
18 #define __deal2__matrix_out_h
19 
20 #include <deal.II/base/config.h>
21 #include <deal.II/base/data_out_base.h>
22 #include <deal.II/lac/sparse_matrix.h>
23 #include <deal.II/lac/block_sparse_matrix.h>
24 
25 DEAL_II_NAMESPACE_OPEN
26 
93 class MatrixOut : public DataOutInterface<2,2>
94 {
95 public:
100 
107  struct Options
108  {
119 
137  unsigned int block_size;
138 
145 
152  Options (const bool show_absolute_values = false,
153  const unsigned int block_size = 1,
154  const bool discontinuous = false);
155  };
156 
161  virtual ~MatrixOut ();
162 
200  template <class Matrix>
201  void build_patches (const Matrix &matrix,
202  const std::string &name,
203  const Options options = Options(false, 1, false));
204 
205 private:
206 
222  typedef ::DataOutBase::Patch<2,2> Patch;
223 
232  std::vector<Patch> patches;
233 
238  std::string name;
239 
246  virtual const std::vector<Patch> &
247  get_patches () const;
248 
255  virtual std::vector<std::string> get_dataset_names () const;
256 
261  template <typename number>
262  static double get_element (const SparseMatrix<number> &matrix,
263  const size_type i,
264  const size_type j);
265 
271  template <typename number>
272  static double get_element (const BlockSparseMatrix<number> &matrix,
273  const size_type i,
274  const size_type j);
275 
284  template <class Matrix>
285  static double get_element (const Matrix &matrix,
286  const size_type i,
287  const size_type j);
288 
302  template <class Matrix>
303  static double get_gridpoint_value (const Matrix &matrix,
304  const size_type i,
305  const size_type j,
306  const Options &options);
307 };
308 
309 
310 /* ---------------------- Template and inline functions ------------- */
311 
312 
313 template <typename number>
314 inline
315 double
317  const size_type i,
318  const size_type j)
319 {
320  return matrix.el(i,j);
321 }
322 
323 
324 
325 
326 template <typename number>
327 inline
328 double
330  const size_type i,
331  const size_type j)
332 {
333  return matrix.el(i,j);
334 }
335 
336 
337 
338 
339 template <class Matrix>
340 inline
341 double
342 MatrixOut::get_element (const Matrix &matrix,
343  const size_type i,
344  const size_type j)
345 {
346  return matrix(i,j);
347 }
348 
349 
350 
351 template <class Matrix>
352 inline
353 double
354 MatrixOut::get_gridpoint_value (const Matrix &matrix,
355  const size_type i,
356  const size_type j,
357  const Options &options)
358 {
359  // special case if block size is
360  // one since we then don't need all
361  // that loop overhead
362  if (options.block_size == 1)
363  {
364  if (options.show_absolute_values == true)
365  return std::fabs(get_element (matrix, i, j));
366  else
367  return get_element (matrix, i, j);
368  }
369 
370  // if blocksize greater than one,
371  // then compute average of elements
372  double average = 0;
373  size_type n_elements = 0;
374  for (size_type row=i*options.block_size;
375  row < std::min(size_type(matrix.m()),
376  size_type((i+1)*options.block_size)); ++row)
377  for (size_type col=j*options.block_size;
378  col < std::min(size_type(matrix.m()),
379  size_type((j+1)*options.block_size)); ++col, ++n_elements)
380  if (options.show_absolute_values == true)
381  average += std::fabs(get_element (matrix, row, col));
382  else
383  average += get_element (matrix, row, col);
384  average /= n_elements;
385  return average;
386 }
387 
388 
389 
390 template <class Matrix>
391 void
392 MatrixOut::build_patches (const Matrix &matrix,
393  const std::string &name,
394  const Options options)
395 {
396  size_type
397  gridpoints_x = (matrix.n() / options.block_size
398  +
399  (matrix.n() % options.block_size != 0 ? 1 : 0)),
400  gridpoints_y = (matrix.m() / options.block_size
401  +
402  (matrix.m() % options.block_size != 0 ? 1 : 0));
403 
404  // If continuous, the number of
405  // plotted patches is matrix size-1
406  if (!options.discontinuous)
407  {
408  --gridpoints_x;
409  --gridpoints_y;
410  }
411 
412  // first clear old data and set it
413  // to virgin state
414  patches.clear ();
415  patches.resize ((gridpoints_x) * (gridpoints_y));
416 
417  // now build the patches
418  size_type index=0;
419  for (size_type i=0; i<gridpoints_y; ++i)
420  for (size_type j=0; j<gridpoints_x; ++j, ++index)
421  {
422  // within each patch, order
423  // the points in such a way
424  // that if some graphical
425  // output program (such as
426  // gnuplot) plots the
427  // quadrilaterals as two
428  // triangles, then the
429  // diagonal of the
430  // quadrilateral which cuts
431  // it into the two printed
432  // triangles is parallel to
433  // the diagonal of the
434  // matrix, rather than
435  // perpendicular to it. this
436  // has the advantage that,
437  // for example, the unit
438  // matrix is plotted as a
439  // straight rim, rather than
440  // as a series of bumps and
441  // valleys along the diagonal
442  patches[index].vertices[0](0) = j;
443  patches[index].vertices[0](1) = static_cast<signed int>(-i);
444  patches[index].vertices[1](0) = j;
445  patches[index].vertices[1](1) = static_cast<signed int>(-i-1);
446  patches[index].vertices[2](0) = j+1;
447  patches[index].vertices[2](1) = static_cast<signed int>(-i);
448  patches[index].vertices[3](0) = j+1;
449  patches[index].vertices[3](1) = static_cast<signed int>(-i-1);
450  // next scale all the patch
451  // coordinates by the block
452  // size, to get original
453  // coordinates
454  for (unsigned int v=0; v<4; ++v)
455  patches[index].vertices[v] *= options.block_size;
456 
457  patches[index].n_subdivisions = 1;
458 
459  patches[index].data.reinit (1,4);
460  if (options.discontinuous)
461  {
462  patches[index].data(0,0) = get_gridpoint_value(matrix, i, j, options);
463  patches[index].data(0,1) = get_gridpoint_value(matrix, i, j, options);
464  patches[index].data(0,2) = get_gridpoint_value(matrix, i, j, options);
465  patches[index].data(0,3) = get_gridpoint_value(matrix, i, j, options);
466  }
467  else
468  {
469  patches[index].data(0,0) = get_gridpoint_value(matrix, i, j, options);
470  patches[index].data(0,1) = get_gridpoint_value(matrix, i+1, j, options);
471  patches[index].data(0,2) = get_gridpoint_value(matrix, i, j+1, options);
472  patches[index].data(0,3) = get_gridpoint_value(matrix, i+1, j+1, options);
473  }
474  };
475 
476  // finally set the name
477  this->name = name;
478 }
479 
480 
481 
482 /*---------------------------- matrix_out.h ---------------------------*/
483 
484 DEAL_II_NAMESPACE_CLOSE
485 
486 #endif
487 /*---------------------------- matrix_out.h ---------------------------*/
virtual const std::vector< Patch > & get_patches() const
virtual ~MatrixOut()
static double get_element(const SparseMatrix< number > &matrix, const size_type i, const size_type j)
Definition: matrix_out.h:316
static double get_gridpoint_value(const Matrix &matrix, const size_type i, const size_type j, const Options &options)
Definition: matrix_out.h:354
std::vector< Patch > patches
Definition: matrix_out.h:232
unsigned int block_size
Definition: matrix_out.h:137
number el(const size_type i, const size_type j) const
unsigned int global_dof_index
Definition: types.h:100
::DataOutBase::Patch< 2, 2 > Patch
Definition: matrix_out.h:222
std::string name
Definition: matrix_out.h:238
void build_patches(const Matrix &matrix, const std::string &name, const Options options=Options(false, 1, false))
Definition: matrix_out.h:392
types::global_dof_index size_type
Definition: matrix_out.h:99
value_type el(const size_type i, const size_type j) const
virtual std::vector< std::string > get_dataset_names() const
Options(const bool show_absolute_values=false, const unsigned int block_size=1, const bool discontinuous=false)