]> https://gitweb.dealii.org/ - dealii.git/commitdiff
New class Function::Polynomial 1067/head
authorangelrca <angel.rcarballo@udc.es>
Wed, 1 Jul 2015 08:44:25 +0000 (10:44 +0200)
committerangelrca <angel.rcarballo@udc.es>
Thu, 2 Jul 2015 07:44:09 +0000 (09:44 +0200)
New class Function::Polynomial

New class Function::Polynomial

New class Function::Polynomial

doc/news/changes.h
include/deal.II/base/function_lib.h
source/base/function_lib.cc
tests/base/functions_12.cc [new file with mode: 0644]
tests/base/functions_12.output [new file with mode: 0644]

index 8f0c051b501f212fe3e379165616d4eefb32ba8a..ad4129e48d4d3b7945198b08e98bb13f1cad06b4 100644 (file)
@@ -385,6 +385,12 @@ inconvenience this causes.
 
 
 <ol>
+  <li> New: Added the class Functions::Polynomial for representation of polynomials. 
+  The new class is derived from the Function class.
+  <br>
+  (Angel Rodriguez, 2015/07/01) 
+  </li>
+
   <li> New: deal.II now supports compilation in C++14 mode, which may be
   enabled with the CMake option <code>DEAL_II_WITH_CXX14</code>.
   <br>
index 42ec3ca5b361af8908d36e12752f5a85ed3f6458..4c023a689d0e77e61dc26ba37626a97ad56fad7b 100644 (file)
@@ -1254,6 +1254,70 @@ namespace Functions
      */
     const Table<dim,double>                     data_values;
   };
+
+
+  /**
+   * A class that represents a function object for a polynomial. A polynomial
+   * is composed by the summation of multiple monomials. If the polynomial has
+   * n monomials and the dimension is equal to dim, the polynomial can be written as
+   * $\sum_{i=1}^{n} a_{i}(\prod_{d=1}^{dim} x_{d}^{\alpha_{i,d}})$, where
+   * $a_{i}$ are the coefficients of the monomials and $\alpha_{i,d}$ are their exponents.
+   * The class's constructor takes a Table<2,double> to describe the set of
+   * exponents and a Vector<double> to describe the set of coefficients.
+   *
+   *  @author Ángel Rodríguez, 2015
+   */
+  template <int dim>
+  class Polynomial : public Function<dim>
+  {
+  public:
+    /**
+     * Constructor. The coefficients and the exponents of the polynomial are
+     * passed as arguments. The Table<2, double> exponents has a number of rows
+     * equal to the number of monomials of the polynomial and a number of columns
+     * equal to dim. The i-th row of the exponents table contains the
+     * ${\alpha_{i,d}}$ exponents of the i-th monomial
+     * $a_{i}\prod_{d=1}^{dim} x_{d}^{\alpha_{i,d}}$. The i-th element of the coefficients
+     * vector contains the coefficient $a_{i}$ for the i-th monomial.
+     */
+    Polynomial (const Table<2,double> &exponents,
+                const Vector<double>   &coefficients);
+
+    /**
+     * 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 std::vector<Point<dim> > &points,
+                             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:
+
+    /**
+     * The set of exponents.
+     */
+    const Table<2,double> exponents;
+
+    /**
+     * The set of coefficients.
+     */
+    const Vector<double> coefficients;
+  };
+
+
+
 }
 DEAL_II_NAMESPACE_CLOSE
 
index ea7cc38dc26466109c12024529a5ef49c26c9765..3553ad153f5abf2729b1b7506203e5b1682e31c6 100644 (file)
@@ -2483,6 +2483,113 @@ namespace Functions
     return interpolate (data_values, ix, p_unit);
   }
 
+  /* ---------------------- Polynomial ----------------------- */
+
+
+
+  template <int dim>
+  Polynomial<dim>::
+  Polynomial(const Table<2,double> &exponents,
+             const Vector<double>   &coefficients)
+    :
+    Function<dim> (1),
+    exponents (exponents),
+    coefficients(coefficients)
+  {
+    Assert(exponents.n_rows() == coefficients.size(),
+           ExcDimensionMismatch(exponents.n_rows(), coefficients.size()));
+    Assert(exponents.n_cols() == dim,
+           ExcDimensionMismatch(exponents.n_cols(), dim));
+  }
+
+
+
+  template <int dim>
+  double
+  Polynomial<dim>::value (const Point<dim> &p,
+                          const unsigned int component) const
+  {
+    (void)component;
+    Assert (component==0, ExcIndexRange(component,0,1));
+
+    double prod;
+    double sum = 0;
+    for (unsigned int monom = 0; monom < exponents.n_rows(); ++monom)
+      {
+        prod = 1;
+        for (unsigned int s=0; s< dim; ++s)
+          {
+            if (p[s] < 0)
+              Assert(std::floor(exponents[monom][s]) == exponents[monom][s],
+                     ExcMessage("Exponentiation of a negative base number with "
+                                "a real exponent can't be performed."));
+            prod *= std::pow(p[s], exponents[monom][s]);
+          }
+        sum += coefficients[monom]*prod;
+      }
+    return sum;
+  }
+
+
+
+  template <int dim>
+  void
+  Polynomial<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()));
+
+    for (unsigned int i=0; i<points.size(); ++i)
+      values[i] = Polynomial<dim>::value (points[i], component);
+  }
+
+
+
+  template <int dim>
+  Tensor<1,dim>
+  Polynomial<dim>::gradient (const Point<dim> &p,
+                             const unsigned int component) const
+  {
+    (void)component;
+    Assert (component==0, ExcIndexRange(component,0,1));
+
+    Tensor<1,dim> r;
+
+    for (unsigned int d=0; d<dim; ++d)
+      {
+        double sum = 0;
+
+        for (unsigned int monom = 0; monom < exponents.n_rows(); ++monom)
+          {
+            double prod = 1;
+            for (unsigned int s=0; s < dim; ++s)
+              {
+                if ((s==d)&&(exponents[monom][s] == 0)&&(p[s] == 0))
+                  {
+                    prod = 0;
+                    break;
+                  }
+                else
+                  {
+                    if (p[s] < 0)
+                      Assert(std::floor(exponents[monom][s]) == exponents[monom][s],
+                             ExcMessage("Exponentiation of a negative base number with "
+                                        "a real exponent can't be performed."));
+                    prod *= (s==d
+                             ?
+                             exponents[monom][s] * std::pow(p[s], exponents[monom][s]-1)
+                             :
+                             std::pow(p[s], exponents[monom][s]));
+                  }
+              }
+            sum += coefficients[monom]*prod;
+          }
+        r[d] = sum;
+      }
+    return r;
+  }
 
 
 // explicit instantiations
