From: kanschat Date: Tue, 16 Nov 2010 15:33:27 +0000 (+0000) Subject: add Bessel function X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ac48d43b4018c483044eecd1a73cfd9564027a29;p=dealii-svn.git add Bessel function git-svn-id: https://svn.dealii.org/trunk@22754 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/configure.in b/deal.II/configure.in index 240b1c98fc..4b9a0bf7fa 100644 --- a/deal.II/configure.in +++ b/deal.II/configure.in @@ -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 index 0000000000..501fe6fa3b --- /dev/null +++ b/deal.II/include/deal.II/base/function_bessel.h @@ -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 +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Functions +{ +/** + * The Bessel functions of first kind or positive integer order. + * + * @author Guido Kanschat + * @date 2010 + */ + template + class Bessel1 : public Function + { + public: + Bessel1(const unsigned int order, + const double wave_number, + const Point center = Point()); + virtual double value (const Point& points, const unsigned int component) const; + virtual void value_list (const std::vector > &points, + std::vector &values, + const unsigned int component = 0) const; + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component = 0) const; + virtual void gradient_list (const std::vector > &points, + std::vector > &gradients, + const unsigned int component = 0) const; + private: + unsigned int order; + double wave_number; + Point center; + }; +} + +DEAL_II_NAMESPACE_CLOSE + +#endif + diff --git a/deal.II/source/base/function_lib.cc b/deal.II/source/base/function_lib.cc index fdcf8fe1aa..6081795600 100644 --- a/deal.II/source/base/function_lib.cc +++ b/deal.II/source/base/function_lib.cc @@ -14,6 +14,7 @@ #include #include #include +#include #include #include @@ -2155,9 +2156,112 @@ namespace Functions for (unsigned int i=0; i::value (points[i], component); } + + + template + Bessel1::Bessel1( + const unsigned int order, + const double wave_number, + const Point center) + : + order(order), + wave_number(wave_number), + center(center) + {} + + template + double + Bessel1::value(const Point& 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 + void + Bessel1::value_list ( + const std::vector > &points, + std::vector &values, + const unsigned int) const + { + Assert(dim==2, ExcNotImplemented()); + AssertDimension(points.size(), values.size()); + for (unsigned int k=0;k + Tensor<1,dim> + Bessel1::gradient (const Point &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 + void + Bessel1::gradient_list ( + const std::vector > &points, + std::vector > &gradients, + const unsigned int) const + { + Assert(dim==2, ExcNotImplemented()); + AssertDimension(points.size(), gradients.size()); + for (unsigned int k=0;k& 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