Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
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
282MatrixOut::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(
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
325MatrixOut::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 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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
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
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
Definition types.h:82
unsigned int block_size
Definition matrix_out.h:100