]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Introduce Utilities::fixed_power. Use it in data_out_base. Move another exception...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 26 Jun 2007 05:04:09 +0000 (05:04 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 26 Jun 2007 05:04:09 +0000 (05:04 +0000)
git-svn-id: https://svn.dealii.org/trunk@14806 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/utilities.h
deal.II/base/source/data_out_base.cc
deal.II/doc/news/changes.html

index 422d07e2b6a50dd0deb77827ce4df40f9698eb45..af5ac3472a59bb68d2625e593ec5563f89c21e70 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2005, 2006 by the deal.II authors
+//    Copyright (C) 2005, 2006, 2007 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -14,6 +14,7 @@
 #define __deal2__utilities_h
 
 #include <base/config.h>
+#include <base/exceptions.h>
 
 #include <vector>
 #include <utility>
@@ -141,6 +142,22 @@ namespace Utilities
                                 const double sigma);
 
 
+                                  /**
+                                   * Calculate a fixed power, provided as a
+                                   * template argument, of a number.
+                                   *
+                                   * This function provides an efficient way
+                                   * to calculate things like
+                                   * <code>t^N</code> where <code>N</code> is
+                                   * a known number at compile time.
+                                   *
+                                   * Use this function as in
+                                   * <code>fixed_power@<dim@> (n)</code>.
+                                   */
+  template <int N, typename T>
+  T
+  fixed_power (const T t);
+  
                                    /**
                                     * A namespace for utility functions that
                                     * probe system properties.
@@ -208,6 +225,37 @@ namespace Utilities
   }
 }
 
+
+// --------------------- inline functions
+
+namespace Utilities
+{
+  template <int N, typename T>
+  inline
+  T fixed_power (const T n)
+  {
+    Assert (N>0, ExcNotImplemented());
+    switch (N)
+      {
+       case 1:
+             return n;
+       case 2:
+             return n*n;
+       case 3:
+             return n*n*n;
+       case 4:
+             return n*n*n*n;
+       default:
+             T result = n;  
+             for (unsigned int d=1;d<N;++d)
+               result *= n;
+             return result;
+      }
+  }
+  
+}
+
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
index 5df8eb3e9fd1d8dea139b7c1487385072f6333aa..f7ef41004898cf31cd90fe974422861941a10e70 100644 (file)
@@ -30,6 +30,7 @@
 
 
 #include <base/data_out_base.h>
+#include <base/utilities.h>
 #include <base/parameter_handler.h>
 #include <base/thread_management.h>
 #include <base/memory_consumption.h>
 
 DEAL_II_NAMESPACE_OPEN
 
-//TODO: Is it reasonable to have undocumented exceptions?
 
-DeclException2 (ExcUnexpectedInput,
-                std::string, std::string,
-                << "Unexpected input: expected line\n  <"
-                << arg1
-                << ">\nbut got\n  <"
-                << arg2 << ">");
+// we need the following exception from a global function, so can't declare it
+// in the usual way inside a class
+namespace
+{
+  DeclException2 (ExcUnexpectedInput,
+                 std::string, std::string,
+                 << "Unexpected input: expected line\n  <"
+                 << arg1
+                 << ">\nbut got\n  <"
+                 << arg2 << ">");
+}
 
-DeclException4 (ExcIncompatibleDimensions,
-                int, int, int, int,
-                << "Either the dimensions <" << arg1 << "> and <"
-                << arg2 << "> or the space dimensions <"
-                << arg3 << "> and <" << arg4
-                << "> do not match!");
 
 
 //----------------------------------------------------------------------//
 //Auxiliary data
 //----------------------------------------------------------------------//
-static const char* gmv_cell_type[4] =
-{
-      "", "line 2", "quad 4", "hex 8"
-};
 
-static const char* ucd_cell_type[4] =
+namespace
 {
-      "", "line", "quad", "hex"
-};
+  
+  static const char* gmv_cell_type[4] =
+  {
+       "", "line 2", "quad 4", "hex 8"
+  };
 
-static const char* tecplot_cell_type[4] =
-{
-      "", "", "quadrilateral", "brick"
-};
+  static const char* ucd_cell_type[4] =
+  {
+       "", "line", "quad", "hex"
+  };
+
+  static const char* tecplot_cell_type[4] =
+  {
+       "", "", "quadrilateral", "brick"
+  };
 
 #ifdef DEAL_II_HAVE_TECPLOT
-static unsigned int tecplot_binary_cell_type[4] =
-{
-      0, 0, 1, 3
-};
+  static unsigned int tecplot_binary_cell_type[4] =
+  {
+       0, 0, 1, 3
+  };
 #endif
 
-static unsigned int vtk_cell_type[4] =
-{
-      0, 3, 9, 12
-};
+  static unsigned int vtk_cell_type[4] =
+  {
+       0, 3, 9, 12
+  };
 
 //----------------------------------------------------------------------//
 //Auxiliary functions
@@ -99,104 +102,80 @@ static unsigned int vtk_cell_type[4] =
 //linearly at the point (xstep, ystep, zstep)*1./n_subdivisions.
 //If the popints are saved in the patch->data member, return the
 //saved point instead
-template <int dim, int spacedim>
-inline
-void
-compute_node(
-  Point<spacedim>& node,
-  const DataOutBase::Patch<dim,spacedim>* patch,
-  const unsigned int xstep,
-  const unsigned int ystep,
-  const unsigned int zstep,
-  const unsigned int n_subdivisions)
-{
-  if (patch->points_are_available)
-    {
-      unsigned int point_no=0;
-                                      // note: switch without break !
-      switch (dim)
-       {
-         case 3:
-               Assert (zstep<n_subdivisions+1, ExcIndexRange(zstep,0,n_subdivisions+1));
-               point_no+=(n_subdivisions+1)*(n_subdivisions+1)*zstep;
-         case 2:
-               Assert (ystep<n_subdivisions+1, ExcIndexRange(ystep,0,n_subdivisions+1));
-               point_no+=(n_subdivisions+1)*ystep;
-         case 1:
-               Assert (xstep<n_subdivisions+1, ExcIndexRange(xstep,0,n_subdivisions+1));
-               point_no+=xstep;
-       }
-      for (unsigned int d=0; d<spacedim; ++d)
-       node[d]=patch->data(patch->data.size(0)-spacedim+d,point_no);
-    }
-  else
-    {
-                                      // perform a dim-linear interpolation
-      const double stepsize=1./n_subdivisions,
-                     xfrac=xstep*stepsize;
-
-      node = (patch->vertices[1] * xfrac) + (patch->vertices[0] * (1-xfrac));
-      if (dim>1)
-       {
-         const double yfrac=ystep*stepsize;
-         node*= 1-yfrac;
-         node += ((patch->vertices[3] * xfrac) + (patch->vertices[2] * (1-xfrac))) * yfrac;
-         if (dim>2)
-           {
-             const double zfrac=zstep*stepsize;
-             node *= (1-zfrac);
-             node += (((patch->vertices[5] * xfrac) + (patch->vertices[4] * (1-xfrac)))
-                      * (1-yfrac) +
-                      ((patch->vertices[7] * xfrac) + (patch->vertices[6] * (1-xfrac)))
-                      * yfrac) * zfrac;
-           }
-       }
-    }
-}
+  template <int dim, int spacedim>
+  inline
+  void
+  compute_node(
+    Point<spacedim>& node,
+    const DataOutBase::Patch<dim,spacedim>* patch,
+    const unsigned int xstep,
+    const unsigned int ystep,
+    const unsigned int zstep,
+    const unsigned int n_subdivisions)
+  {
+    if (patch->points_are_available)
+      {
+       unsigned int point_no=0;
+                                        // note: switch without break !
+       switch (dim)
+         {
+           case 3:
+                 Assert (zstep<n_subdivisions+1, ExcIndexRange(zstep,0,n_subdivisions+1));
+                 point_no+=(n_subdivisions+1)*(n_subdivisions+1)*zstep;
+           case 2:
+                 Assert (ystep<n_subdivisions+1, ExcIndexRange(ystep,0,n_subdivisions+1));
+                 point_no+=(n_subdivisions+1)*ystep;
+           case 1:
+                 Assert (xstep<n_subdivisions+1, ExcIndexRange(xstep,0,n_subdivisions+1));
+                 point_no+=xstep;
+         }
+       for (unsigned int d=0; d<spacedim; ++d)
+         node[d]=patch->data(patch->data.size(0)-spacedim+d,point_no);
+      }
+    else
+      {
+                                        // perform a dim-linear interpolation
+       const double stepsize=1./n_subdivisions,
+                       xfrac=xstep*stepsize;
 
+       node = (patch->vertices[1] * xfrac) + (patch->vertices[0] * (1-xfrac));
+       if (dim>1)
+         {
+           const double yfrac=ystep*stepsize;
+           node*= 1-yfrac;
+           node += ((patch->vertices[3] * xfrac) + (patch->vertices[2] * (1-xfrac))) * yfrac;
+           if (dim>2)
+             {
+               const double zfrac=zstep*stepsize;
+               node *= (1-zfrac);
+               node += (((patch->vertices[5] * xfrac) + (patch->vertices[4] * (1-xfrac)))
+                        * (1-yfrac) +
+                        ((patch->vertices[7] * xfrac) + (patch->vertices[6] * (1-xfrac)))
+                        * yfrac) * zfrac;
+             }
+         }
+      }
+  }
 
-//TODO: Should this go to the utilities namespace?
-
-// Compute n^dim, where dim is a template parameter
-template <int dim>
-inline
-unsigned int template_power (const unsigned int n)
-{
-  Assert (dim>0, ExcNotImplemented());
-  switch (dim)
-    {
-      case 1:
-           return n;
-      case 2:
-           return n*n;
-      case 3:
-           return n*n*n;
-      case 4:
-           return n*n*n*n;
-      default:
-           unsigned int result = n;  
-           for (unsigned int d=1;d<dim;++d)
-             result *= n;
-           return result;
-    }
-}
 
 
-template<int dim, int spacedim>
-static
-void
-compute_sizes(const std::vector<DataOutBase::Patch<dim, spacedim> >& patches,
-             unsigned int& n_nodes,
-             unsigned int& n_cells)
-{
-  n_nodes = 0;
-  n_cells = 0;
-  for (typename std::vector<DataOutBase::Patch<dim,spacedim> >::const_iterator patch=patches.begin();
-       patch!=patches.end(); ++patch)
-    {
-      n_nodes += template_power<dim>(patch->n_subdivisions+1);
-      n_cells += template_power<dim>(patch->n_subdivisions);
-    }
+  template<int dim, int spacedim>
+  static
+  void
+  compute_sizes(const std::vector<DataOutBase::Patch<dim, spacedim> >& patches,
+               unsigned int& n_nodes,
+               unsigned int& n_cells)
+  {
+    n_nodes = 0;
+    n_cells = 0;
+    for (typename std::vector<DataOutBase::Patch<dim,spacedim> >::const_iterator patch=patches.begin();
+        patch!=patches.end(); ++patch)
+      {
+       n_nodes += Utilities::fixed_power<dim>(patch->n_subdivisions+1);
+       n_cells += Utilities::fixed_power<dim>(patch->n_subdivisions);
+      }
+  }
+  
 }
 
   
@@ -1305,7 +1284,7 @@ void DataOutBase::write_cells(
            }
                                       // finally update the number
                                       // of the first vertex of this patch
-      first_vertex_of_patch += template_power<dim>(n_subdivisions+1);
+      first_vertex_of_patch += Utilities::fixed_power<dim>(n_subdivisions+1);
     }
 }
 
@@ -1331,7 +1310,7 @@ DataOutBase::write_data (
       Assert ((patch->data.n_rows() == n_data_sets && !patch->points_are_available) ||
              (patch->data.n_rows() == n_data_sets+spacedim && patch->points_are_available),
              ExcDimensionMismatch (patch->points_are_available ? (patch->data.n_rows() + spacedim) : patch->data.n_rows(), n_data_sets));
-      Assert (patch->data.n_cols() == template_power<dim>(n),
+      Assert (patch->data.n_cols() == Utilities::fixed_power<dim>(n),
              ExcInvalidDatasetSize (patch->data.n_cols(), n));
       
       std::vector<float>  floats(n_data_sets);
@@ -1339,7 +1318,7 @@ DataOutBase::write_data (
 
                                       // Data is already in
                                       // lexicographic ordering
-      for (unsigned int i=0; i<template_power<dim>(n); ++i)
+      for (unsigned int i=0; i<Utilities::fixed_power<dim>(n); ++i)
        if (double_precision)
          {
            for (unsigned int data_set=0; data_set<n_data_sets; ++data_set)
@@ -1523,7 +1502,7 @@ void DataOutBase::write_dx (const std::vector<Patch<dim,spacedim> > &patches,
          const unsigned int n1 = (dim>0) ? n : 1;
          const unsigned int n2 = (dim>1) ? n : 1;
          const unsigned int n3 = (dim>2) ? n : 1;
-         unsigned int cells_per_patch = template_power<dim>(n);
+         unsigned int cells_per_patch = Utilities::fixed_power<dim>(n);
          unsigned int dx = 1;
          unsigned int dy = n;
          unsigned int dz = n*n;
@@ -1775,7 +1754,7 @@ void DataOutBase::write_gnuplot (const std::vector<Patch<dim,spacedim> > &patche
       Assert ((patch->data.n_rows() == n_data_sets && !patch->points_are_available) ||
              (patch->data.n_rows() == n_data_sets+spacedim && patch->points_are_available),
              ExcDimensionMismatch (patch->points_are_available ? (patch->data.n_rows() + spacedim) : patch->data.n_rows(), n_data_sets));
-      Assert (patch->data.n_cols() == template_power<dim>(n),
+      Assert (patch->data.n_cols() == Utilities::fixed_power<dim>(n),
              ExcInvalidDatasetSize (patch->data.n_cols(), n_subdivisions+1));
 
       Point<spacedim> this_point;
@@ -1993,7 +1972,7 @@ void DataOutBase::write_povray (const std::vector<Patch<dim,spacedim> > &patches
       Assert ((patch->data.n_rows() == n_data_sets && !patch->points_are_available) ||
              (patch->data.n_rows() == n_data_sets+spacedim && patch->points_are_available),
              ExcDimensionMismatch (patch->points_are_available ? (patch->data.n_rows() + spacedim) : patch->data.n_rows(), n_data_sets));
-      Assert (patch->data.n_cols() == template_power<dim>(n_subdivisions+1),
+      Assert (patch->data.n_cols() == Utilities::fixed_power<dim>(n_subdivisions+1),
              ExcInvalidDatasetSize (patch->data.n_cols(), n_subdivisions+1));
       
       for (unsigned int i=0; i<n_subdivisions+1; ++i)
@@ -2043,7 +2022,7 @@ void DataOutBase::write_povray (const std::vector<Patch<dim,spacedim> > &patches
       Assert ((patch->data.n_rows() == n_data_sets && !patch->points_are_available) ||
              (patch->data.n_rows() == n_data_sets+spacedim && patch->points_are_available),
              ExcDimensionMismatch (patch->points_are_available ? (patch->data.n_rows() + spacedim) : patch->data.n_rows(), n_data_sets));
-      Assert (patch->data.n_cols() == template_power<dim>(n),
+      Assert (patch->data.n_cols() == Utilities::fixed_power<dim>(n),
              ExcInvalidDatasetSize (patch->data.n_cols(), n_subdivisions+1));
 
 
@@ -3240,7 +3219,7 @@ void DataOutBase::write_tecplot_binary (const std::vector<Patch<dim,spacedim> >
 
                                        // finally update the number
                                        // of the first vertex of this patch
-       first_vertex_of_patch += template_power<dim>(n);
+       first_vertex_of_patch += Utilities::fixed_power<dim>(n);
      }
 
 
@@ -3522,7 +3501,7 @@ DataOutBase::write_gmv_reorder_data_vectors (const std::vector<Patch<dim,spacedi
       Assert ((patch->data.n_rows() == n_data_sets && !patch->points_are_available) ||
              (patch->data.n_rows() == n_data_sets+spacedim && patch->points_are_available),
              ExcDimensionMismatch (patch->points_are_available ? (patch->data.n_rows() + spacedim) : patch->data.n_rows(), n_data_sets));
-      Assert (patch->data.n_cols() == template_power<dim>(n_subdivisions+1),
+      Assert (patch->data.n_cols() == Utilities::fixed_power<dim>(n_subdivisions+1),
              ExcInvalidDatasetSize (patch->data.n_cols(), n_subdivisions+1));
       
       for (unsigned int i=0;i<patch->data.n_cols();++i, ++next_value)
index cc516410f3fdd9703b80036a6e882715e1c5a8f0..48d5cf736eae9b55747309fab9e4fe3053fc7d28 100644 (file)
@@ -437,6 +437,15 @@ inconvenience this causes.
 <h3>base</h3>
 
 <ol>
+  <li> <p> New: The function <code>Utilities::fixed_power&lt;n&gt;(q)</code>
+       calculates <code>q</code> to the power of <code>n</code> where
+       <code>n</code> is a constant known at compile time. It allows to
+       calculate powers efficiently at compile time, most often things like a
+       number to the power <code>dim</code>.
+       <br> 
+       (WB 2007/06/23)
+       </p>
+
   <li> <p> New: The deal.II intermediate format has been
        changed. Since files in this format are only meant to be
        processed by exactly the same version of deal.II, this should

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.