From 97447bddf9d0d28ae7a83e925ce60aafc47ea5e6 Mon Sep 17 00:00:00 2001
From: Daniel Arndt <d.arndt@math.uni-goettingen.de>
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 <deal.II/base/config.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/polynomial.h>
+#include <deal.II/base/tensor_product_polynomials.h>
+#include <deal.II/base/utilities.h>
+
+#include <vector>
+
+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 <int dim>
+class TensorProductPolynomialsBubbles : public TensorProductPolynomials<dim>
+{
+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. <tt>pols</tt> 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 <class Pol>
+  TensorProductPolynomialsBubbles (const std::vector<Pol> &pols);
+
+  /**
+   * Computes the value and the first and second derivatives of each tensor
+   * product polynomial at <tt>unit_point</tt>.
+   *
+   * 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<dim>            &unit_point,
+                std::vector<double>         &values,
+                std::vector<Tensor<1,dim> > &grads,
+                std::vector<Tensor<2,dim> > &grad_grads) const;
+
+  /**
+   * Computes the value of the <tt>i</tt>th tensor product polynomial at
+   * <tt>unit_point</tt>. Here <tt>i</tt> 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
+   * <tt>values.size()==</tt>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<dim> &p) const;
+
+  /**
+   * Computes the grad of the <tt>i</tt>th tensor product polynomial at
+   * <tt>unit_point</tt>. Here <tt>i</tt> 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
+   * <tt>grads.size()==</tt>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<dim> &p) const;
+
+  /**
+   * Computes the second derivative (grad_grad) of the <tt>i</tt>th tensor
+   * product polynomial at <tt>unit_point</tt>. Here <tt>i</tt> 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
+   * <tt>grad_grads.size()==</tt>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<dim> &p) const;
+
+  /**
+   * Returns the number of tensor product polynomials plus the bubble enrichments.
+   * For <i>n</i> 1d polynomials this is <i>n<sup>dim</sup>+1</i> if the maximum
+   * degree of the polynomials is one and <i>n<sup>dim</sup>+dim</i> otherwise.
+   */
+  unsigned int n () const;
+};
+
+/** @} */
+
+
+/* ---------------- template and inline functions ---------- */
+
+#ifndef DOXYGEN
+
+template <int dim>
+template <class Pol>
+inline
+TensorProductPolynomialsBubbles<dim>::
+TensorProductPolynomialsBubbles(const std::vector<Pol> &pols)
+  :
+  TensorProductPolynomials<dim>(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; i<n_bubbles; ++i)
+    {
+      this->index_map.push_back(i+this->n_tensor_pols);
+      this->index_map_inverse.push_back(i+this->n_tensor_pols);
+    }
+}
+
+
+
+template <int dim>
+inline
+unsigned int
+TensorProductPolynomialsBubbles<dim>::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 <deal.II/base/tensor_product_polynomials_bubbles.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/table.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+
+/* ------------------- TensorProductPolynomialsBubbles -------------- */
+
+
+template <int dim>
+double
+TensorProductPolynomialsBubbles<dim>::compute_value (const unsigned int i,
+                                                     const Point<dim> &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 (i<max_q_indices+n_bubbles, ExcInternalError());
+
+  // treat the regular basis functions
+  if (i<max_q_indices)
+    return this->TensorProductPolynomials<dim>::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<dim; ++j)
+    value*=4*p(j)*(1-p(j));
+  // and multiply with (2x_i-1)^{r-1}
+  for (unsigned int i=0; i<q_degree-1; ++i)
+    value*=2*p(comp)-1;
+  return value;
+}
+
+
+
+template <>
+double
+TensorProductPolynomialsBubbles<0>::compute_value (const unsigned int ,
+                                                   const Point<0> &) const
+{
+  Assert (false, ExcNotImplemented());
+  return 0.;
+}
+
+
+template <int dim>
+Tensor<1,dim>
+TensorProductPolynomialsBubbles<dim>::compute_grad (const unsigned int i,
+                                                    const Point<dim> &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 (i<max_q_indices+n_bubbles, ExcInternalError());
+
+  // treat the regular basis functions
+  if (i<max_q_indices)
+    return this->TensorProductPolynomials<dim>::compute_grad(i,p);
+
+  const unsigned int comp = i - this->n_tensor_pols;
+  Tensor<1,dim> grad;
+
+  for (unsigned int d=0; d<dim ; ++d)
+    {
+      grad[d] = 1.;
+      //compute grad(4*\prod_{i=1}^d (x_i(1-x_i)))(p)
+      for (unsigned j=0; j<dim; ++j)
+        grad[d] *= (d==j ? 4*(1-2*p(j)) : 4*p(j)*(1-p(j)));
+      // and multiply with (2*x_i-1)^{r-1}
+      for (unsigned int i=0; i<q_degree-1; ++i)
+        grad[d]*=2*p(comp)-1;
+    }
+
+  if (q_degree>=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<q_degree-2; ++i)
+        tmp*=2*p(comp)-1;
+      grad[comp]+=tmp;
+    }
+
+  return grad;
+}
+
+
+
+template <int dim>
+Tensor<2,dim>
+TensorProductPolynomialsBubbles<dim>::compute_grad_grad (const unsigned int i,
+                                                         const Point<dim> &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 (i<max_q_indices+n_bubbles, ExcInternalError());
+
+  // treat the regular basis functions
+  if (i<max_q_indices)
+    return this->TensorProductPolynomials<dim>::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<dim; ++c)
+      {
+        v[c][0] = 4*p(c)*(1-p(c));
+        v[c][1] = 4*(1-2*p(c));
+        v[c][2] = -8;
+      }
+
+    double tmp=1.;
+    for (unsigned int i=0; i<q_degree-1; ++i)
+      tmp *= 2*p(comp)-1;
+    v[dim][0] = tmp;
+
+    if (q_degree>=2)
+      {
+        double tmp = 2*(q_degree-1);
+        for (unsigned int i=0; i<q_degree-2; ++i)
+          tmp *= 2*p(comp)-1;
+        v[dim][1] = tmp;
+      }
+    else
+      v[dim][1] = 0.;
+
+    if (q_degree>=3)
+      {
+        double tmp=4*(q_degree-2)*(q_degree-1);
+        for (unsigned int i=0; i<q_degree-3; ++i)
+          tmp *= 2*p(comp)-1;
+        v[dim][2] = tmp;
+      }
+    else
+      v[dim][2] = 0.;
+  }
+
+  //calculate (\partial_j \partial_k \psi) * monomial
+  Tensor<2,dim> grad_grad_1;
+  for (unsigned int d1=0; d1<dim; ++d1)
+    for (unsigned int d2=0; d2<dim; ++d2)
+      {
+        grad_grad_1[d1][d2] = v[dim][0];
+        for (unsigned int x=0; x<dim; ++x)
+          {
+            unsigned int derivative=0;
+            if (d1==x || d2==x)
+              {
+                if (d1==d2)
+                  derivative=2;
+                else
+                  derivative=1;
+              }
+            grad_grad_1[d1][d2] *= v[x][derivative];
+          }
+      }
+
+  //calculate (\partial_j  \psi) *(\partial_k monomial)
+  // and (\partial_k  \psi) *(\partial_j monomial)
+  Tensor<2,dim> grad_grad_2;
+  Tensor<2,dim> grad_grad_3;
+  for (unsigned int d=0; d<dim; ++d)
+    {
+      grad_grad_2[d][comp] = v[dim][1];
+      grad_grad_3[comp][d] = v[dim][1];
+      for (unsigned int x=0; x<dim; ++x)
+        {
+          grad_grad_2[d][comp] *= v[x][d==x];
+          grad_grad_3[comp][d] *= v[x][d==x];
+        }
+    }
+
+  //calculate \psi *(\partial j \partial_k monomial) and sum
+  Tensor<2,dim> grad_grad;
+  double psi_value = 1.;
+  for (unsigned int x=0; x<dim; ++x)
+    psi_value *= v[x][0];
+
+  for (unsigned int d1=0; d1<dim; ++d1)
+    for (unsigned int d2=0; d2<dim; ++d2)
+      grad_grad[d1][d2] = grad_grad_1[d1][d2]
+                          +grad_grad_2[d1][d2]
+                          +grad_grad_3[d1][d2];
+  grad_grad[comp][comp]+=psi_value*v[dim][2];
+
+  return grad_grad;
+}
+
+template <int dim>
+void
+TensorProductPolynomialsBubbles<dim>::
+compute (const Point<dim>            &p,
+         std::vector<double>         &values,
+         std::vector<Tensor<1,dim> > &grads,
+         std::vector<Tensor<2,dim> > &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<dim>::compute(p, values, grads, grad_grads);
+
+  for (unsigned int i=this->n_tensor_pols; i<this->n_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