Loading [MathJax]/extensions/TeX/AMSsymbols.js
 Reference documentation for deal.II version 9.3.3
\(\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\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
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
71class MatrixOut : public DataOutInterface<2, 2>
72{
73public:
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
144private:
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
196namespace internal
197{
198 namespace MatrixOutImplementation
199 {
203 template <typename number>
204 double
205 get_element(const ::SparseMatrix<number> &matrix,
208 {
209 return matrix.el(i, j);
210 }
211
212
213
217 template <typename number>
218 double
219 get_element(const ::BlockSparseMatrix<number> &matrix,
222 {
223 return matrix.el(i, j);
224 }
225
226
227# ifdef DEAL_II_WITH_TRILINOS
231 inline double
235 {
236 return matrix.el(i, j);
237 }
238
239
240
245 inline double
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,
272 {
273 return matrix(i, j);
274 }
275 } // namespace MatrixOutImplementation
276} // namespace internal
277
278
279
280template <class Matrix>
281inline 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
323template <class Matrix>
324void
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 // within each patch, order the points in such a way that if some
353 // graphical output program (such as gnuplot) plots the quadrilaterals
354 // as two triangles, then the diagonal of the quadrilateral which cuts
355 // it into the two printed triangles is parallel to the diagonal of the
356 // matrix, rather than perpendicular to it. this has the advantage that,
357 // for example, the unit matrix is plotted as a straight rim, rather
358 // than as a series of bumps and valleys along the diagonal
359 patches[index].vertices[0](0) = j;
360 patches[index].vertices[0](1) = -static_cast<signed int>(i);
361 patches[index].vertices[1](0) = j;
362 patches[index].vertices[1](1) = -static_cast<signed int>(i + 1);
363 patches[index].vertices[2](0) = j + 1;
364 patches[index].vertices[2](1) = -static_cast<signed int>(i);
365 patches[index].vertices[3](0) = j + 1;
366 patches[index].vertices[3](1) = -static_cast<signed int>(i + 1);
367 // next scale all the patch
368 // coordinates by the block
369 // size, to get original
370 // coordinates
371 for (auto &vertex : patches[index].vertices)
372 vertex *= options.block_size;
373
374 patches[index].n_subdivisions = 1;
375
376 patches[index].data.reinit(1, 4);
377 if (options.discontinuous)
378 {
379 patches[index].data(0, 0) =
380 get_gridpoint_value(matrix, i, j, options);
381 patches[index].data(0, 1) =
382 get_gridpoint_value(matrix, i, j, options);
383 patches[index].data(0, 2) =
384 get_gridpoint_value(matrix, i, j, options);
385 patches[index].data(0, 3) =
386 get_gridpoint_value(matrix, i, j, options);
387 }
388 else
389 {
390 patches[index].data(0, 0) =
391 get_gridpoint_value(matrix, i, j, options);
392 patches[index].data(0, 1) =
393 get_gridpoint_value(matrix, i + 1, j, options);
394 patches[index].data(0, 2) =
395 get_gridpoint_value(matrix, i, j + 1, options);
396 patches[index].data(0, 3) =
397 get_gridpoint_value(matrix, i + 1, j + 1, options);
398 }
399 };
400
401 // finally set the name
402 this->name = name;
403}
404
405
406
407/*---------------------------- matrix_out.h ---------------------------*/
408
410
411#endif
412/*---------------------------- 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
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:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
Point< 3 > vertices[4]
Expression fabs(const Expression &x)
@ matrix
Contents is actually a matrix.
double get_element(const ::SparseMatrix< number > &matrix, const types::global_dof_index i, const types::global_dof_index j)
Definition: matrix_out.h:205
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
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