deal.II version GIT relicensing-3083-g7b89508ac7 2025-04-18 12:50:00+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\}}\)
Loading...
Searching...
No Matches
matrix_out.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2001 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_matrix_out_h
16# define dealii_matrix_out_h
17
18# include <deal.II/base/config.h>
19
21
24
25# ifdef DEAL_II_WITH_TRILINOS
28# endif
29
30
32
73class MatrixOut : public DataOutInterface<2, 2>
74{
75public:
80
85 struct Options
86 {
96
106 unsigned int block_size;
107
112
117 Options(const bool show_absolute_values = false,
118 const unsigned int block_size = 1,
119 const bool discontinuous = false);
120 };
121
125 virtual ~MatrixOut() override = default;
126
144 template <class Matrix>
145 void
146 build_patches(const Matrix &matrix,
147 const std::string &name,
148 const Options options = Options(false, 1, false));
149
150private:
156
162 std::vector<Patch> patches;
163
167 std::string name;
168
173 virtual const std::vector<Patch> &
174 get_patches() const override;
175
180 virtual std::vector<std::string>
181 get_dataset_names() const override;
182
190 template <class Matrix>
191 static double
192 get_gridpoint_value(const Matrix &matrix,
193 const size_type i,
194 const size_type j,
195 const Options &options);
196};
197
198
199/* ---------------------- Template and inline functions ------------- */
200
201
202namespace internal
203{
204 namespace MatrixOutImplementation
205 {
209 template <typename number>
210 double
211 get_element(const ::SparseMatrix<number> &matrix,
214 {
215 return matrix.el(i, j);
216 }
217
218
219
223 template <typename number>
224 double
225 get_element(const ::BlockSparseMatrix<number> &matrix,
228 {
229 return matrix.el(i, j);
230 }
231
232
233# ifdef DEAL_II_WITH_TRILINOS
237 inline double
241 {
242 return matrix.el(i, j);
243 }
244
245
246
251 inline double
255 {
256 return matrix.el(i, j);
257 }
258# endif
259
260
261# ifdef DEAL_II_WITH_PETSC
262 // no need to do anything: PETSc matrix objects do not distinguish
263 // between operator() and el(i,j), so we can safely access elements
264 // through the generic function below
265# endif
266
267
273 template <class Matrix>
274 double
275 get_element(const Matrix &matrix,
278 {
279 return matrix(i, j);
280 }
281 } // namespace MatrixOutImplementation
282} // namespace internal
283
284
285
286template <class Matrix>
287inline double
289 const size_type i,
290 const size_type j,
291 const Options &options)
292{
293 // special case if block size is
294 // one since we then don't need all
295 // that loop overhead
296 if (options.block_size == 1)
297 {
298 if (options.show_absolute_values == true)
299 return std::fabs(
301 else
303 }
304
305 // if blocksize greater than one,
306 // then compute average of elements
307 double average = 0;
308 size_type n_elements = 0;
309 for (size_type row = i * options.block_size;
310 row <
311 std::min(size_type(matrix.m()), size_type((i + 1) * options.block_size));
312 ++row)
313 for (size_type col = j * options.block_size;
314 col < std::min(size_type(matrix.m()),
315 size_type((j + 1) * options.block_size));
316 ++col, ++n_elements)
317 if (options.show_absolute_values == true)
318 average += std::fabs(
320 else
321 average +=
323 average /= n_elements;
324 return average;
325}
326
327
328
329template <class Matrix>
330void
331MatrixOut::build_patches(const Matrix &matrix,
332 const std::string &name,
333 const Options options)
334{
335 size_type gridpoints_x = (matrix.n() / options.block_size +
336 (matrix.n() % options.block_size != 0 ? 1 : 0)),
337 gridpoints_y = (matrix.m() / options.block_size +
338 (matrix.m() % options.block_size != 0 ? 1 : 0));
339
340 // If continuous, the number of
341 // plotted patches is matrix size-1
342 if (!options.discontinuous)
343 {
344 --gridpoints_x;
345 --gridpoints_y;
346 }
347
348 // first clear old data and re-set the object to a correctly sized state:
349 patches.clear();
350 try
351 {
352 patches.resize(gridpoints_x * gridpoints_y);
353 }
354 catch (const std::bad_alloc &)
355 {
356 AssertThrow(false,
357 ExcMessage("You are trying to create a graphical "
358 "representation of a matrix that would "
359 "requiring outputting " +
360 std::to_string(gridpoints_x) + "x" +
361 std::to_string(gridpoints_y) +
362 " patches. There is not enough memory to output " +
363 "this many patches."));
364 }
365
366 // now build the patches
367 size_type index = 0;
368 for (size_type i = 0; i < gridpoints_y; ++i)
369 for (size_type j = 0; j < gridpoints_x; ++j, ++index)
370 {
371 patches[index].n_subdivisions = 1;
372 patches[index].reference_cell = ReferenceCells::Quadrilateral;
373
374 // within each patch, order the points in such a way that if some
375 // graphical output program (such as gnuplot) plots the quadrilaterals
376 // as two triangles, then the diagonal of the quadrilateral which cuts
377 // it into the two printed triangles is parallel to the diagonal of the
378 // matrix, rather than perpendicular to it. this has the advantage that,
379 // for example, the unit matrix is plotted as a straight ridge, rather
380 // than as a series of bumps and valleys along the diagonal
381 patches[index].vertices[0][0] = j;
382 patches[index].vertices[0][1] = -static_cast<signed int>(i);
383 patches[index].vertices[1][0] = j;
384 patches[index].vertices[1][1] = -static_cast<signed int>(i + 1);
385 patches[index].vertices[2][0] = j + 1;
386 patches[index].vertices[2][1] = -static_cast<signed int>(i);
387 patches[index].vertices[3][0] = j + 1;
388 patches[index].vertices[3][1] = -static_cast<signed int>(i + 1);
389 // next scale all the patch
390 // coordinates by the block
391 // size, to get original
392 // coordinates
393 for (auto &vertex : patches[index].vertices)
394 vertex *= options.block_size;
395
396 patches[index].n_subdivisions = 1;
397
398 patches[index].data.reinit(1, 4);
399 if (options.discontinuous)
400 {
401 patches[index].data(0, 0) =
402 get_gridpoint_value(matrix, i, j, options);
403 patches[index].data(0, 1) =
404 get_gridpoint_value(matrix, i, j, options);
405 patches[index].data(0, 2) =
406 get_gridpoint_value(matrix, i, j, options);
407 patches[index].data(0, 3) =
408 get_gridpoint_value(matrix, i, j, options);
409 }
410 else
411 {
412 patches[index].data(0, 0) =
413 get_gridpoint_value(matrix, i, j, options);
414 patches[index].data(0, 1) =
415 get_gridpoint_value(matrix, i + 1, j, options);
416 patches[index].data(0, 2) =
417 get_gridpoint_value(matrix, i, j + 1, options);
418 patches[index].data(0, 3) =
419 get_gridpoint_value(matrix, i + 1, j + 1, options);
420 }
421 };
422
423 // finally set the name
424 this->name = name;
425}
426
427
428
429/*---------------------------- matrix_out.h ---------------------------*/
430
432
433#endif
434/*---------------------------- matrix_out.h ---------------------------*/
virtual const std::vector< Patch > & get_patches() const override
Definition matrix_out.cc:35
void build_patches(const Matrix &matrix, const std::string &name, const Options options=Options(false, 1, false))
Definition matrix_out.h:331
static double get_gridpoint_value(const Matrix &matrix, const size_type i, const size_type j, const Options &options)
Definition matrix_out.h:288
std::vector< Patch > patches
Definition matrix_out.h:162
types::global_dof_index size_type
Definition matrix_out.h:79
virtual ~MatrixOut() override=default
std::string name
Definition matrix_out.h:167
virtual std::vector< std::string > get_dataset_names() const override
Definition matrix_out.cc:43
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
constexpr 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:211
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
Definition types.h:94
unsigned int block_size
Definition matrix_out.h:106