Add gradient of Cinfty cutoff function. Reindent slightly. Add TODO.

git-svn-id: https://svn.dealii.org/trunk@5150 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 941e52b523922c662f5035d0a09844a4506e2418..0cd930c318fdb5a60c5524515564a29ebf2a8bd6 100644 (file)
--- a/deal.II/base/include/base/function_lib.h
+++ b/deal.II/base/include/base/function_lib.h
@@ -877,7 +877,8 @@ namespace Functions
 /**
  * Cut-off function for an arbitrary ball. This is the traditional
  * cut-off function in C-infinity for a ball of certain @p{radius}
- * around @p{center}.
+ * around @p{center}, $f(r)=exp(1-1/(1-r**2/s**2))$, where $r$ is the
+ * distance to the center, and $s$ is the radius of the sphere.
  *
  * @author Guido Kanschat, 2001
  */
@@ -906,6 +907,12 @@ namespace Functions
                               std::vector<double>            &values,
                               const unsigned int              component = 0) const;
 
+                                      /**
+                                       * Function gradient at one point.
+                                       */
+      virtual Tensor<1,dim> gradient (const Point<dim>   &p,
+                                     const unsigned int  component = 0) const;
+      
     private:
                                       /**
                                        * Center of the integration ball.

diff --git a/deal.II/base/source/function_lib_cutoff.cc b/deal.II/base/source/function_lib_cutoff.cc
index 771234705058aa7bbe239328c550d2ef946a1781..e5a2b3ab12fb2ff17a8eb8339f9887ab6da2c248 100644 (file)
--- a/deal.II/base/source/function_lib_cutoff.cc
+++ b/deal.II/base/source/function_lib_cutoff.cc
@@ -21,6 +21,7 @@
 
 // in strict ANSI C mode, the following constants are not defined by
 // default, so we do it ourselves
+//TODO:[?] Unify the various places where PI is defined to a central instance 
 #ifndef M_PI
 #  define      M_PI            3.14159265358979323846
 #endif
@@ -41,9 +42,10 @@ namespace Functions
   template<int dim>
   CutOffFunctionLinfty<dim>::CutOffFunctionLinfty (const double r,
                                                   const Point<dim> p)
-    : Function<dim> (1),
-    center(p),
-    radius(r)
+                 :
+                 Function<dim> (1),
+                  center(p),
+                  radius(r)
   {}
 
 
@@ -74,16 +76,17 @@ namespace Functions
   template<int dim>
   CutOffFunctionW1<dim>::CutOffFunctionW1 (const double     r,
                                           const Point<dim> p)
-    : Function<dim> (1),
-    center(p),
-    radius(r)
+                 :
+                 Function<dim> (1),
+                  center(p),
+                  radius(r)
   {}
 
 
   template<int dim>
   double
   CutOffFunctionW1<dim>::value (const Point<dim>   &p,
-                             const unsigned int) const
+                               const unsigned int) const
   {
     const double d = center.distance(p);
     return ((d<radius) ? (radius-d) : 0.);
@@ -111,16 +114,17 @@ namespace Functions
   template<int dim>
   CutOffFunctionCinfty<dim>::CutOffFunctionCinfty (const double     r,
                                                   const Point<dim> p)
-    : Function<dim> (1),
-    center(p),
-    radius(r)
+                 :
+                 Function<dim> (1),
+                  center(p),
+                  radius(r)
   {}
 
 
   template<int dim>
   double
   CutOffFunctionCinfty<dim>::value (const Point<dim>   &p,
-                             const unsigned int) const
+                                   const unsigned int) const
   {
     const double d = center.distance(p);
     const double r = radius;
@@ -156,6 +160,24 @@ namespace Functions
   }
 
 
+
+  template<int dim>
+  Tensor<1,dim>
+  CutOffFunctionCinfty<dim>::gradient (const Point<dim>   &p,
+                                      const unsigned int) const
+  {
+    const double d = center.distance(p);
+    const double r = radius;
+    if (d>=r)
+      return Tensor<1,dim>();
+    const double e = -d*d/(r-d)/(r+d);
+    return  ((e<-50) ?
+            Point<dim>() :
+            (p-center)/d*(-2.0*r*r/pow(-r*r+d*d,2.0)*d*exp(e)));
+  }
+  
+
+// explicit instantiations  
   template class CutOffFunctionLinfty <1>;
   template class CutOffFunctionLinfty <2>;
   template class CutOffFunctionLinfty <3>;
@@ -167,6 +189,4 @@ namespace Functions
   template class CutOffFunctionCinfty <1>;
   template class CutOffFunctionCinfty <2>;
   template class CutOffFunctionCinfty <3>;
-
-
 }