deal.II version GIT relicensing-1982-gbb401394a0 2024-10-13 08:40:01+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
70class MatrixOut : public DataOutInterface<2, 2>
71{
72public:
77
82 struct Options
83 {
89
99 unsigned int block_size;
100
105
110 Options(const bool show_absolute_values = false,
111 const unsigned int block_size = 1,
112 const bool discontinuous = false);
113 };
114
118 virtual ~MatrixOut() override = default;
119
137 template <class Matrix>
138 void
139 build_patches(const Matrix &matrix,
140 const std::string &name,
141 const Options options = Options(false, 1, false));
142
143private:
149
155 std::vector<Patch> patches;
156
160 std::string name;
161
166 virtual const std::vector<Patch> &
167 get_patches() const override;
168
173 virtual std::vector<std::string>
174 get_dataset_names() const override;
175
183 template <class Matrix>
184 static double
185 get_gridpoint_value(const Matrix &matrix,
186 const size_type i,
187 const size_type j,
188 const Options &options);
189};
190
191
192/* ---------------------- Template and inline functions ------------- */
193
194
195namespace internal
196{
197 namespace MatrixOutImplementation
198 {
202 template <typename number>
203 double
204 get_element(const ::SparseMatrix<number> &matrix,
207 {
208 return matrix.el(i, j);
209 }
210
211
212
216 template <typename number>
217 double
218 get_element(const ::BlockSparseMatrix<number> &matrix,
221 {
222 return matrix.el(i, j);
223 }
224
225
226# ifdef DEAL_II_WITH_TRILINOS
230 inline double
234 {
235 return matrix.el(i, j);
236 }
237
238
239
244 inline double
248 {
249 return matrix.el(i, j);
250 }
251# endif
252
253
254# ifdef DEAL_II_WITH_PETSC
255 // no need to do anything: PETSc matrix objects do not distinguish
256 // between operator() and el(i,j), so we can safely access elements
257 // through the generic function below
258# endif
259
260
266 template <class Matrix>
267 double
268 get_element(const Matrix &matrix,
271 {
272 return matrix(i, j);
273 }
274 } // namespace MatrixOutImplementation
275} // namespace internal
276
277
278
279template <class Matrix>
280inline double
282 const size_type i,
283 const size_type j,
284 const Options &options)
285{
286 // special case if block size is
287 // one since we then don't need all
288 // that loop overhead
289 if (options.block_size == 1)
290 {
291 if (options.show_absolute_values == true)
292 return std::fabs(
294 else
296 }
297
298 // if blocksize greater than one,
299 // then compute average of elements
300 double average = 0;
301 size_type n_elements = 0;
302 for (size_type row = i * options.block_size;
303 row <
304 std::min(size_type(matrix.m()), size_type((i + 1) * options.block_size));
305 ++row)
306 for (size_type col = j * options.block_size;
307 col < std::min(size_type(matrix.m()),
308 size_type((j + 1) * options.block_size));
309 ++col, ++n_elements)
310 if (options.show_absolute_values == true)
311 average += std::fabs(
313 else
314 average +=
316 average /= n_elements;
317 return average;
318}
319
320
321
322template <class Matrix>
323void
324MatrixOut::build_patches(const Matrix &matrix,
325 const std::string &name,
326 const Options options)
327{
328 size_type gridpoints_x = (matrix.n() / options.block_size +
329 (matrix.n() % options.block_size != 0 ? 1 : 0)),
330 gridpoints_y = (matrix.m() / options.block_size +
331 (matrix.m() % options.block_size != 0 ? 1 : 0));
332
333 // If continuous, the number of
334 // plotted patches is matrix size-1
335 if (!options.discontinuous)
336 {
337 --gridpoints_x;
338 --gridpoints_y;
339 }
340
341 // first clear old data and re-set the object to a correctly sized state:
342 patches.clear();
343 patches.resize((gridpoints_x) * (gridpoints_y));
344
345 // now build the patches
346 size_type index = 0;
347 for (size_type i = 0; i < gridpoints_y; ++i)
348 for (size_type j = 0; j < gridpoints_x; ++j, ++index)
349 {
350 patches[index].n_subdivisions = 1;
351 patches[index].reference_cell = ReferenceCells::Quadrilateral;
352
353 // within each patch, order the points in such a way that if some
354 // graphical output program (such as gnuplot) plots the quadrilaterals
355 // as two triangles, then the diagonal of the quadrilateral which cuts
356 // it into the two printed triangles is parallel to the diagonal of the
357 // matrix, rather than perpendicular to it. this has the advantage that,
358 // for example, the unit matrix is plotted as a straight ridge, rather
359 // than as a series of bumps and valleys along the diagonal
360 patches[index].vertices[0][0] = j;
361 patches[index].vertices[0][1] = -static_cast<signed int>(i);
362 patches[index].vertices[1][0] = j;
363 patches[index].vertices[1][1] = -static_cast<signed int>(i + 1);
364 patches[index].vertices[2][0] = j + 1;
365 patches[index].vertices[2][1] = -static_cast<signed int>(i);
366 patches[index].vertices[3][0] = j + 1;
367 patches[index].vertices[3][1] = -static_cast<signed int>(i + 1);
368 // next scale all the patch
369 // coordinates by the block
370 // size, to get original
371 // coordinates
372 for (auto &vertex : patches[index].vertices)
373 vertex *= options.block_size;
374
375 patches[index].n_subdivisions = 1;
376
377 patches[index].data.reinit(1, 4);
378 if (options.discontinuous)
379 {
380 patches[index].data(0, 0) =
381 get_gridpoint_value(matrix, i, j, options);
382 patches[index].data(0, 1) =
383 get_gridpoint_value(matrix, i, j, options);
384 patches[index].data(0, 2) =
385 get_gridpoint_value(matrix, i, j, options);
386 patches[index].data(0, 3) =
387 get_gridpoint_value(matrix, i, j, options);
388 }
389 else
390 {
391 patches[index].data(0, 0) =
392 get_gridpoint_value(matrix, i, j, options);
393 patches[index].data(0, 1) =
394 get_gridpoint_value(matrix, i + 1, j, options);
395 patches[index].data(0, 2) =
396 get_gridpoint_value(matrix, i, j + 1, options);
397 patches[index].data(0, 3) =
398 get_gridpoint_value(matrix, i + 1, j + 1, options);
399 }
400 };
401
402 // finally set the name
403 this->name = name;
404}
405
406
407
408/*---------------------------- matrix_out.h ---------------------------*/
409
411
412#endif
413/*---------------------------- matrix_out.h ---------------------------*/
virtual const std::vector< Patch > & get_patches() const override
Definition matrix_out.cc:31
void build_patches(const Matrix &matrix, const std::string &name, const Options options=Options(false, 1, false))
Definition matrix_out.h:324
static double get_gridpoint_value(const Matrix &matrix, const size_type i, const size_type j, const Options &options)
Definition matrix_out.h:281
std::vector< Patch > patches
Definition matrix_out.h:155
types::global_dof_index size_type
Definition matrix_out.h:76
virtual ~MatrixOut() override=default
std::string name
Definition matrix_out.h:160
virtual std::vector< std::string > get_dataset_names() const override
Definition matrix_out.cc:39
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
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:204
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
Definition types.h:81
unsigned int block_size
Definition matrix_out.h:99