From 97447bddf9d0d28ae7a83e925ce60aafc47ea5e6 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Wed, 12 Aug 2015 04:47:51 -0500 Subject: [PATCH] Implement a tensor product polynomial space enriched by bubble polynomials --- .../base/tensor_product_polynomials_bubbles.h | 187 ++++++++++++ source/base/CMakeLists.txt | 1 + .../tensor_product_polynomials_bubbles.cc | 273 ++++++++++++++++++ 3 files changed, 461 insertions(+) create mode 100644 include/deal.II/base/tensor_product_polynomials_bubbles.h create mode 100644 source/base/tensor_product_polynomials_bubbles.cc diff --git a/include/deal.II/base/tensor_product_polynomials_bubbles.h b/include/deal.II/base/tensor_product_polynomials_bubbles.h new file mode 100644 index 0000000000..deddfbd08d --- /dev/null +++ b/include/deal.II/base/tensor_product_polynomials_bubbles.h @@ -0,0 +1,187 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2012 - 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. +// +// --------------------------------------------------------------------- + +#ifndef __deal2__tensor_product_polynomials_bubbles_h +#define __deal2__tensor_product_polynomials_bubbles_h + + +#include +#include +#include +#include +#include +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + + +/** + * @addtogroup Polynomials + * @{ + */ + +/** + * Tensor product of given polynomials and bubble functions of form + * $(2*x_j-1)^{degree-1}\prod_{i=0}^{dim-1}(x_i(1-x_i))$. This class inherits most of its + * functionality from TensorProductPolynomials. The bubble enrichments + * are added for the last indices. + * index. + * + * @author Daniel Arndt, 2015 + */ +template +class TensorProductPolynomialsBubbles : public TensorProductPolynomials +{ +public: + /** + * Access to the dimension of + * this object, for checking and + * automatic setting of dimension + * in other classes. + */ + static const unsigned int dimension = dim; + + /** + * Constructor. pols is a vector of objects that should be derived + * or otherwise convertible to one-dimensional polynomial objects. It will + * be copied element by element into a private variable. + */ + template + TensorProductPolynomialsBubbles (const std::vector &pols); + + /** + * Computes the value and the first and second derivatives of each tensor + * product polynomial at unit_point. + * + * The size of the vectors must either be equal 0 or equal n(). In the first + * case, the function will not compute these values. + * + * If you need values or derivatives of all tensor product polynomials then + * use this function, rather than using any of the compute_value(), + * compute_grad() or compute_grad_grad() functions, see below, in a loop + * over all tensor product polynomials. + */ + void compute (const Point &unit_point, + std::vector &values, + std::vector > &grads, + std::vector > &grad_grads) const; + + /** + * Computes the value of the ith tensor product polynomial at + * unit_point. Here i is given in tensor product + * numbering. + * + * Note, that using this function within a loop over all tensor product + * polynomials is not efficient, because then each point value of the + * underlying (one-dimensional) polynomials is (unnecessarily) computed + * several times. Instead use the compute() function with + * values.size()==n() to get the point values of all tensor + * polynomials all at once and in a much more efficient way. + */ + double compute_value (const unsigned int i, + const Point &p) const; + + /** + * Computes the grad of the ith tensor product polynomial at + * unit_point. Here i is given in tensor product + * numbering. + * + * Note, that using this function within a loop over all tensor product + * polynomials is not efficient, because then each derivative value of the + * underlying (one-dimensional) polynomials is (unnecessarily) computed + * several times. Instead use the compute() function, see above, with + * grads.size()==n() to get the point value of all tensor + * polynomials all at once and in a much more efficient way. + */ + Tensor<1,dim> compute_grad (const unsigned int i, + const Point &p) const; + + /** + * Computes the second derivative (grad_grad) of the ith tensor + * product polynomial at unit_point. Here i is given in + * tensor product numbering. + * + * Note, that using this function within a loop over all tensor product + * polynomials is not efficient, because then each derivative value of the + * underlying (one-dimensional) polynomials is (unnecessarily) computed + * several times. Instead use the compute() function, see above, with + * grad_grads.size()==n() to get the point value of all tensor + * polynomials all at once and in a much more efficient way. + */ + Tensor<2,dim> compute_grad_grad (const unsigned int i, + const Point &p) const; + + /** + * Returns the number of tensor product polynomials plus the bubble enrichments. + * For n 1d polynomials this is ndim+1 if the maximum + * degree of the polynomials is one and ndim+dim otherwise. + */ + unsigned int n () const; +}; + +/** @} */ + + +/* ---------------- template and inline functions ---------- */ + +#ifndef DOXYGEN + +template +template +inline +TensorProductPolynomialsBubbles:: +TensorProductPolynomialsBubbles(const std::vector &pols) + : + TensorProductPolynomials(pols) +{ + const unsigned int q_degree = this->polynomials.size()-1; + const unsigned int n_bubbles = ((q_degree<=1)?1:dim); + // append index for renumbering + for (unsigned int i=0; iindex_map.push_back(i+this->n_tensor_pols); + this->index_map_inverse.push_back(i+this->n_tensor_pols); + } +} + + + +template +inline +unsigned int +TensorProductPolynomialsBubbles::n() const +{ + return this->n_tensor_pols+dim; +} + + + +template <> +inline +unsigned int +TensorProductPolynomialsBubbles<0>::n() const +{ + return numbers::invalid_unsigned_int; +} + + +#endif // DOXYGEN +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/base/CMakeLists.txt b/source/base/CMakeLists.txt index f59d76cd10..9102a1e0ec 100644 --- a/source/base/CMakeLists.txt +++ b/source/base/CMakeLists.txt @@ -61,6 +61,7 @@ SET(_src table_handler.cc tensor_function.cc tensor_product_polynomials.cc + tensor_product_polynomials_bubbles.cc tensor_product_polynomials_const.cc thread_management.cc timer.cc diff --git a/source/base/tensor_product_polynomials_bubbles.cc b/source/base/tensor_product_polynomials_bubbles.cc new file mode 100644 index 0000000000..8aa55d796e --- /dev/null +++ b/source/base/tensor_product_polynomials_bubbles.cc @@ -0,0 +1,273 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2012 - 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. +// +// --------------------------------------------------------------------- + + +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + + + +/* ------------------- TensorProductPolynomialsBubbles -------------- */ + + +template +double +TensorProductPolynomialsBubbles::compute_value (const unsigned int i, + const Point &p) const +{ + const unsigned int q_degree = this->polynomials.size()-1; + const unsigned int max_q_indices = this->n_tensor_pols; + const unsigned int n_bubbles = ((q_degree<=1)?1:dim); + (void)n_bubbles; + Assert (iTensorProductPolynomials::compute_value(i,p); + + const unsigned int comp = i - this->n_tensor_pols; + + //compute \prod_{i=1}^d 4*(1-x_i^2)(p) + double value=1.; + for (unsigned int j=0; j +double +TensorProductPolynomialsBubbles<0>::compute_value (const unsigned int , + const Point<0> &) const +{ + Assert (false, ExcNotImplemented()); + return 0.; +} + + +template +Tensor<1,dim> +TensorProductPolynomialsBubbles::compute_grad (const unsigned int i, + const Point &p) const +{ + const unsigned int q_degree = this->polynomials.size()-1; + const unsigned int max_q_indices = this->n_tensor_pols; + const unsigned int n_bubbles = ((q_degree<=1)?1:dim); + (void)n_bubbles; + Assert (iTensorProductPolynomials::compute_grad(i,p); + + const unsigned int comp = i - this->n_tensor_pols; + Tensor<1,dim> grad; + + for (unsigned int d=0; d=2) + { + //add \prod_{i=1}^d 4*(x_i(1-x_i))(p) + double value=1.; + for (unsigned int j=0; j < dim; ++j) + value*=4*p(j)*(1-p(j)); + //and multiply with grad(2*x_i-1)^{r-1} + double tmp=value*2*(q_degree-1); + for (unsigned int i=0; i +Tensor<2,dim> +TensorProductPolynomialsBubbles::compute_grad_grad (const unsigned int i, + const Point &p) const +{ + const unsigned int q_degree = this->polynomials.size()-1; + const unsigned int max_q_indices = this->n_tensor_pols; + const unsigned int n_bubbles = ((q_degree<=1)?1:dim); + (void)n_bubbles; + Assert (iTensorProductPolynomials::compute_grad_grad(i,p); + + const unsigned int comp = i - this->n_tensor_pols; + + double v [dim+1][3]; + { + for (unsigned int c=0; c=2) + { + double tmp = 2*(q_degree-1); + for (unsigned int i=0; i=3) + { + double tmp=4*(q_degree-2)*(q_degree-1); + for (unsigned int i=0; i grad_grad_1; + for (unsigned int d1=0; d1 grad_grad_2; + Tensor<2,dim> grad_grad_3; + for (unsigned int d=0; d grad_grad; + double psi_value = 1.; + for (unsigned int x=0; x +void +TensorProductPolynomialsBubbles:: +compute (const Point &p, + std::vector &values, + std::vector > &grads, + std::vector > &grad_grads) const +{ + const unsigned int q_degree = this->polynomials.size()-1; + const unsigned int max_q_indices = this->n_tensor_pols; + (void) max_q_indices; + const unsigned int n_bubbles = ((q_degree<=1)?1:dim); + Assert (values.size()==max_q_indices+n_bubbles || values.size()==0, + ExcDimensionMismatch2(values.size(), max_q_indices+n_bubbles, 0)); + Assert (grads.size()==max_q_indices+n_bubbles || grads.size()==0, + ExcDimensionMismatch2(grads.size(), max_q_indices+n_bubbles, 0)); + Assert (grad_grads.size()==max_q_indices+n_bubbles || grad_grads.size()==0, + ExcDimensionMismatch2(grad_grads.size(), max_q_indices+n_bubbles, 0)); + + bool do_values = false, do_grads = false, do_grad_grads = false; + if (values.empty() == false) + { + values.resize(this->n_tensor_pols); + do_values = true; + } + if (grads.empty() == false) + { + grads.resize(this->n_tensor_pols); + do_grads = true; + } + if (grad_grads.empty() == false) + { + grad_grads.resize(this->n_tensor_pols); + do_grad_grads = true; + } + + this->TensorProductPolynomials::compute(p, values, grads, grad_grads); + + for (unsigned int i=this->n_tensor_pols; in_tensor_pols+n_bubbles; ++i) + { + if (do_values) + values.push_back(compute_value(i,p)); + if (do_grads) + grads.push_back(compute_grad(i,p)); + if (do_grad_grads) + grad_grads.push_back(compute_grad_grad(i,p)); + } +} + + +/* ------------------- explicit instantiations -------------- */ +template class TensorProductPolynomialsBubbles<1>; +template class TensorProductPolynomialsBubbles<2>; +template class TensorProductPolynomialsBubbles<3>; + +DEAL_II_NAMESPACE_CLOSE -- 2.39.5