Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
matrix_out.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_matrix_out_h
17 # define dealii_matrix_out_h
18 
19 # include <deal.II/base/config.h>
20 
21 # include <deal.II/base/data_out_base.h>
22 
23 # include <deal.II/lac/block_sparse_matrix.h>
24 # include <deal.II/lac/sparse_matrix.h>
25 
26 # ifdef DEAL_II_WITH_TRILINOS
27 # include <deal.II/lac/trilinos_block_sparse_matrix.h>
28 # include <deal.II/lac/trilinos_sparse_matrix.h>
29 # endif
30 
31 
32 DEAL_II_NAMESPACE_OPEN
33 
72 class MatrixOut : public DataOutInterface<2, 2>
73 {
74 public:
79 
84  struct Options
85  {
91 
101  unsigned int block_size;
102 
107 
112  Options(const bool show_absolute_values = false,
113  const unsigned int block_size = 1,
114  const bool discontinuous = false);
115  };
116 
120  virtual ~MatrixOut() override = default;
121 
139  template <class Matrix>
140  void
141  build_patches(const Matrix & matrix,
142  const std::string &name,
143  const Options options = Options(false, 1, false));
144 
145 private:
151 
157  std::vector<Patch> patches;
158 
162  std::string name;
163 
168  virtual const std::vector<Patch> &
169  get_patches() const override;
170 
175  virtual std::vector<std::string>
176  get_dataset_names() const override;
177 
185  template <class Matrix>
186  static double
187  get_gridpoint_value(const Matrix & matrix,
188  const size_type i,
189  const size_type j,
190  const Options & options);
191 };
192 
193 
194 /* ---------------------- Template and inline functions ------------- */
195 
196 
197 namespace internal
198 {
199  namespace MatrixOutImplementation
200  {
204  template <typename number>
205  double
206  get_element(const ::SparseMatrix<number> &matrix,
207  const types::global_dof_index i,
208  const types::global_dof_index j)
209  {
210  return matrix.el(i, j);
211  }
212 
213 
214 
218  template <typename number>
219  double
220  get_element(const ::BlockSparseMatrix<number> &matrix,
221  const types::global_dof_index i,
222  const types::global_dof_index j)
223  {
224  return matrix.el(i, j);
225  }
226 
227 
228 # ifdef DEAL_II_WITH_TRILINOS
229 
232  inline double
233  get_element(const TrilinosWrappers::SparseMatrix &matrix,
234  const types::global_dof_index i,
235  const types::global_dof_index j)
236  {
237  return matrix.el(i, j);
238  }
239 
240 
241 
246  inline double
247  get_element(const TrilinosWrappers::BlockSparseMatrix &matrix,
248  const types::global_dof_index i,
249  const types::global_dof_index j)
250  {
251  return matrix.el(i, j);
252  }
253 # endif
254 
255 
256 # ifdef DEAL_II_WITH_PETSC
257  // no need to do anything: PETSc matrix objects do not distinguish
258  // between operator() and el(i,j), so we can safely access elements
259  // through the generic function below
260 # endif
261 
262 
268  template <class Matrix>
269  double
270  get_element(const Matrix & matrix,
271  const types::global_dof_index i,
272  const types::global_dof_index j)
273  {
274  return matrix(i, j);
275  }
276  } // namespace MatrixOutImplementation
277 } // namespace internal
278 
279 
280 
281 template <class Matrix>
282 inline double
283 MatrixOut::get_gridpoint_value(const Matrix & matrix,
284  const size_type i,
285  const size_type j,
286  const Options & options)
287 {
288  // special case if block size is
289  // one since we then don't need all
290  // that loop overhead
291  if (options.block_size == 1)
292  {
293  if (options.show_absolute_values == true)
294  return std::fabs(
295  internal::MatrixOutImplementation::get_element(matrix, i, j));
296  else
297  return internal::MatrixOutImplementation::get_element(matrix, i, j);
298  }
299 
300  // if blocksize greater than one,
301  // then compute average of elements
302  double average = 0;
303  size_type n_elements = 0;
304  for (size_type row = i * options.block_size;
305  row <
306  std::min(size_type(matrix.m()), size_type((i + 1) * options.block_size));
307  ++row)
308  for (size_type col = j * options.block_size;
309  col < std::min(size_type(matrix.m()),
310  size_type((j + 1) * options.block_size));
311  ++col, ++n_elements)
312  if (options.show_absolute_values == true)
313  average += std::fabs(
314  internal::MatrixOutImplementation::get_element(matrix, row, col));
315  else
316  average +=
317  internal::MatrixOutImplementation::get_element(matrix, row, col);
318  average /= n_elements;
319  return average;
320 }
321 
322 
323 
324 template <class Matrix>
325 void
326 MatrixOut::build_patches(const Matrix & matrix,
327  const std::string &name,
328  const Options options)
329 {
330  size_type gridpoints_x = (matrix.n() / options.block_size +
331  (matrix.n() % options.block_size != 0 ? 1 : 0)),
332  gridpoints_y = (matrix.m() / options.block_size +
333  (matrix.m() % options.block_size != 0 ? 1 : 0));
334 
335  // If continuous, the number of
336  // plotted patches is matrix size-1
337  if (!options.discontinuous)
338  {
339  --gridpoints_x;
340  --gridpoints_y;
341  }
342 
343  // first clear old data and set it
344  // to virgin state
345  patches.clear();
346  patches.resize((gridpoints_x) * (gridpoints_y));
347 
348  // now build the patches
349  size_type index = 0;
350  for (size_type i = 0; i < gridpoints_y; ++i)
351  for (size_type j = 0; j < gridpoints_x; ++j, ++index)
352  {
353  // within each patch, order the points in such a way that if some
354  // graphical output program (such as gnuplot) plots the quadrilaterals
355  // as two triangles, then the diagonal of the quadrilateral which cuts
356  // it into the two printed triangles is parallel to the diagonal of the
357  // matrix, rather than perpendicular to it. this has the advantage that,
358  // for example, the unit matrix is plotted as a straight rim, rather
359  // than as a series of bumps and valleys along the diagonal
360  patches[index].vertices[0](0) = j;
361  patches[index].vertices[0](1) = -static_cast<signed int>(i);
362  patches[index].vertices[1](0) = j;
363  patches[index].vertices[1](1) = -static_cast<signed int>(i + 1);
364  patches[index].vertices[2](0) = j + 1;
365  patches[index].vertices[2](1) = -static_cast<signed int>(i);
366  patches[index].vertices[3](0) = j + 1;
367  patches[index].vertices[3](1) = -static_cast<signed int>(i + 1);
368  // next scale all the patch
369  // coordinates by the block
370  // size, to get original
371  // coordinates
372  for (auto &vertex : patches[index].vertices)
373  vertex *= options.block_size;
374 
375  patches[index].n_subdivisions = 1;
376 
377  patches[index].data.reinit(1, 4);
378  if (options.discontinuous)
379  {
380  patches[index].data(0, 0) =
381  get_gridpoint_value(matrix, i, j, options);
382  patches[index].data(0, 1) =
383  get_gridpoint_value(matrix, i, j, options);
384  patches[index].data(0, 2) =
385  get_gridpoint_value(matrix, i, j, options);
386  patches[index].data(0, 3) =
387  get_gridpoint_value(matrix, i, j, options);
388  }
389  else
390  {
391  patches[index].data(0, 0) =
392  get_gridpoint_value(matrix, i, j, options);
393  patches[index].data(0, 1) =
394  get_gridpoint_value(matrix, i + 1, j, options);
395  patches[index].data(0, 2) =
396  get_gridpoint_value(matrix, i, j + 1, options);
397  patches[index].data(0, 3) =
398  get_gridpoint_value(matrix, i + 1, j + 1, options);
399  }
400  };
401 
402  // finally set the name
403  this->name = name;
404 }
405 
406 
407 
408 /*---------------------------- matrix_out.h ---------------------------*/
409 
410 DEAL_II_NAMESPACE_CLOSE
411 
412 #endif
413 /*---------------------------- matrix_out.h ---------------------------*/
types::global_dof_index size_type
Definition: matrix_out.h:78
static double get_gridpoint_value(const Matrix &matrix, const size_type i, const size_type j, const Options &options)
Definition: matrix_out.h:283
Contents is actually a matrix.
virtual std::vector< std::string > get_dataset_names() const override
Definition: matrix_out.cc:40
std::vector< Patch > patches
Definition: matrix_out.h:157
virtual ~MatrixOut() override=default
unsigned int block_size
Definition: matrix_out.h:101
std::string name
Definition: matrix_out.h:162
void build_patches(const Matrix &matrix, const std::string &name, const Options options=Options(false, 1, false))
Definition: matrix_out.h:326
bool show_absolute_values
Definition: matrix_out.h:90
unsigned int global_dof_index
Definition: types.h:89
Options(const bool show_absolute_values=false, const unsigned int block_size=1, const bool discontinuous=false)
Definition: matrix_out.cc:21
virtual const std::vector< Patch > & get_patches() const override
Definition: matrix_out.cc:32