FourierSineSum and FourierCosineSum

git-svn-id: https://svn.dealii.org/trunk@4853 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/base/include/base/function_lib.h b/deal.II/base/include/base/function_lib.h
index 9c08e2ae2571b2a78dfb0b0113e9859a8d272a97..858ffdd8507b3c68d7992eff1ebd0a74811f77ca 100644 (file)
--- a/deal.II/base/include/base/function_lib.h
+++ b/deal.II/base/include/base/function_lib.h
@@ -652,6 +652,137 @@ namespace Functions
       const Point<dim> fourier_coefficients;
   };
   
+
+/**
+ * Given a sequence of wavenumber vectors and weights generate a sum
+ * of sine functions. Each wavenumber coefficient is given as a
+ * @p{d}-dimensional point @p{k} in Fourier space, and the entire
+ * function is then recovered as
+ * @p{f(x) = \sum_j w_j \prod_i sin(k_i x_i) = Im(\sum_j w_j \exp(i k.x))}.
+ *
+ * @author Wolfgang Bangerth, 2001
+ */
+  template <int dim>
+  class FourierSineSum : public Function<dim> 
+  {
+    public:
+                                      /**
+                                       * Constructor. Take the Fourier
+                                       * coefficients in each space
+                                       * direction as argument.
+                                       */
+      FourierSineSum (const std::vector<Point<dim> > &fourier_coefficients,
+                     const std::vector<double>      &weights);
+      
+                                      /**
+                                       * Return the value of the
+                                       * function at the given
+                                       * point. Unless there is only
+                                       * one component (i.e. the
+                                       * function is scalar), you
+                                       * should state the component you
+                                       * want to have evaluated; it
+                                       * defaults to zero, i.e. the
+                                       * first component.
+                                       */
+      virtual double value (const Point<dim>   &p,
+                           const unsigned int  component = 0) const;
+      
+                                      /**
+                                       * Return the gradient of the
+                                       * specified component of the
+                                       * function at the given point.
+                                       */
+      virtual Tensor<1,dim> gradient (const Point<dim>   &p,
+                                     const unsigned int  component = 0) const;
+      
+                                      /**
+                                       * Compute the Laplacian of a
+                                       * given component at point @p{p}.
+                                       */
+      virtual double laplacian (const Point<dim>   &p,
+                               const unsigned int  component = 0) const;
+
+                                      /**
+                                       * Exception
+                                       */
+      DeclException0 (ExcInvalidArraySize);
+      
+    private:
+                                      /**
+                                       * Stored Fourier coefficients
+                                       * and weights.
+                                       */
+      const std::vector<Point<dim> > fourier_coefficients;
+      const std::vector<double>      weights;
+  };
+
+
+
+/**
+ * Given a sequence of wavenumber vectors and weights generate a sum
+ * of cosine functions. Each wavenumber coefficient is given as a
+ * @p{d}-dimensional point @p{k} in Fourier space, and the entire
+ * function is then recovered as
+ * @p{f(x) = \sum_j w_j \prod_i cos(k_i x_i) = Re(\sum_j w_j \exp(i k.x))}.
+ *
+ * @author Wolfgang Bangerth, 2001
+ */
+  template <int dim>
+  class FourierCosineSum : public Function<dim> 
+  {
+    public:
+                                      /**
+                                       * Constructor. Take the Fourier
+                                       * coefficients in each space
+                                       * direction as argument.
+                                       */
+      FourierCosineSum (const std::vector<Point<dim> > &fourier_coefficients,
+                       const std::vector<double>      &weights);
+      
+                                      /**
+                                       * Return the value of the
+                                       * function at the given
+                                       * point. Unless there is only
+                                       * one component (i.e. the
+                                       * function is scalar), you
+                                       * should state the component you
+                                       * want to have evaluated; it
+                                       * defaults to zero, i.e. the
+                                       * first component.
+                                       */
+      virtual double value (const Point<dim>   &p,
+                           const unsigned int  component = 0) const;
+      
+                                      /**
+                                       * Return the gradient of the
+                                       * specified component of the
+                                       * function at the given point.
+                                       */
+      virtual Tensor<1,dim> gradient (const Point<dim>   &p,
+                                     const unsigned int  component = 0) const;
+      
+                                      /**
+                                       * Compute the Laplacian of a
+                                       * given component at point @p{p}.
+                                       */
+      virtual double laplacian (const Point<dim>   &p,
+                               const unsigned int  component = 0) const;
+
+                                      /**
+                                       * Exception
+                                       */
+      DeclException0 (ExcInvalidArraySize);
+      
+    private:
+                                      /**
+                                       * Stored Fourier coefficients
+                                       * and weights.
+                                       */
+      const std::vector<Point<dim> > fourier_coefficients;
+      const std::vector<double>      weights;
+  };
+  
   
 };
 

