]> https://gitweb.dealii.org/ - dealii.git/commitdiff
gradient of InterpolatedTensorProductGridData
authordanshapero <shapero.daniel@gmail.com>
Sun, 2 Aug 2015 23:10:33 +0000 (16:10 -0700)
committerMatthias Maier <tamiko@43-1.org>
Wed, 16 Sep 2015 22:15:07 +0000 (17:15 -0500)
doc/news/changes.h
include/deal.II/base/function_lib.h
source/base/function_lib.cc
tests/base/functions_10.cc

index 6f5c18bc5eb05eba62d2c6b6eed53dc57f2cc134..29a4d4626ea1baf0e7df43675828c2b1754ff05d 100644 (file)
@@ -150,6 +150,11 @@ inconvenience this causes.
 <a name="general"></a>
 <h3>General</h3>
 
+  <li> New: implemented the gradient method for
+  InterpolatedTensorProductGridData
+  <br>
+  (Daniel Shapero, 2015/08/12)
+  </li>
 
 <ol>
   <li> Changed: All doxygen-generated pages now contain a link to the
index d05b71e1b84fe53c4ebbc096d31e2d8275be5958..67d34f3b408da1f35c840481a41b6cadb6bf6e9f 100644 (file)
@@ -1154,6 +1154,22 @@ namespace Functions
     value (const Point<dim> &p,
            const unsigned int component = 0) const;
 
