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
34 /**
35 * Output a matrix in graphical form using the generic format independent
36 * output routines of the base class. The matrix is converted into a list of
37 * patches on a 2d domain where the height is given by the elements of the
38 * matrix. The functions of the base class can then write this "mountain
39 * representation" of the matrix in a variety of graphical output formats. The
40 * coordinates of the matrix output are that the columns run with increasing
41 * x-axis, as usual, starting from zero, while the rows run into the negative
42 * y-axis, also starting from zero. Note that due to some internal
43 * restrictions, this class can only output one matrix at a time, i.e. it can
44 * not take advantage of the multiple dataset capabilities of the base class.
45 *
46 * A typical usage of this class would be as follows:
47 * @code
48 * FullMatrix<double> M;
49 * // fill matrix M with some values
50 * ...
51 *
52 * // now write out M:
53 * MatrixOut matrix_out;
54 * std::ofstream out ("M.gnuplot");
55 * matrix_out.build_patches (M, "M");
56 * matrix_out.write_gnuplot (out);
57 * @endcode
58 * Of course, you can as well choose a different graphical output format.
59 * Also, this class supports any matrix, not only of type FullMatrix, as long
60 * as it satisfies a number of requirements, stated with the member functions
61 * of this class.
62 *
63 * The generation of patches through the build_patches() function can be
64 * modified by giving it an object holding certain flags. See the
65 * documentation of the members of the Options class for a description of
66 * these flags.
67 *
68 *
69 * @ingroup output
70 */
71 class MatrixOut : public DataOutInterface<2, 2>
72 {
73 public:
74 /**
75 * Declare type for container size.
76 */
77 using size_type = types::global_dof_index;
78
79 /**
80 * Class holding various variables which are used to modify the output of
81 * the MatrixOut class.
82 */
83 struct Options
84 {
85 /**
86 * If @p true, only show the absolute values of the matrix entries, rather
87 * than their true values including the sign. Default value is @p false.
88 */
89 bool show_absolute_values;
90
91 /**
92 * If larger than one, do not show each element of the matrix, but rather
93 * an average over a number of entries. The number of output patches is
94 * accordingly smaller. This flag determines how large each shown block
95 * shall be (in rows/columns). For example, if it is two, then always four
96 * entries are collated into one.
97 *
98 * Default value is one.
99 */
100 unsigned int block_size;
101
102 /**
103 * If true, plot discontinuous patches, one for each entry.
104 */
105 bool discontinuous;
106
107 /**
108 * Default constructor. Set all elements of this structure to their
109 * default values.
110 */
111 Options(const bool show_absolute_values = false,
112 const unsigned int block_size = 1,
113 const bool discontinuous = false);
114 };
115
116 /**
117 * Destructor. Declared in order to make it virtual.
118 */
119 virtual ~MatrixOut() override = default;
120
121 /**
122 * Generate a list of patches from the given matrix and use the given string
123 * as the name of the data set upon writing to a file. Once patches have
124 * been built, you can use the functions of the base class to write the data
125 * into a files, using one of the supported output formats.
126 *
127 * You may give a structure holding various options. See the description of
128 * the fields of this structure for more information.
129 *
130 * Note that this function requires that we can extract elements of the
131 * matrix, which is done using the get_element() function declared in an
132 * internal namespace. By adding specializations, you can extend this class
133 * to other matrix classes which are not presently supported. Furthermore,
134 * we need to be able to extract the size of the matrix, for which we assume
135 * that the matrix type offers member functions <tt>m()</tt> and
136 * <tt>n()</tt>, which return the number of rows and columns, respectively.
137 */
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:
145 /**
146 * Abbreviate the somewhat lengthy name for the dealii::DataOutBase::Patch
147 * class.
148 */
149 using Patch = DataOutBase::Patch<2, 2>;
150
151 /**
152 * This is a list of patches that is created each time build_patches() is
153 * called. These patches are used in the output routines of the base
154 * classes.
155 */
156 std::vector<Patch> patches;
157
158 /**
159 * Name of the matrix to be written.
160 */
161 std::string name;
162
163 /**
164 * Function by which the base class's functions get to know what patches
165 * they shall write to a file.
166 */
167 virtual const std::vector<Patch> &
168 get_patches() const override;
169
170 /**
171 * Virtual function through which the names of data sets are obtained by the
172 * output functions of the base class.
173 */
174 virtual std::vector<std::string>
175 get_dataset_names() const override;
176
177 /**
178 * Get the value of the matrix at gridpoint <tt>(i,j)</tt>. Depending on the
179 * given flags, this can mean different things, for example if only absolute
180 * values shall be shown then the absolute value of the matrix entry is
181 * taken. If the block size is larger than one, then an average of several
182 * matrix entries is taken.
183 */
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 {
200 /**
201 * Return the element with given indices of a sparse matrix.
202 */
203 template <typename number>
204 double
get_element(const dealii::SparseMatrix<number> & matrix,const types::global_dof_index i,const types::global_dof_index j)205 get_element(const dealii::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
214 /**
215 * Return the element with given indices of a block sparse matrix.
216 */
217 template <typename number>
218 double
get_element(const dealii::BlockSparseMatrix<number> & matrix,const types::global_dof_index i,const types::global_dof_index j)219 get_element(const dealii::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
228 /**
229 * Return the element with given indices of a Trilinos sparse matrix.
230 */
231 inline double
get_element(const TrilinosWrappers::SparseMatrix & matrix,const types::global_dof_index i,const types::global_dof_index j)232 get_element(const TrilinosWrappers::SparseMatrix &matrix,
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
241 /**
242 * Return the element with given indices of a Trilinos block sparse
243 * matrix.
244 */
245 inline double
get_element(const TrilinosWrappers::BlockSparseMatrix & matrix,const types::global_dof_index i,const types::global_dof_index j)246 get_element(const TrilinosWrappers::BlockSparseMatrix &matrix,
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
262 /**
263 * Return the element with given indices from any matrix type for which
264 * no specialization of this function was declared above. This will call
265 * <tt>operator()</tt> on the matrix.
266 */
267 template <class Matrix>
268 double
get_element(const Matrix & matrix,const types::global_dof_index i,const types::global_dof_index j)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
get_gridpoint_value(const Matrix & matrix,const size_type i,const size_type j,const Options & options)282 MatrixOut::get_gridpoint_value(const Matrix & matrix,
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(
294 internal::MatrixOutImplementation::get_element(matrix, i, j));
295 else
296 return internal::MatrixOutImplementation::get_element(matrix, i, j);
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(
313 internal::MatrixOutImplementation::get_element(matrix, row, col));
314 else
315 average +=
316 internal::MatrixOutImplementation::get_element(matrix, row, col);
317 average /= n_elements;
318 return average;
319 }
320
321
322
323 template <class Matrix>
324 void
build_patches(const Matrix & matrix,const std::string & name,const Options options)325 MatrixOut::build_patches(const Matrix & matrix,
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
409 DEAL_II_NAMESPACE_CLOSE
410
411 #endif
412 /*---------------------------- matrix_out.h ---------------------------*/
413