diff --git a/deal.II/base/source/function_lib.cc b/deal.II/base/source/function_lib.cc
index 042e0362dfd86478353e841387662c4be16e970f..1917c57624b1984bb5e363b37cd58e66abfce549 100644 (file)
--- a/deal.II/base/source/function_lib.cc
+++ b/deal.II/base/source/function_lib.cc
@@ -1255,9 +1255,209 @@ namespace Functions
     return val;
   };
   
+
+
+
+/* ---------------------- FourierSineSum ----------------------- */
+  
+  
+  
+  template <int dim>
+  FourierSineSum<dim>::
+  FourierSineSum (const std::vector<Point<dim> > &fourier_coefficients,
+                 const std::vector<double>      &weights)
+                 :
+                 Function<dim> (1),
+    fourier_coefficients (fourier_coefficients),
+    weights (weights)
+  {
+    Assert (fourier_coefficients.size() > 0, ExcInvalidArraySize());
+    Assert (fourier_coefficients.size() == weights.size(),
+           ExcInvalidArraySize());
+  };
+  
+  
+  
+  template <int dim>
+  double
+  FourierSineSum<dim>::value (const Point<dim>   &p,
+                             const unsigned int  component) const
+  {
+    Assert (component==0, ExcIndexRange(component,0,1));
+
+    const unsigned int n = weights.size();
+    double sum = 0;
+    for (unsigned int s=0; s<n; ++s)
+      {
+       double val=1;
+       for (unsigned int i=0; i<dim; ++i)
+         val *= std::sin(fourier_coefficients[s][i] * p[i]);
+       sum += weights[s] * val;
+      };
+    
+    return sum;
+  };
+  
+  
+  
+  template <int dim>
+  Tensor<1,dim>
+  FourierSineSum<dim>::gradient (const Point<dim>   &p,
+                                const unsigned int  component) const
+  {
+    Assert (component==0, ExcIndexRange(component,0,1));
+
+    const unsigned int n = weights.size();
+    Tensor<1,dim> sum;
+    for (unsigned int s=0; s<n; ++s)
+      {
+       Tensor<1,dim> grad;
+       for (unsigned int i=0; i<dim; ++i)
+         grad[i] = 1;
+    
+       for (unsigned int i=0; i<dim; ++i)
+         {
+           const double cos_i = std::cos(fourier_coefficients[s][i] * p[i]);
+           const double sin_i = std::sin(fourier_coefficients[s][i] * p[i]);
+           
+           for (unsigned int d=0; d<dim; ++d)
+             if (d==i)
+               grad[d] *= fourier_coefficients[s][i] * cos_i;
+             else
+               grad[d] *= sin_i;
+         };
+       
+       grad *= weights[s];
+       sum += grad;
+      };
+    
+    return sum;
+  };
+  
+  
+  
+  template <int dim>
+  double
+  FourierSineSum<dim>::laplacian (const Point<dim>   &p,
+                                 const unsigned int  component) const
+  {
+    Assert (component==0, ExcIndexRange(component,0,1));
+
+    const unsigned int n = weights.size();
+    double sum = 0;
+    for (unsigned int s=0; s<n; ++s)
+      {
+       double val = -(fourier_coefficients[s]*fourier_coefficients[s]);
+       for (unsigned int i=0; i<dim; ++i)
+         val *= std::sin(fourier_coefficients[s][i] * p[i]);
+       sum += val * weights[s];
+      };
+    
+    return sum;
+  };
+
+
+
+/* ---------------------- FourierCosineSum ----------------------- */
+  
+  
+  
+  template <int dim>
+  FourierCosineSum<dim>::
+  FourierCosineSum (const std::vector<Point<dim> > &fourier_coefficients,
+                 const std::vector<double>      &weights)
+                 :
+                 Function<dim> (1),
+    fourier_coefficients (fourier_coefficients),
+    weights (weights)
+  {
+    Assert (fourier_coefficients.size() > 0, ExcInvalidArraySize());
+    Assert (fourier_coefficients.size() == weights.size(),
+           ExcInvalidArraySize());
+  };
+  
+  
+  
+  template <int dim>
+  double
+  FourierCosineSum<dim>::value (const Point<dim>   &p,
+                               const unsigned int  component) const
+  {
+    Assert (component==0, ExcIndexRange(component,0,1));
+
+    const unsigned int n = weights.size();
+    double sum = 0;
+    for (unsigned int s=0; s<n; ++s)
+      {
+       double val=1;
+       for (unsigned int i=0; i<dim; ++i)
+         val *= std::cos(fourier_coefficients[s][i] * p[i]);
+       sum += weights[s] * val;
+      };
+    
+    return sum;
+  };
+  
+  
+  
+  template <int dim>
+  Tensor<1,dim>
+  FourierCosineSum<dim>::gradient (const Point<dim>   &p,
+                                  const unsigned int  component) const
+  {
+    Assert (component==0, ExcIndexRange(component,0,1));
+
+    const unsigned int n = weights.size();
+    Tensor<1,dim> sum;
+    for (unsigned int s=0; s<n; ++s)
+      {
+       Tensor<1,dim> grad;
+       for (unsigned int i=0; i<dim; ++i)
+         grad[i] = 1;
+    
+       for (unsigned int i=0; i<dim; ++i)
+         {
+           const double cos_i = std::cos(fourier_coefficients[s][i] * p[i]);
+           const double sin_i = std::sin(fourier_coefficients[s][i] * p[i]);
+           
+           for (unsigned int d=0; d<dim; ++d)
+             if (d==i)
+               grad[d] *= -fourier_coefficients[s][i] * sin_i;
+             else
+               grad[d] *= cos_i;
+         };
+       
+       grad *= weights[s];
+       sum += grad;
+      };
+    
+    return sum;
+  };
   
   
   
