From ea56c9a8e6b03d6407ca011c98c86be3f274bf6b Mon Sep 17 00:00:00 2001
From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Wed, 18 Jul 2001 10:44:57 +0000
Subject: [PATCH] FourierSineSum and FourierCosineSum

git-svn-id: https://svn.dealii.org/trunk@4853 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/base/include/base/function_lib.h | 131 ++++++++++++++
 deal.II/base/source/function_lib.cc      | 208 ++++++++++++++++++++++-
 deal.II/doc/news/2001/c-3-1.html         |  11 +-
 3 files changed, 344 insertions(+), 6 deletions(-)

diff --git a/deal.II/base/include/base/function_lib.h b/deal.II/base/include/base/function_lib.h
index 9c08e2ae25..858ffdd850 100644
--- 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 042e0362df..1917c57624 100644
--- 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 5355aa3020..3004913a0a 100644
--- 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)
 
-- 
2.39.5