<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>
*/
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
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
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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>();
+
+}
+
+
+
--- /dev/null
+
+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::
projects / dealii.git / commitdiff