+  template <int dim>
+  double
+  FourierCosineSum<dim>::laplacian (const Point<dim>   &p,
+                                   const unsigned int  component) const
+  {
+    Assert (component==0, ExcIndexRange(component,0,1));
+
+    const unsigned int n = weights.size();
+    double sum = 0;
+    for (unsigned int s=0; s<n; ++s)
+      {
+       double val = -(fourier_coefficients[s]*fourier_coefficients[s]);
+       for (unsigned int i=0; i<dim; ++i)
+         val *= std::cos(fourier_coefficients[s][i] * p[i]);
+       sum += val * weights[s];
+      };
+    
+    return sum;
+  };
+  
+  
+// explicit instantiations  
   template class SquareFunction<1>;
   template class SquareFunction<2>;
   template class SquareFunction<3>;
@@ -1282,6 +1482,10 @@ namespace Functions
   template class FourierSineFunction<1>;
   template class FourierSineFunction<2>;
   template class FourierSineFunction<3>;
-  
-  
+  template class FourierCosineSum<1>;
+  template class FourierCosineSum<2>;
+  template class FourierCosineSum<3>;
+  template class FourierSineSum<1>;
+  template class FourierSineSum<2>;
+  template class FourierSineSum<3>;
 };

diff --git a/deal.II/doc/news/2001/c-3-1.html b/deal.II/doc/news/2001/c-3-1.html
index 5355aa30201ed8c114ff0594623fd6f51fcdee15..3004913a0a10013f35cc6b0625ec486566ac18f5 100644 (file)
--- a/deal.II/doc/news/2001/c-3-1.html
+++ b/deal.II/doc/news/2001/c-3-1.html
@@ -168,10 +168,13 @@ documentation, etc</a>.
        <br>
        (WB 2001/07/18)
 
-  <li> <p>
-       New: classes <code class="class">FourierSineFunction</code> and
-       <code class="class">FourierCosineFunction</code>, resembling
-       one mode of a Fourier decomposition.
+  <li> <p> 
+       New: classes <code class="class">FourierSineFunction</code>
+       and <code class="class">FourierCosineFunction</code>,
+       resembling one mode of a Fourier decomposition. Classes <code
+       class="class">FourierSineSum</code> and <code
+       class="class">FourierCosineSum</code>, resembling sums of such
+       modes of a Fourier decomposition.
        <br>
        (WB 2001/07/18)