@@ -2532,6 +2639,9 @@ namespace Functions
   template class InterpolatedUniformGridData<1>;
   template class InterpolatedUniformGridData<2>;
   template class InterpolatedUniformGridData<3>;
+  template class Polynomial<1>;
+  template class Polynomial<2>;
+  template class Polynomial<3>;
 }
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/tests/base/functions_12.cc b/tests/base/functions_12.cc
new file mode 100644 (file)
index 0000000..8c19144
--- /dev/null
@@ -0,0 +1,84 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2007 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Test Polynomial
+
+#include "../tests.h"
+#include <deal.II/base/function_lib.h>
+
+template <int dim>
+void check()
+{
+  unsigned int n_mon = 3;
+
+  Table<2,double> exponents(n_mon,dim);
+
+  for (unsigned int i = 0; i < n_mon; ++i)
+    for (unsigned int d = 0; d < dim; ++d)
+      exponents[i][d] = i + d;
+
+  Vector<double> coeffs(n_mon);
+  for (unsigned int i = 0; i < n_mon; ++i)
+    coeffs[i] = std::pow(-1,i)*(i+1);
+
+
+  CustomFunctions::Polynomial<dim> poly(exponents,coeffs);
+
+  Point<dim> p;
+  for (unsigned int d=0; d<dim; ++d)
+    p[d] = d;
+
+  deallog << dim << "-D check" << std::endl;
+  deallog << "Polynomial: ";
+  for (unsigned int i = 0; i < n_mon; ++i)
+    {
+      deallog << coeffs[i];
+      for (unsigned int d = 0; d < dim; ++d)
+        deallog << " x" << d << "^" << exponents[i][d];
+      if (i < n_mon-1)
+        deallog << " + ";
+    }
+  deallog << std::endl;
+  deallog << "Point: ";
+  for (unsigned int d = 0; d < dim; ++d)
+    deallog << p[d] << " ";
+  deallog << std::endl;
+
+  deallog << "Value: " << poly.value(p) << std::endl;
+  deallog << " values checked" << std::endl;
+
+  deallog << "Gradient: " << poly.gradient(p) << std::endl;
+  deallog << " gradients checked" << std::endl;
+  deallog << std::endl;
+
+}
+
+int main()
+{
+  std::string logname = "output";
+  std::ofstream logfile(logname.c_str());
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  check<1>();
+  check<2>();
+  check<3>();
+
+}
+
+
+
diff --git a/tests/base/functions_12.output b/tests/base/functions_12.output
new file mode 100644 (file)
index 0000000..da19317
--- /dev/null
@@ -0,0 +1,25 @@
+
+DEAL::1-D check
+DEAL::Polynomial: 1.00000 x0^0 + -2.00000 x0^1.00000 + 3.00000 x0^2.00000
+DEAL::Point: 0 
+DEAL::Value: 1.00000
+DEAL:: values checked
+DEAL::Gradient: -2.00000
+DEAL:: gradients checked
+DEAL::
+DEAL::2-D check
+DEAL::Polynomial: 1.00000 x0^0 x1^1.00000 + -2.00000 x0^1.00000 x1^2.00000 + 3.00000 x0^2.00000 x1^3.00000
+DEAL::Point: 0 1.00000 
+DEAL::Value: 1.00000
+DEAL:: values checked
+DEAL::Gradient: -2.00000 1.00000
+DEAL:: gradients checked
+DEAL::
+DEAL::3-D check
+DEAL::Polynomial: 1.00000 x0^0 x1^1.00000 x2^2.00000 + -2.00000 x0^1.00000 x1^2.00000 x2^3.00000 + 3.00000 x0^2.00000 x1^3.00000 x2^4.00000
+DEAL::Point: 0 1.00000 2.00000 
+DEAL::Value: 4.00000
+DEAL:: values checked
+DEAL::Gradient: -16.0000 4.00000 4.00000
+DEAL:: gradients checked
+DEAL::

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.