Reference documentation for deal.II version GIT b6bf1e606d 2022-08-11 15:25:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
matrix_out.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2021 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 
22 
25 
26 # ifdef DEAL_II_WITH_TRILINOS
29 # endif
30 
31 
33 
71 class MatrixOut : public DataOutInterface<2, 2>
72 {
73 public:
78 
83  struct Options
84  {
90 
100  unsigned int block_size;
101 
106 
111  Options(const bool show_absolute_values = false,
112  const unsigned int block_size = 1,
113  const bool discontinuous = false);
114  };
115 
119  virtual ~MatrixOut() override = default;
120 
138  template <class Matrix>
139  void
140  build_patches(const Matrix & matrix,
141  const std::string &name,
142  const Options options = Options(false, 1, false));
143 
144 private:
150 
156  std::vector<Patch> patches;
157 
161  std::string name;
162 
167  virtual const std::vector<Patch> &
168  get_patches() const override;
169 
174  virtual std::vector<std::string>
175  get_dataset_names() const override;
176 
184  template <class Matrix>
185  static double
186  get_gridpoint_value(const Matrix & matrix,
187  const size_type i,
188  const size_type j,
189  const Options & options);
190 };
191 
192 
193 /* ---------------------- Template and inline functions ------------- */
194 
195 
196 namespace internal
197 {
198  namespace MatrixOutImplementation
199  {
203  template <typename number>
204  double
205  get_element(const ::SparseMatrix<number> &matrix,
206  const types::global_dof_index i,
207  const types::global_dof_index j)
208  {
209  return matrix.el(i, j);
210  }
211 
212 
213 
217  template <typename number>
218  double
219  get_element(const ::BlockSparseMatrix<number> &matrix,
220  const types::global_dof_index i,
221  const types::global_dof_index j)
222  {
223  return matrix.el(i, j);
224  }
225 
226 
227 # ifdef DEAL_II_WITH_TRILINOS
231  inline double
233  const types::global_dof_index i,
234  const types::global_dof_index j)
235  {
236  return matrix.el(i, j);
237  }
238 
239 
240 
245  inline double
247  const types::global_dof_index i,
248  const types::global_dof_index j)
249  {
250  return matrix.el(i, j);
251  }
252 # endif
253 
254 
255 # ifdef DEAL_II_WITH_PETSC
256  // no need to do anything: PETSc matrix objects do not distinguish
257  // between operator() and el(i,j), so we can safely access elements
258  // through the generic function below
259 # endif
260 
261 
267  template <class Matrix>
268  double
269  get_element(const Matrix & matrix,
270  const types::global_dof_index i,
271  const types::global_dof_index j)
272  {
273  return matrix(i, j);
274  }
275  } // namespace MatrixOutImplementation
276 } // namespace internal
277 
278 
279 
280 template <class Matrix>
281 inline double
283  const size_type i,
284  const size_type j,
285  const Options & options)
286 {
287  // special case if block size is
288  // one since we then don't need all
289  // that loop overhead
290  if (options.block_size == 1)
291  {
292  if (options.show_absolute_values == true)
293  return std::fabs(
295  else
297  }
298 
299  // if blocksize greater than one,
300  // then compute average of elements
301  double average = 0;
302  size_type n_elements = 0;
303  for (size_type row = i * options.block_size;
304  row <
305  std::min(size_type(matrix.m()), size_type((i + 1) * options.block_size));
306  ++row)
307  for (size_type col = j * options.block_size;
308  col < std::min(size_type(matrix.m()),
309  size_type((j + 1) * options.block_size));
310  ++col, ++n_elements)
311  if (options.show_absolute_values == true)
312  average += std::fabs(
314  else
315  average +=
317  average /= n_elements;
318  return average;
319 }
320 
321 
322 
323 template <class Matrix>
324 void
326  const std::string &name,
327  const Options options)
328 {
329  size_type gridpoints_x = (matrix.n() / options.block_size +
330  (matrix.n() % options.block_size != 0 ? 1 : 0)),
331  gridpoints_y = (matrix.m() / options.block_size +
332  (matrix.m() % options.block_size != 0 ? 1 : 0));
333 
334  // If continuous, the number of
335  // plotted patches is matrix size-1
336  if (!options.discontinuous)
337  {
338  --gridpoints_x;
339  --gridpoints_y;
340  }
341 
342  // first clear old data and set it
343  // to virgin state
344  patches.clear();
345  patches.resize((gridpoints_x) * (gridpoints_y));
346 
347  // now build the patches
348  size_type index = 0;
349  for (size_type i = 0; i < gridpoints_y; ++i)
350  for (size_type j = 0; j < gridpoints_x; ++j, ++index)
351  {
352  patches[index].n_subdivisions = 1;
353  patches[index].reference_cell = ReferenceCells::Quadrilateral;
354 
355  // within each patch, order the points in such a way that if some
356  // graphical output program (such as gnuplot) plots the quadrilaterals
357  // as two triangles, then the diagonal of the quadrilateral which cuts
358  // it into the two printed triangles is parallel to the diagonal of the
359  // matrix, rather than perpendicular to it. this has the advantage that,
360  // for example, the unit matrix is plotted as a straight ridge, rather
361  // than as a series of bumps and valleys along the diagonal
362  patches[index].vertices[0](0) = j;
363  patches[index].vertices[0](1) = -static_cast<signed int>(i);
364  patches[index].vertices[1](0) = j;
365  patches[index].vertices[1](1) = -static_cast<signed int>(i + 1);
366  patches[index].vertices[2](0) = j + 1;
367  patches[index].vertices[2](1) = -static_cast<signed int>(i);
368  patches[index].vertices[3](0) = j + 1;
369  patches[index].vertices[3](1) = -static_cast<signed int>(i + 1);
370  // next scale all the patch
371  // coordinates by the block
372  // size, to get original
373  // coordinates
374  for (auto &vertex : patches[index].vertices)
375  vertex *= options.block_size;
376 
377  patches[index].n_subdivisions = 1;
378 
379  patches[index].data.reinit(1, 4);
380  if (options.discontinuous)
381  {
382  patches[index].data(0, 0) =
383  get_gridpoint_value(matrix, i, j, options);
384  patches[index].data(0, 1) =
385  get_gridpoint_value(matrix, i, j, options);
386  patches[index].data(0, 2) =
387  get_gridpoint_value(matrix, i, j, options);
388  patches[index].data(0, 3) =
389  get_gridpoint_value(matrix, i, j, options);
390  }
391  else
392  {
393  patches[index].data(0, 0) =
394  get_gridpoint_value(matrix, i, j, options);
395  patches[index].data(0, 1) =
396  get_gridpoint_value(matrix, i + 1, j, options);
397  patches[index].data(0, 2) =
398  get_gridpoint_value(matrix, i, j + 1, options);
399  patches[index].data(0, 3) =
400  get_gridpoint_value(matrix, i + 1, j + 1, options);
401  }
402  };
403 
404  // finally set the name
405  this->name = name;
406 }
407 
408 
409 
410 /*---------------------------- matrix_out.h ---------------------------*/
411 
413 
414 #endif
415 /*---------------------------- matrix_out.h ---------------------------*/
virtual const std::vector< Patch > & get_patches() const override
Definition: matrix_out.cc:32
void build_patches(const Matrix &matrix, const std::string &name, const Options options=Options(false, 1, false))
Definition: matrix_out.h:325
static double get_gridpoint_value(const Matrix &matrix, const size_type i, const size_type j, const Options &options)
Definition: matrix_out.h:282
std::vector< Patch > patches
Definition: matrix_out.h:156
types::global_dof_index size_type
Definition: matrix_out.h:77
virtual ~MatrixOut() override=default
std::string name
Definition: matrix_out.h:161
virtual std::vector< std::string > get_dataset_names() const override
Definition: matrix_out.cc:40
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
Point< 3 > vertices[4]
Expression fabs(const Expression &x)
@ matrix
Contents is actually a matrix.
constexpr const ReferenceCell Quadrilateral
double get_element(const ::SparseMatrix< number > &matrix, const types::global_dof_index i, const types::global_dof_index j)
Definition: matrix_out.h:205
unsigned int global_dof_index
Definition: types.h:76
bool show_absolute_values
Definition: matrix_out.h:89
Options(const bool show_absolute_values=false, const unsigned int block_size=1, const bool discontinuous=false)
Definition: matrix_out.cc:21
unsigned int block_size
Definition: matrix_out.h:100