From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Tue, 16 Oct 2001 12:46:22 +0000 (+0000)
Subject: Add gradient of Cinfty cutoff function. Reindent slightly. Add TODO.
X-Git-Tag: v8.0.0~18673
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4a7d35221bf735aebe236fdf148f91bc84d00db0;p=dealii.git

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 941e52b523..0cd930c318 100644
--- 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 7712347050..e5a2b3ab12 100644
--- 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>;
-
-
 }