From 5c5556a1124ff7311d6431ea2ef89c4f25c552b8 Mon Sep 17 00:00:00 2001 From: angelrca Date: Wed, 1 Jul 2015 10:44:25 +0200 Subject: [PATCH] New class Function::Polynomial New class Function::Polynomial New class Function::Polynomial New class Function::Polynomial --- doc/news/changes.h | 6 ++ include/deal.II/base/function_lib.h | 64 ++++++++++++++++ source/base/function_lib.cc | 110 ++++++++++++++++++++++++++++ tests/base/functions_12.cc | 84 +++++++++++++++++++++ tests/base/functions_12.output | 25 +++++++ 5 files changed, 289 insertions(+) create mode 100644 tests/base/functions_12.cc create mode 100644 tests/base/functions_12.output diff --git a/doc/news/changes.h b/doc/news/changes.h index 8f0c051b50..ad4129e48d 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -385,6 +385,12 @@ inconvenience this causes.
    +
  1. New: Added the class Functions::Polynomial for representation of polynomials. + The new class is derived from the Function class. +
    + (Angel Rodriguez, 2015/07/01) +
    2. +
  2. New: deal.II now supports compilation in C++14 mode, which may be enabled with the CMake option DEAL_II_WITH_CXX14.
    diff --git a/include/deal.II/base/function_lib.h b/include/deal.II/base/function_lib.h index 42ec3ca5b3..4c023a689d 100644 --- a/include/deal.II/base/function_lib.h +++ b/include/deal.II/base/function_lib.h @@ -1254,6 +1254,70 @@ namespace Functions */ const Table 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 to describe the set of coefficients. + * + * @author Ángel Rodríguez, 2015 + */ + template + class Polynomial : public Function + { + 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 &coefficients); + + /** + * Function value at one point. + */ + virtual double value (const Point &p, + const unsigned int component = 0) const; + + + /** + * Function values at multiple points. + */ + virtual void value_list (const std::vector > &points, + std::vector &values, + const unsigned int component = 0) const; + + /** + * Function gradient at one point. + */ + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component = 0) const; + + private: + + /** + * The set of exponents. + */ + const Table<2,double> exponents; + + /** + * The set of coefficients. + */ + const Vector coefficients; + }; + + + } DEAL_II_NAMESPACE_CLOSE diff --git a/source/base/function_lib.cc b/source/base/function_lib.cc index ea7cc38dc2..3553ad153f 100644 --- a/source/base/function_lib.cc +++ b/source/base/function_lib.cc @@ -2483,6 +2483,113 @@ namespace Functions return interpolate (data_values, ix, p_unit); } + /* ---------------------- Polynomial ----------------------- */ + + + + template + Polynomial:: + Polynomial(const Table<2,double> &exponents, + const Vector &coefficients) + : + Function (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 + double + Polynomial::value (const Point &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 + void + Polynomial::value_list (const std::vector > &points, + std::vector &values, + const unsigned int component) const + { + Assert (values.size() == points.size(), + ExcDimensionMismatch(values.size(), points.size())); + + for (unsigned int i=0; i::value (points[i], component); + } + + + + template + Tensor<1,dim> + Polynomial::gradient (const Point &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; 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 index 0000000000..8c19144441 --- /dev/null +++ b/tests/base/functions_12.cc @@ -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 + +template +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 coeffs(n_mon); + for (unsigned int i = 0; i < n_mon; ++i) + coeffs[i] = std::pow(-1,i)*(i+1); + + + CustomFunctions::Polynomial poly(exponents,coeffs); + + Point p; + for (unsigned int d=0; d(); + check<2>(); + check<3>(); + +} + + + diff --git a/tests/base/functions_12.output b/tests/base/functions_12.output new file mode 100644 index 0000000000..da19317833 --- /dev/null +++ b/tests/base/functions_12.output @@ -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:: -- 2.39.5