From: guido Date: Tue, 15 May 2001 12:22:59 +0000 (+0000) Subject: derivatives of singularity functions X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=044f6525485f8f6cc1f17a4100b24387a6d33074;p=dealii-svn.git derivatives of singularity functions git-svn-id: https://svn.dealii.org/trunk@4606 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/source/function_derivative.cc b/deal.II/base/source/function_derivative.cc index b7fefecbf3..6ee9b01cc9 100644 --- a/deal.II/base/source/function_derivative.cc +++ b/deal.II/base/source/function_derivative.cc @@ -56,7 +56,7 @@ FunctionDerivative::set_formula (DifferenceFormula form) formula = form; } -//TODO:[GK] Discussion on an efficient implementation of Point additions. +//TODO:[?] Discussion on an efficient implementation of Point additions. template double @@ -88,7 +88,7 @@ FunctionDerivative::value (const Point &p, using namespace std; #endif -//TODO:[GK] Optimize construction of vectors! +//TODO:[?] Optimize construction of vectors thread-safe template void diff --git a/deal.II/base/source/function_lib.cc b/deal.II/base/source/function_lib.cc index 10f4e43749..426ec13377 100644 --- a/deal.II/base/source/function_lib.cc +++ b/deal.II/base/source/function_lib.cc @@ -847,30 +847,19 @@ LSingularityFunction::laplacian_list (const std::vector > &points, } -//TODO:[?] Implement derivatives - Tensor<1,2> LSingularityFunction::gradient (const Point<2> &p, const unsigned int) const { - Assert(false, ExcNotImplemented()); double x = p(0); double y = p(1); - - if ((x>=0) && (y>=0)) - { - // in this region, the gradient - // does not exist, so return - // infinity - static const double infty[2] = { atof("Inf"), atof("Inf") }; - return Tensor<1,2>(infty); - } - -//TODO:[GK] this is no reasonable return value -// double phi = atan2(y,-x)+M_PI; -// double r2 = x*x+y*y; - - return Tensor<1,2>(); + double phi = atan2(y,-x)+M_PI; + double r43 = pow(x*x+y*y,2./3.); + + Tensor<1,2> result; + result[0] = 2./3.*(sin(2./3.*phi)*x + cos(2./3.*phi)*y)/r43; + result[1] = 2./3.*(sin(2./3.*phi)*y - cos(2./3.*phi)*x)/r43; + return result; } @@ -881,7 +870,18 @@ LSingularityFunction::gradient_list (const std::vector > &points, { Assert (gradients.size() == points.size(), ExcDimensionMismatch(gradients.size(), points.size())); - Assert(false, ExcNotImplemented()); + + for (unsigned int i=0;i& p = points[i]; + double x = p(0); + double y = p(1); + double phi = atan2(y,-x)+M_PI; + double r43 = pow(x*x+y*y,2./3.); + + gradients[i][0] = 2./3.*(sin(2./3.*phi)*x + cos(2./3.*phi)*y)/r43; + gradients[i][1] = 2./3.*(sin(2./3.*phi)*y - cos(2./3.*phi)*x)/r43; + } } ////////////////////////////////////////////////////////////////////// @@ -894,7 +894,7 @@ SlitSingularityFunction::value (const Point<2> &p, double x = p(0); double y = p(1); - double phi = atan2(x,-y)+M_PI; + double phi = atan2(x,y)+M_PI; double r2 = x*x+y*y; return pow(r2,.25) * sin(.5*phi); @@ -914,7 +914,7 @@ SlitSingularityFunction::value_list (const std::vector > &points, double x = points[i](0); double y = points[i](1); - double phi = atan2(x,-y)+M_PI; + double phi = atan2(x,y)+M_PI; double r2 = x*x+y*y; values[i] = pow(r2,.25) * sin(.5*phi); @@ -943,14 +943,20 @@ SlitSingularityFunction::laplacian_list (const std::vector > &points, } -//TODO:[?] Implement derivatives - Tensor<1,2> -SlitSingularityFunction::gradient (const Point<2> &/*p*/, +SlitSingularityFunction::gradient (const Point<2> &p, const unsigned int) const { - Assert(false, ExcNotImplemented()); - return Tensor<1,2>(); + double x = p(0); + double y = p(1); + double phi = atan2(x,y)+M_PI; + double r64 = pow(x*x+y*y,3./4.); + + + Tensor<1,2> result; + result[0] = 1./2.*(sin(1./2.*phi)*x + cos(1./2.*phi)*y)/r64; + result[1] = 1./2.*(sin(1./2.*phi)*y - cos(1./2.*phi)*x)/r64; + return result; } @@ -961,7 +967,18 @@ SlitSingularityFunction::gradient_list (const std::vector > &points, { Assert (gradients.size() == points.size(), ExcDimensionMismatch(gradients.size(), points.size())); - Assert(false, ExcNotImplemented()); + + for (unsigned int i=0;i& p = points[i]; + double x = p(0); + double y = p(1); + double phi = atan2(x,y)+M_PI; + double r64 = pow(x*x+y*y,3./4.); + + gradients[i][0] = 1./2.*(sin(1./2.*phi)*x + cos(1./2.*phi)*y)/r64; + gradients[i][1] = 1./2.*(sin(1./2.*phi)*y - cos(1./2.*phi)*x)/r64; + } } //////////////////////////////////////////////////////////////////////