00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__matrix_out_h
00013 #define __deal2__matrix_out_h
00014
00015 #include <deal.II/base/config.h>
00016 #include <deal.II/base/data_out_base.h>
00017 #include <deal.II/lac/sparse_matrix.h>
00018 #include <deal.II/lac/block_sparse_matrix.h>
00019
00020 DEAL_II_NAMESPACE_OPEN
00021
00088 class MatrixOut : public DataOutInterface<2,2>
00089 {
00090 public:
00097 struct Options
00098 {
00108 bool show_absolute_values;
00109
00127 unsigned int block_size;
00128
00134 bool discontinuous;
00135
00142 Options (const bool show_absolute_values = false,
00143 const unsigned int block_size = 1,
00144 const bool discontinuous = false);
00145 };
00146
00151 virtual ~MatrixOut ();
00152
00190 template <class Matrix>
00191 void build_patches (const Matrix &matrix,
00192 const std::string &name,
00193 const Options options = Options(false, 1, false));
00194
00195 private:
00196
00212 typedef ::DataOutBase::Patch<2,2> Patch;
00213
00222 std::vector<Patch> patches;
00223
00228 std::string name;
00229
00236 virtual const std::vector<Patch> &
00237 get_patches () const;
00238
00245 virtual std::vector<std::string> get_dataset_names () const;
00246
00251 template <typename number>
00252 static double get_element (const SparseMatrix<number> &matrix,
00253 const unsigned int i,
00254 const unsigned int j);
00255
00261 template <typename number>
00262 static double get_element (const BlockSparseMatrix<number> &matrix,
00263 const unsigned int i,
00264 const unsigned int j);
00265
00274 template <class Matrix>
00275 static double get_element (const Matrix &matrix,
00276 const unsigned int i,
00277 const unsigned int j);
00278
00292 template <class Matrix>
00293 static double get_gridpoint_value (const Matrix &matrix,
00294 const unsigned int i,
00295 const unsigned int j,
00296 const Options &options);
00297 };
00298
00299
00300
00301
00302
00303 template <typename number>
00304 inline
00305 double
00306 MatrixOut::get_element (const SparseMatrix<number> &matrix,
00307 const unsigned int i,
00308 const unsigned int j)
00309 {
00310 return matrix.el(i,j);
00311 }
00312
00313
00314
00315
00316 template <typename number>
00317 inline
00318 double
00319 MatrixOut::get_element (const BlockSparseMatrix<number> &matrix,
00320 const unsigned int i,
00321 const unsigned int j)
00322 {
00323 return matrix.el(i,j);
00324 }
00325
00326
00327
00328
00329 template <class Matrix>
00330 inline
00331 double
00332 MatrixOut::get_element (const Matrix &matrix,
00333 const unsigned int i,
00334 const unsigned int j)
00335 {
00336 return matrix(i,j);
00337 }
00338
00339
00340
00341 template <class Matrix>
00342 inline
00343 double
00344 MatrixOut::get_gridpoint_value (const Matrix &matrix,
00345 const unsigned int i,
00346 const unsigned int j,
00347 const Options &options)
00348 {
00349
00350
00351
00352 if (options.block_size == 1)
00353 {
00354 if (options.show_absolute_values == true)
00355 return std::fabs(get_element (matrix, i, j));
00356 else
00357 return get_element (matrix, i, j);
00358 };
00359
00360
00361
00362 double average = 0;
00363 unsigned int n_elements = 0;
00364 for (unsigned int row=i*options.block_size;
00365 row < std::min(matrix.m(), (i+1)*options.block_size); ++row)
00366 for (unsigned int col=j*options.block_size;
00367 col < std::min(matrix.m(), (j+1)*options.block_size); ++col, ++n_elements)
00368 if (options.show_absolute_values == true)
00369 average += std::fabs(get_element (matrix, row, col));
00370 else
00371 average += get_element (matrix, row, col);
00372 average /= n_elements;
00373 return average;
00374 }
00375
00376
00377
00378 template <class Matrix>
00379 void
00380 MatrixOut::build_patches (const Matrix &matrix,
00381 const std::string &name,
00382 const Options options)
00383 {
00384 unsigned int
00385 gridpoints_x = (matrix.n() / options.block_size
00386 +
00387 (matrix.n() % options.block_size != 0 ? 1 : 0)),
00388 gridpoints_y = (matrix.m() / options.block_size
00389 +
00390 (matrix.m() % options.block_size != 0 ? 1 : 0));
00391
00392
00393
00394 if (!options.discontinuous)
00395 {
00396 --gridpoints_x;
00397 --gridpoints_y;
00398 }
00399
00400
00401
00402 patches.clear ();
00403 patches.resize ((gridpoints_x) * (gridpoints_y));
00404
00405
00406 unsigned int index=0;
00407 for (unsigned int i=0; i<gridpoints_y; ++i)
00408 for (unsigned int j=0; j<gridpoints_x; ++j, ++index)
00409 {
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430 patches[index].vertices[0](0) = j;
00431 patches[index].vertices[0](1) = static_cast<signed int>(-i);
00432 patches[index].vertices[1](0) = j;
00433 patches[index].vertices[1](1) = static_cast<signed int>(-i-1);
00434 patches[index].vertices[2](0) = j+1;
00435 patches[index].vertices[2](1) = static_cast<signed int>(-i);
00436 patches[index].vertices[3](0) = j+1;
00437 patches[index].vertices[3](1) = static_cast<signed int>(-i-1);
00438
00439
00440
00441
00442 for (unsigned int v=0; v<4; ++v)
00443 patches[index].vertices[v] *= options.block_size;
00444
00445 patches[index].n_subdivisions = 1;
00446
00447 patches[index].data.reinit (1,4);
00448 if (options.discontinuous)
00449 {
00450 patches[index].data(0,0) = get_gridpoint_value(matrix, i, j, options);
00451 patches[index].data(0,1) = get_gridpoint_value(matrix, i, j, options);
00452 patches[index].data(0,2) = get_gridpoint_value(matrix, i, j, options);
00453 patches[index].data(0,3) = get_gridpoint_value(matrix, i, j, options);
00454 } else {
00455 patches[index].data(0,0) = get_gridpoint_value(matrix, i, j, options);
00456 patches[index].data(0,1) = get_gridpoint_value(matrix, i+1, j, options);
00457 patches[index].data(0,2) = get_gridpoint_value(matrix, i, j+1, options);
00458 patches[index].data(0,3) = get_gridpoint_value(matrix, i+1, j+1, options);
00459 }
00460 };
00461
00462
00463 this->name = name;
00464 }
00465
00466
00467
00468
00469
00470 DEAL_II_NAMESPACE_CLOSE
00471
00472 #endif
00473
00474