]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added implementation of CutOffFunctionC1
authorLuca Heltai <luca.heltai@sissa.it>
Thu, 18 Apr 2019 19:42:08 +0000 (21:42 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Tue, 23 Apr 2019 13:00:02 +0000 (15:00 +0200)
include/deal.II/base/function_lib.h
source/base/function_lib_cutoff.cc
tests/base/function_cutoff_01.cc
tests/base/function_cutoff_01.output
tests/base/function_cutoff_02.cc
tests/base/function_cutoff_02.output

index ff38a07bc863b6d69768c9959b4c9bd3c6902d05..ab7635e89aa265e066086d45fd307b9e7561468c 100644 (file)
@@ -915,7 +915,7 @@ namespace Functions
    *
    * This class can also be used for approximated Dirac delta functions. These
    * are special cut-off functions whose integral is always equal to one,
-   * independently on the radius of the supporting ball.
+   * independently of the radius of the supporting ball.
    *
    * @ingroup functions
    * @author Guido Kanschat, 2002, Luca Heltai, 2019.
@@ -931,17 +931,23 @@ namespace Functions
     static const unsigned int no_component = numbers::invalid_unsigned_int;
 
     /**
-     * Constructor. Arguments are the center of the ball and its radius.
-     *
-     * If an argument <tt>select</tt> is given and not -1, the cut-off
-     * function will be non-zero for this component only.
+     * Constructor.
      *
-     * If the argument @p integrate_to_one is set to true, then the value of
-     * the function is rescaled whenever a new radius is set.
+     * @param[in] radius Radius of the ball
+     * @param[in] center Center of the ball
+     * @param[in] n_components Number of components of this function object
+     * @param[in] select If this is different from
+     * CutOffFunctionBase<dim>::no_component, then the function will be non-zero
+     * for this component only
+     * @param[in] integrate_to_one Rescale the value of the function whenever a
+     * new radius is set, to guarantee that the integral is equal to one
+     * @param[in] unitary_integral_value Value of the integral when the radius
+     * is equal to 1.0. Derived classes will need to supply this value, to
+     * guarantee that the rescaling is performed correctly.
      */
     CutOffFunctionBase(
-      const double radius             = 1.,
-      const Point<dim>                = Point<dim>(),
+      const double       radius       = 1.,
+      const Point<dim>   center       = Point<dim>(),
       const unsigned int n_components = 1,
       const unsigned int select       = CutOffFunctionBase<dim>::no_component,
       const bool         integrate_to_one       = false,
@@ -976,8 +982,6 @@ namespace Functions
 
     /**
      * Set the radius of the ball to @p r
-     *
-     * @deprecated Use set_radius() instead.
      */
     virtual void
     set_radius(const double r);
@@ -995,7 +999,7 @@ namespace Functions
     get_radius() const;
 
     /**
-     * Return a boolean indicating if this function integrates to one.
+     * Return a boolean indicating whether this function integrates to one.
      */
     bool
     integrates_to_one() const;
@@ -1204,6 +1208,63 @@ namespace Functions
   };
 
 
+  /**
+   * A cut-off function for an arbitrarily-sized ball that is in the space $C^1$
+   * (i.e., continuously differentiable). This is a cut-off function that is
+   * often used in the literature of the Immersed Boundary Method.
+   *
+   * The expression of the function in radial coordinates is given by
+   * $f(r)=1/2(cos(\pi r/s)+1)$ where $r<s$ is the distance to the center, and
+   * $s$ is the radius of the sphere. If vector valued, it can be restricted to
+   * a single component.
+   *
+   * @ingroup functions
+   * @author Luca Heltai, 2019
+   */
+  template <int dim>
+  class CutOffFunctionC1 : public CutOffFunctionBase<dim>
+  {
+  public:
+    /**
+     * Constructor.
+     */
+    CutOffFunctionC1(
+      const double radius             = 1.,
+      const Point<dim>                = Point<dim>(),
+      const unsigned int n_components = 1,
+      const unsigned int select       = CutOffFunctionBase<dim>::no_component,
+      bool               integrate_to_one = false);
+
+    /**
+     * Function value at one point.
+     */
+    virtual double
+    value(const Point<dim> &p, const unsigned int component = 0) const override;
+
+    /**
+     * Function values at multiple points.
+     */
+    virtual void
+    value_list(const std::vector<Point<dim>> &points,
+               std::vector<double> &          values,
+               const unsigned int             component = 0) const override;
+
+    /**
+     * Function values at multiple points.
+     */
+    virtual void
+    vector_value_list(const std::vector<Point<dim>> &points,
+                      std::vector<Vector<double>> &  values) const override;
+
+    /**
+     * Function gradient at one point.
+     */
+    virtual Tensor<1, dim>
+    gradient(const Point<dim> & p,
+             const unsigned int component = 0) const override;
+  };
+
+
   /**
    * Cut-off function for an arbitrary ball. This is the traditional cut-off
    * function in C-infinity for a ball of certain <tt>radius</tt> around
index ff6ca15b0e884dcf00b2a40bfff418121dd80f0d..87d79d9195670767ff5c44122b1ea752112116ad 100644 (file)
@@ -189,16 +189,29 @@ namespace Functions
   //////////////////////////////////////////////////////////////////////
   namespace
   {
-    double integral_Linfty[] = {2.0,
-                                3.14159265358979323846264338328,
-                                4.18879020478639098461685784437};
-    double integral_W1[]     = {1.0,
-                            1.04719755119659774615421446109,
-                            1.04719755119659774615421446109};
-
-    double integral_Cinfty[] = {1.20690032243787617533623799633,
-                                1.26811216112759608094632335664,
-                                1.1990039070192139033798473858};
+    // Integral of CutOffFunctionLinfty in dimension 1, 2, and 3 when the radius
+    // is one
+    const double integral_Linfty[] = {2.0,
+                                      3.14159265358979323846264338328,
+                                      4.18879020478639098461685784437};
+
+    // Integral of CutOffFunctionW1 in dimension 1, 2, and 3 when the radius
+    // is one
+    const double integral_W1[] = {1.0,
+                                  1.04719755119659774615421446109,
+                                  1.04719755119659774615421446109};
+
+    // Integral of CutOffFunctionCinfty in dimension 1, 2, and 3 when the radius
+    // is one
+    const double integral_Cinfty[] = {1.20690032243787617533623799633,
+                                      1.26811216112759608094632335664,
+                                      1.1990039070192139033798473858};
+
+    // Integral of CutOffFunctionC1 in dimension 1, 2, and 3 when the radius
+    // is one
+    const double integral_C1[] = {1.0,
+                                  0.93417655442731527615578663815,
+                                  0.821155557658032806157358815206};
   } // namespace
 
 
@@ -478,6 +491,115 @@ namespace Functions
   }
 
 
+
+  template <int dim>
+  CutOffFunctionC1<dim>::CutOffFunctionC1(const double       r,
+                                          const Point<dim>   p,
+                                          const unsigned int n_components,
+                                          const unsigned int select,
+                                          bool               integrate_to_one)
+    : CutOffFunctionBase<dim>(r,
+                              p,
+                              n_components,
+                              select,
+                              integrate_to_one,
+                              integral_C1[dim - 1])
+  {}
+
+
+  template <int dim>
+  double
+  CutOffFunctionC1<dim>::value(const Point<dim> & p,
+                               const unsigned int component) const
+  {
+    if (this->selected == CutOffFunctionBase<dim>::no_component ||
+        component == this->selected)
+      {
+        const double d = this->center.distance(p);
+        const double r = this->radius;
+        if (d >= r)
+          return 0.;
+        return .5 * (std::cos(numbers::PI * d / r) + 1) * this->rescaling;
+      }
+    return 0.;
+  }
+
+
+  template <int dim>
+  void
+  CutOffFunctionC1<dim>::value_list(const std::vector<Point<dim>> &points,
+                                    std::vector<double> &          values,
+                                    const unsigned int component) const
+  {
+    Assert(values.size() == points.size(),
+           ExcDimensionMismatch(values.size(), points.size()));
+
+    const double r = this->radius;
+
+    if (this->selected == CutOffFunctionBase<dim>::no_component ||
+        component == this->selected)
+      for (unsigned int i = 0; i < values.size(); ++i)
+        {
+          const double d = this->center.distance(points[i]);
+          if (d >= r)
+            {
+              values[i] = 0.;
+            }
+          else
+            {
+              values[i] =
+                .5 * (std::cos(numbers::PI * d / r) + 1) * this->rescaling;
+            }
+        }
+    else
+      std::fill(values.begin(), values.end(), 0.);
+  }
+
+
+  template <int dim>
+  void
+  CutOffFunctionC1<dim>::vector_value_list(
+    const std::vector<Point<dim>> &points,
+    std::vector<Vector<double>> &  values) const
+  {
+    Assert(values.size() == points.size(),
+           ExcDimensionMismatch(values.size(), points.size()));
+
+    for (unsigned int k = 0; k < values.size(); ++k)
+      {
+        const double d   = this->center.distance(points[k]);
+        const double r   = this->radius;
+        double       val = 0.;
+        if (d < this->radius)
+          {
+            val = .5 * (std::cos(numbers::PI * d / r) + 1) * this->rescaling;
+          }
+
+        if (this->selected == CutOffFunctionBase<dim>::no_component)
+          values[k] = val;
+        else
+          {
+            values[k]                 = 0;
+            values[k](this->selected) = val;
+          }
+      }
+  }
+
+
+
+  template <int dim>
+  Tensor<1, dim>
+  CutOffFunctionC1<dim>::gradient(const Point<dim> &p, const unsigned int) const
+  {
+    const double d = this->center.distance(p);
+    const double r = this->radius;
+    if (d >= r)
+      return Tensor<1, dim>();
+    return (-0.5 * numbers::PI * std::sin(numbers::PI * d / r) / r) *
+           (p - this->center) / d * this->rescaling;
+  }
+
+
   // explicit instantiations
   template class CutOffFunctionBase<1>;
   template class CutOffFunctionBase<2>;
@@ -494,6 +616,14 @@ namespace Functions
   template class CutOffFunctionCinfty<1>;
   template class CutOffFunctionCinfty<2>;
   template class CutOffFunctionCinfty<3>;
+
+  template class CutOffFunctionC1<1>;
+  template class CutOffFunctionC1<2>;
+  template class CutOffFunctionC1<3>;
+
+  template class CutOffFunctionTensorProduct<1>;
+  template class CutOffFunctionTensorProduct<2>;
+  template class CutOffFunctionTensorProduct<3>;
 } // namespace Functions
 
 DEAL_II_NAMESPACE_CLOSE
index c51ac0f797477fd50be107c1c168d38e66733a53..2757d5bcf82a57451cf5cc6a75adaa395b2db768 100644 (file)
@@ -117,4 +117,8 @@ main()
   test<1, Functions::CutOffFunctionCinfty>();
   test<2, Functions::CutOffFunctionCinfty>();
   test<3, Functions::CutOffFunctionCinfty>();
+
+  test<1, Functions::CutOffFunctionC1>();
+  test<2, Functions::CutOffFunctionC1>();
+  test<3, Functions::CutOffFunctionC1>();
 }
index 2554c237b72e0c50617fbbdbf025b0d18514680e..f36ce814d264369dfcc87971e34ba1a442675a6d 100644 (file)
@@ -62,3 +62,24 @@ DEAL::Integrating dealii::Functions::CutOffFunctionCinfty<3> -2,2 cube: 0.999614
 DEAL::Center: 0.500000 0.500000 0.500000
 DEAL::Radius: 0.500000
 DEAL::Integrating dealii::Functions::CutOffFunctionCinfty<3> -2,2 cube: 1.00165
+DEAL::Testing dim = 1
+DEAL::Center: 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionC1<1> -2,2 cube: 1.00000
+DEAL::Center: 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionC1<1> -2,2 cube: 1.00000
+DEAL::Testing dim = 2
+DEAL::Center: 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionC1<2> -2,2 cube: 0.999972
+DEAL::Center: 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionC1<2> -2,2 cube: 0.999764
+DEAL::Testing dim = 3
+DEAL::Center: 0.00000 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionC1<3> -2,2 cube: 0.999478
+DEAL::Center: 0.500000 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionC1<3> -2,2 cube: 1.00089
index 5540e56745dfc35fd7af234d49bfa8d22f0f9516..d2997c14d69a3187f98ba36426646d54fd31749b 100644 (file)
@@ -117,4 +117,8 @@ main()
   test<1, Functions::CutOffFunctionCinfty>();
   test<2, Functions::CutOffFunctionCinfty>();
   test<3, Functions::CutOffFunctionCinfty>();
+
+  test<1, Functions::CutOffFunctionC1>();
+  test<2, Functions::CutOffFunctionC1>();
+  test<3, Functions::CutOffFunctionC1>();
 }
index 2f3fe3da81d82578e8fbf810d988c002e9612490..d6943309948f666dc6bdf07354c154dafe9c80bd 100644 (file)
@@ -62,3 +62,24 @@ DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<3> -2,2 cube: 1
 DEAL::Center: 0.500000 0.500000 0.500000
 DEAL::Radius: 0.500000
 DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<3> -2,2 cube: 1.00779
+DEAL::Testing dim = 1
+DEAL::Center: 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<1> -2,2 cube: 1.00000
+DEAL::Center: 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<1> -2,2 cube: 1.00000
+DEAL::Testing dim = 2
+DEAL::Center: 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<2> -2,2 cube: 1.00000
+DEAL::Center: 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<2> -2,2 cube: 1.00000
+DEAL::Testing dim = 3
+DEAL::Center: 0.00000 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<3> -2,2 cube: 1.00000
+DEAL::Center: 0.500000 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<3> -2,2 cube: 1.00009

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.