new JumpFunction

git-svn-id: https://svn.dealii.org/trunk@3177 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 85c68bc4f3827fb530f207b411dbd136fe6f9e0c..3b6c838659a7b2a0c55e0f233f53f46e09aa46f6 100644 (file)
--- a/deal.II/base/include/base/function_lib.h
+++ b/deal.II/base/include/base/function_lib.h
@@ -32,39 +32,39 @@ class SquareFunction : public Function<dim>
 {
   public:
                                     /**
-                                     * @reimplemented
+                                     * Function value at one point.
                                      */
     virtual double value (const Point<dim>   &p,
                          const unsigned int  component = 0) const;
 
                                     /**
-                                     * @reimplemented
+                                     * Function values at multiple points.
                                      */
     virtual void value_list (const vector<Point<dim> > &points,
                             vector<double>            &values,
                             const unsigned int         component = 0) const;
 
                                     /**
-                                     * @reimplemented
+                                     * Gradient at one point.
                                      */
     virtual Tensor<1,dim> gradient (const Point<dim>   &p,
                                    const unsigned int  component = 0) const;
 
                                     /**
-                                     * @reimplemented
+                                       Gradients at multiple points.
                                      */
     virtual void gradient_list (const vector<Point<dim> > &points,
                                vector<Tensor<1,dim> >    &gradients,
                                const unsigned int         component = 0) const;
 
                                     /**
-                                     * @reimplemented
+                                     * Laplacian of the function at one point.
                                      */
     double laplacian (const Point<dim>   &p,
                      const unsigned int  component = 0) const;
 
                                     /**
-                                     * @reimplemented
+                                     * Laplacian of the function at multiple points.
                                      */
     void laplacian_list (const vector<Point<dim> > &points,
                         vector<double>            &values,
@@ -329,4 +329,97 @@ class SlitSingularityFunction : public Function<2>
 };
 
 
+/**
+ * A jump in x-direction transported into some direction.
+ *
+ * If the advection is parallel to the y-axis, the function is
+ * @p{-atan(sx)}, where @p{s} is the steepness parameter provided in
+ * the constructor.
+ *
+ * For different advection directions, this function will be turned in
+ * the parameter space.
+ *
+ * Together with the function, its derivatives and Laplacian are defined.
+ *
+ * @author: Guido Kanschat, 2000
+ */
+template<int dim>
+class JumpFunction : public Function<dim>
+{
+  public:
+                                    /**
+                                     * Constructor. Provide the
+                                     * advection direction here and
+                                     * the steepness of the slope.
+                                     */
+    JumpFunction (const Point<dim>& direction, double steepness);
+    
+                                    /**
+                                     * Function value at one point.
+                                     */
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;
+
+                                    /**
+                                     * Function values at multiple points.
+                                     */
+    virtual void value_list (const vector<Point<dim> > &points,
+                            vector<double>            &values,
+                            const unsigned int         component = 0) const;
+
+                                    /**
+                                     * Gradient at one point.
+                                     */
+    virtual Tensor<1,dim> gradient (const Point<dim>   &p,
+                                   const unsigned int  component = 0) const;
+
+                                    /**
+                                       Gradients at multiple points.
+                                     */
+    virtual void gradient_list (const vector<Point<dim> > &points,
+                               vector<Tensor<1,dim> >    &gradients,
+                               const unsigned int         component = 0) const;
+
+                                    /**
+                                     * Laplacian of the function at one point.
+                                     */
+    double laplacian (const Point<dim>   &p,
+                     const unsigned int  component = 0) const;
+
+                                    /**
+                                     * Laplacian of the function at multiple points.
+                                     */
+    void laplacian_list (const vector<Point<dim> > &points,
+                        vector<double>            &values,
+                        const unsigned int         component = 0) const;
+
+  protected:
+                                    /**
+                                     * Advection vector.
+                                     */
+    Point<dim> direction;
+                                    /**
+                                     * Steepness (maximal derivative) of the slope.
+                                     */
+    double steepness;
+
+                                    /**
+                                     * Advection angle.
+                                     */
+    double angle;
+
+                                    /**
+                                     * Sine of @p{angle}.
+                                     */
+    double sine;
+
+                                    /**
+                                     * Cosine of @p{angle}.
+                                     */
+    double cosine;
+};
+
+
+
+
 #endif

diff --git a/deal.II/base/source/function_lib.cc b/deal.II/base/source/function_lib.cc
index 1f3845a959dbc9ca1295872655f2627a15cff2be..eef1c4d20e5f7ee880294d5dd1c0e4610cb07f83 100644 (file)
--- a/deal.II/base/source/function_lib.cc
+++ b/deal.II/base/source/function_lib.cc
@@ -763,6 +763,124 @@ SlitSingularityFunction::gradient_list (const vector<Point<2> > &points,
 
 //////////////////////////////////////////////////////////////////////
 
+template<int dim>
+JumpFunction<dim>::JumpFunction(const Point<dim>& direction,
+                               double steepness)
+               :
+               direction(direction),
+               steepness(steepness)
+{
+  switch (dim)
+    {
+      case 1:
+           angle = 0;
+           break;
+      case 2:
+           angle = atan2(direction(0),direction(1));
+           break;
+      case 3:
+           Assert(false, ExcNotImplemented());
+    }
+  sine = sin(angle);
+  cosine = cos(angle);
+}
+
+           
+
+template<int dim>
+double
+JumpFunction<dim>::value (const Point<dim>   &p,
+                           const unsigned int) const
+{
+  double x = steepness*(cosine*p(0)+sine*p(1));
+  return -atan(x);
+}
+
+
+
+template<int dim>
+void
+JumpFunction<dim>::value_list (const vector<Point<dim> > &p,
+                                vector<double>          &values,
+                                const unsigned int) const
+{
+  Assert (values.size() == p.size(),
+         ExcDimensionMismatch(values.size(), p.size()));
+
+  for (unsigned int i=0;i<p.size();++i)
+    {
+      double x = steepness*(cosine*p[i](0)+sine*p[i](1));
+      values[i] = -atan(x);
+    }
+}
+
+
+template<int dim>
+double
+JumpFunction<dim>::laplacian (const Point<dim>   &p,
+                           const unsigned int) const
+{
+  double x = steepness*(cosine*p(0)+sine*p(1));
+  double r = 1+x*x;
+  return 2*steepness*steepness*x/(r*r);
+}
+
+
+template<int dim>
+void
+JumpFunction<dim>::laplacian_list (const vector<Point<dim> > &p,
+                                    vector<double>          &values,
+                                    const unsigned int) const
+{
+  Assert (values.size() == p.size(),
+         ExcDimensionMismatch(values.size(), p.size()));
+
+  double f = 2*steepness*steepness;
+  
+  for (unsigned int i=0;i<p.size();++i)
+    {
+      double x = steepness*(cosine*p[i](0)+sine*p[i](1));
+      double r = 1+x*x;
+      values[i] = f*x/(r*r);
+    }
+}
+
+
+
+template<int dim>
+Tensor<1,dim>
+JumpFunction<dim>::gradient (const Point<dim>   &p,
+                              const unsigned int) const
+{
+  double x = steepness*(cosine*p(0)+sine*p(1));
+  double r = -steepness*(1+x*x);
+  Tensor<1,dim> erg;
+  erg[0] = cosine*r;
+  erg[1] = sine*r;
+  return erg;
+}
+
+
+
+template<int dim>
+void
+JumpFunction<dim>::gradient_list (const vector<Point<dim> > &p,
+                                   vector<Tensor<1,dim> >  &gradients,
+                                   const unsigned int) const
+{
+  Assert (gradients.size() == p.size(),
+         ExcDimensionMismatch(gradients.size(), p.size()));
+
+  for (unsigned int i=0; i<p.size(); ++i)
+    {
+      double x = steepness*(cosine*p[i](0)+sine*p[i](1));
+      double r = -steepness*(1+x*x);
+      gradients[i][0] = cosine*r;
+      gradients[i][1] = sine*r;
+    }
+}
+
+
 template SquareFunction<1>;
 template SquareFunction<2>;
 template SquareFunction<3>;
@@ -775,3 +893,6 @@ template CosineFunction<3>;
 template ExpFunction<1>;
 template ExpFunction<2>;
 template ExpFunction<3>;
+template JumpFunction<1>;
+template JumpFunction<2>;
+template JumpFunction<3>;