include/deal.II/lac/matrix_out.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: matrix_out.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
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 /* ---------------------- Template and inline functions ------------- */
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                                    // special case if block size is
00350                                    // one since we then don't need all
00351                                    // that loop overhead
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                                    // if blocksize greater than one,
00361                                    // then compute average of elements
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                                    // If continuous, the number of
00393                                    // plotted patches is matrix size-1
00394   if (!options.discontinuous)
00395     {
00396       --gridpoints_x;
00397       --gridpoints_y;
00398     }
00399 
00400                                    // first clear old data and set it
00401                                    // to virgin state
00402   patches.clear ();
00403   patches.resize ((gridpoints_x) * (gridpoints_y));
00404 
00405                                    // now build the patches
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                                          // within each patch, order
00411                                          // the points in such a way
00412                                          // that if some graphical
00413                                          // output program (such as
00414                                          // gnuplot) plots the
00415                                          // quadrilaterals as two
00416                                          // triangles, then the
00417                                          // diagonal of the
00418                                          // quadrilateral which cuts
00419                                          // it into the two printed
00420                                          // triangles is parallel to
00421                                          // the diagonal of the
00422                                          // matrix, rather than
00423                                          // perpendicular to it. this
00424                                          // has the advantage that,
00425                                          // for example, the unit
00426                                          // matrix is plotted as a
00427                                          // straight rim, rather than
00428                                          // as a series of bumps and
00429                                          // valleys along the diagonal
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                                          // next scale all the patch
00439                                          // coordinates by the block
00440                                          // size, to get original
00441                                          // coordinates
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                                    // finally set the name
00463   this->name = name;
00464 }
00465 
00466 
00467 
00468 /*----------------------------   matrix_out.h     ---------------------------*/
00469 
00470 DEAL_II_NAMESPACE_CLOSE
00471 
00472 #endif
00473 /*----------------------------   matrix_out.h     ---------------------------*/
00474 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Tue May 22 2012 12:06:11 by doxygen 1.7.3