]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
add Bessel function
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 16 Nov 2010 15:33:27 +0000 (15:33 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 16 Nov 2010 15:33:27 +0000 (15:33 +0000)
git-svn-id: https://svn.dealii.org/trunk@22754 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/configure.in
deal.II/include/deal.II/base/function_bessel.h [new file with mode: 0644]
deal.II/source/base/function_lib.cc

index 240b1c98fcda00c763a551e80129c45ae698c44b..4b9a0bf7fabc16290cced90779c414486d4de6dc 100644 (file)
@@ -624,6 +624,8 @@ AC_CHECK_FUNCS([dgeqrf_ sgeqrf_ dormqr_ sormqr_ dorgqr_ sorgqr_ dtrtrs_ strtrs_]
 AC_CHECK_FUNCS([dgetrs_ sgetrs_ dstev_ sstev_ dsygv_ ssygv_])
 dnl Singular value decomposition
 AC_CHECK_FUNCS([dgesvd_ sgesvd_ dgesdd_ sgesdd_ dgelsd_ sgelsd_])
+dnl Check Bessel functions in GNU libc
+AC_CHECK_FUNCS([jn])
 
 dnl -------------------------------------------------------------
 dnl       set include paths of several libraries
diff --git a/deal.II/include/deal.II/base/function_bessel.h b/deal.II/include/deal.II/base/function_bessel.h
new file mode 100644 (file)
index 0000000..501fe6f
--- /dev/null
@@ -0,0 +1,56 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//
+//    Copyright (C) 2010 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__function_bessel_h
+#define __deal2__function_bessel_h
+
+
+#include <base/config.h>
+#include <base/function.h>
+#include <base/point.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Functions
+{
+/**
+ * The Bessel functions of first kind or positive integer order.
+ *
+ * @author Guido Kanschat
+ * @date 2010
+ */
+  template <int dim>
+  class Bessel1 : public Function<dim>
+  {
+    public:
+      Bessel1(const unsigned int order,
+             const double wave_number,
+             const Point<dim> center = Point<dim>());
+      virtual double value (const Point<dim>& points, const unsigned int component) const;
+      virtual void value_list (const std::vector<Point<dim> > &points,
+                              std::vector<double>            &values,
+                              const unsigned int              component = 0) const;
+      virtual Tensor<1,dim> gradient (const Point<dim>   &p,
+                                     const unsigned int  component = 0) const;
+      virtual void gradient_list (const std::vector<Point<dim> > &points,
+                                 std::vector<Tensor<1,dim> >    &gradients,
+                                 const unsigned int              component = 0) const;
+    private:
+      unsigned int order;
+      double wave_number;
+      Point<dim> center;
+  };
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
+
index fdcf8fe1aafcd79fab2820596bbf7bf086b9ab57..608179560025764c4a9f8ecb84a1c07c056f66a2 100644 (file)
@@ -14,6 +14,7 @@
 #include <base/tensor.h>
 #include <base/point.h>
 #include <base/function_lib.h>
+#include <base/function_bessel.h>
 #include <lac/vector.h>
 
 #include <cmath>
@@ -2155,9 +2156,112 @@ namespace Functions
     for (unsigned int i=0; i<points.size(); ++i)
       values[i] = Monomial<dim>::value (points[i], component);
   }
+  
+  
+  template <int dim>
+  Bessel1<dim>::Bessel1(
+    const unsigned int order,
+    const double wave_number,
+    const Point<dim> center)
+                 :
+                 order(order),
+                 wave_number(wave_number),
+                 center(center)
+  {}
+  
+  template <int dim>
+  double
+  Bessel1<dim>::value(const Point<dim>& p, const unsigned int) const
+  {
+    Assert(dim==2, ExcNotImplemented());
+    const double r = p.distance(center);
+#ifdef HAVE_JN_
+    return jn(order, r*wave_number);
+#else
+    Assert(false, ExcMessage("Bessel function jn was not found by configure"));
+    return r;
+#endif
+  }
+  
+  
+  template <int dim>
+  void
+  Bessel1<dim>::value_list (
+    const std::vector<Point<dim> > &points,
+    std::vector<double>            &values,
+    const unsigned int) const
+  {
+    Assert(dim==2, ExcNotImplemented());
+    AssertDimension(points.size(), values.size());
+    for (unsigned int k=0;k<points.size();++k)
+      {
+#ifdef HAVE_JN_
+       const double r = points[k].distance(center);
+       values[k] = jn(order, r*wave_number);
+#else
+    Assert(false, ExcMessage("Bessel function jn was not found by configure"));
+#endif    
+      }
+  }
 
   
+  template <int dim>
+  Tensor<1,dim>
+  Bessel1<dim>::gradient (const Point<dim>   &p,
+                         const unsigned int) const
+  {
+    Assert(dim==2, ExcNotImplemented());
+    const double r = p.distance(center);
+    const double co = (r==0.) ? 0. : (p(0)-center(0))/r;
+    const double si = (r==0.) ? 0. : (p(1)-center(1))/r;
+    
+    const double dJn = (order==0)
+                      ? (-jn(1, r*wave_number))
+                      : (.5*(jn(order-1, wave_number*r) -jn(order+1, wave_number*r)));
+    Tensor<1,dim> result;
+    result[0] = wave_number * co * dJn;
+    result[1] = wave_number * si * dJn;
+#ifdef HAVE_JN_
+    return result;
+#else
+    Assert(false, ExcMessage("Bessel function jn was not found by configure"));
+    return Tensor<1,dim>();
+#endif
+  }
+  
+
+  
+  template <int dim>
+  void
+  Bessel1<dim>::gradient_list (
+    const std::vector<Point<dim> > &points,
+    std::vector<Tensor<1,dim> >    &gradients,
+    const unsigned int) const
+  {
+    Assert(dim==2, ExcNotImplemented());
+    AssertDimension(points.size(), gradients.size());
+    for (unsigned int k=0;k<points.size();++k)
+      {
+       const Point<dim>& p = points[k];
+       const double r = p.distance(center);
+       const double co = (r==0.) ? 0. : (p(0)-center(0))/r;
+       const double si = (r==0.) ? 0. : (p(1)-center(1))/r;
+       
+#ifdef HAVE_JN_
+       const double dJn = (order==0)
+                          ? (-jn(1, r*wave_number))
+                          : (.5*(jn(order-1, wave_number*r) -jn(order+1, wave_number*r)));
+#else
+       const double dJn = 0.;
+       Assert(false, ExcMessage("Bessel function jn was not found by configure"));
+#endif
+       Tensor<1,dim>& result = gradients[k];
+       result[0] = wave_number * co * dJn;
+       result[1] = wave_number * si * dJn;
+      }
+  }
   
+
 // explicit instantiations  
   template class SquareFunction<1>;
   template class SquareFunction<2>;
@@ -2197,6 +2301,8 @@ namespace Functions
   template class Monomial<1>;
   template class Monomial<2>;
   template class Monomial<3>;
+  template class Bessel1<2>;
+  template class Bessel1<3>;
 }
 
 DEAL_II_NAMESPACE_CLOSE

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.