+    /**
+     * Compute the gradient of the function defined by bilinear interpolation
+     * of the given data set.
+     *
+     * @param p The point at which the function gradient is to be evaluated.
+     * @param component The vector component. Since this function is scalar,
+     * only zero is a valid argument here.
+     * @return The value of the gradient of the interpolated function at this
+     * point. If the point lies outside the set of coordinates, the function
+     * is extended by a constant and so its gradient is extended by 0.
+     */
+    virtual
+    Tensor<1, dim>
+    gradient (const Point<dim>    &p,
+              const unsigned int component = 0) const;
+
   private:
     /**
      * The set of coordinate values in each of the coordinate directions.
index e86263e0f6d71ee14340e77587e6f77f3085dc03..ab271f81b603d805b148f2ec00fcc4a6ad44dbfe 100644 (file)
@@ -2340,9 +2340,87 @@ namespace Functions
                 +
                 p_unit[0]*data_values[ix[0]+1][ix[1]+1][ix[2]+1])*p_unit[1]) * p_unit[2]);
     }
-  }
 
 
+    // Interpolate the gradient of a data value from a table where ix
+    // denotes the lower left endpoint of the interval to interpolate
+    // in, p_unit denotes the point in unit coordinates, and dx
+    // denotes the width of the interval in each dimension.
+    Tensor<1,1> gradient_interpolate (const Table<1,double> &data_values,
+                                      const TableIndices<1> &ix,
+                                      const Point<1>        &p_unit,
+                                      const Point<1>        &dx)
+    {
+      (void)p_unit;
+      Tensor<1,1> grad;
+      grad[0] = (data_values[ix[0]+1] - data_values[ix[0]]) / dx[0];
+      return grad;
+    }
+
+
+    Tensor<1,2> gradient_interpolate (const Table<2,double> &data_values,
+                                      const TableIndices<2> &ix,
+                                      const Point<2>        &p_unit,
+                                      const Point<2>        &dx)
+    {
+      Tensor<1,2> grad;
+      double
+        u00 = data_values[ix[0]][ix[1]],
+        u01 = data_values[ix[0]+1][ix[1]],
+        u10 = data_values[ix[0]][ix[1]+1],
+        u11 = data_values[ix[0]+1][ix[1]+1];
+
+      grad[0] = ((1-p_unit[1])*(u01-u00) + p_unit[1]*(u11-u10))/dx[0];
+      grad[1] = ((1-p_unit[0])*(u10-u00) + p_unit[0]*(u11-u01))/dx[1];
+      return grad;
+    }
+
+
+    Tensor<1,3> gradient_interpolate (const Table<3,double> &data_values,
+                                      const TableIndices<3> &ix,
+                                      const Point<3>        &p_unit,
+                                      const Point<3>        &dx)
+    {
+      Tensor<1,3> grad;
+      double
+        u000 = data_values[ix[0]][ix[1]][ix[2]],
+        u001 = data_values[ix[0]+1][ix[1]][ix[2]],
+        u010 = data_values[ix[0]][ix[1]+1][ix[2]],
+        u100 = data_values[ix[0]][ix[1]][ix[2]+1],
+        u011 = data_values[ix[0]+1][ix[1]+1][ix[2]],
+        u101 = data_values[ix[0]+1][ix[1]][ix[2]+1],
+        u110 = data_values[ix[0]][ix[1]+1][ix[2]+1],
+        u111 = data_values[ix[0]+1][ix[1]+1][ix[2]+1];
+
+      grad[0] = ((1-p_unit[2])
+                 *
+                 ((1-p_unit[1])*(u001-u000) + p_unit[1]*(u011-u010))
+                 +
+                 p_unit[2]
+                 *
+                 ((1-p_unit[1])*(u101-u100) + p_unit[1]*(u111-u110))
+                 )/dx[0];
+      grad[1] = ((1-p_unit[2])
+                 *
+                 ((1-p_unit[0])*(u010-u000) + p_unit[0]*(u011-u001))
+                 +
+                 p_unit[2]
+                 *
+                 ((1-p_unit[0])*(u110-u100) + p_unit[0]*(u111-u101))
+                 )/dx[1];
+      grad[2] = ((1-p_unit[1])
+                 *
+                 ((1-p_unit[0])*(u100-u000) + p_unit[0]*(u101-u001))
+                 +
+                 p_unit[1]
+                 *
+                 ((1-p_unit[0])*(u110-u010) + p_unit[0]*(u111-u011))
+                 )/dx[2];
+
+      return grad;
+    }
+  }
+
   template <int dim>
   InterpolatedTensorProductGridData<dim>::
   InterpolatedTensorProductGridData(const std_cxx11::array<std::vector<double>,dim> &coordinate_values,
@@ -2415,6 +2493,45 @@ namespace Functions
 
 
 
+  template <int dim>
+  Tensor<1, dim>
+  InterpolatedTensorProductGridData<dim>::gradient(const Point<dim> &p,
+                                                   const unsigned int component) const
+  {
+    (void)component;
+    Assert (component == 0,
+            ExcMessage ("This is a scalar function object, the component can only be zero."));
+
+    // find out where this data point lies
+    TableIndices<dim> ix;
+    for (unsigned int d=0; d<dim; ++d)
+      {
+        ix[d] = (std::lower_bound (coordinate_values[d].begin(),
+                                   coordinate_values[d].end(),
+                                   p[d])
+                 - coordinate_values[d].begin());
+
+        if (ix[d] == coordinate_values[d].size())
+          ix[d] = coordinate_values[d].size()-2;
+        else if (ix[d] > 0)
+          --ix[d];
+      }
+
+    Point<dim> dx;
+    for (unsigned int d=0; d<dim; ++d)
+      dx[d] = coordinate_values[d][ix[d]+1]-coordinate_values[d][ix[d]];
+
+    Point<dim> p_unit;
+    for (unsigned int d=0; d<dim; ++d)
+      p_unit[d] = std::max(std::min((p[d]-coordinate_values[d][ix[d]]) / dx[d],
+                                    1.),
+                           0.0);
+
+    return gradient_interpolate(data_values, ix, p_unit, dx);
+  }
+
+
+
   template <int dim>
   InterpolatedUniformGridData<dim>::
   InterpolatedUniformGridData(const std_cxx11::array<std::pair<double,double>,dim> &interval_endpoints,
index 9ffb6d3be6682f0e65c1ab7105db3cbaaadaee81..7cb1ac0e472c528a024a41e3298a73d9955c2e6e 100644 (file)
@@ -70,7 +70,7 @@ void check ()
   Functions::InterpolatedTensorProductGridData<dim> f(coordinates, data);
 
   // now choose a number of randomly chosen points inside the box and
-  // verify that the functions returned are correct
+  // verify that the function values and gradients returned are correct
   for (unsigned int i=0; i<10; ++i)
     {
       Point<dim> p;
@@ -85,7 +85,18 @@ void check ()
 
       AssertThrow (std::fabs (exact_value - f.value(p)) < 1e-12,
                    ExcInternalError());
-    }
+
+      Tensor<1,dim> exact_gradient;
+      for (unsigned int d=0; d<dim; ++d)
+        {
+          exact_gradient[d] = 1.0;
+          for (unsigned int k=0; k<dim; ++k)
+            exact_gradient[d] *= (k == d) ? 1.0 : p[k];
+        }
+
+      AssertThrow ((exact_gradient - f.gradient(p)).norm() < 1e-12,
+                   ExcInternalError());
+   }
 
   // now also verify that it computes values outside the box correctly, as
   // documented

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.