From: guido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Sun, 4 Jan 2004 19:00:08 +0000 (+0000)
Subject: new BDM polynomial space
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=52af1de61432bd8243ea3e67aff3ea6be5914555;p=dealii-svn.git

new BDM polynomial space


git-svn-id: https://svn.dealii.org/trunk@8275 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/base/include/base/polynomial_space.h b/deal.II/base/include/base/polynomial_space.h
index 2f06766044..a361da5637 100644
--- a/deal.II/base/include/base/polynomial_space.h
+++ b/deal.II/base/include/base/polynomial_space.h
@@ -132,6 +132,18 @@ class PolynomialSpace
 				      */
     unsigned int n () const;
 
+				     /**
+				      * Degree of the space. This is
+				      * by definition the number of
+				      * polynomials given to the
+				      * constructor, NOT the maximal
+				      * degree of a polynomial in this
+				      * vector. The latter value is
+				      * never checked and therefore
+				      * left to the application.
+				      */
+    unsigned int degree () const;
+    
 				     /**
 				      * Exception.
 				      */
@@ -202,5 +214,24 @@ PolynomialSpace (const std::vector<Pol> &pols)
 {}
 
 
+template<int dim>
+inline
+unsigned int
+PolynomialSpace<dim>::n() const
+{
+  return n_pols;
+}
+
+  
+
+template<int dim>
+inline
+unsigned int
+PolynomialSpace<dim>::degree() const
+{
+  return polynomials.size();
+}
+
+  
 
 #endif
diff --git a/deal.II/base/include/base/polynomials_bdm.h b/deal.II/base/include/base/polynomials_bdm.h
new file mode 100644
index 0000000000..1d5381a84b
--- /dev/null
+++ b/deal.II/base/include/base/polynomials_bdm.h
@@ -0,0 +1,246 @@
+//-------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000, 2001, 2002, 2003 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-------------------------------------------------------------------
+#ifndef __deal2__polynomials_BDM_h
+#define __deal2__polynomials_BDM_h
+
+
+#include <base/config.h>
+#include <base/exceptions.h>
+#include <base/tensor.h>
+#include <base/point.h>
+#include <base/polynomial.h>
+#include <base/polynomial_space.h>
+#include <base/table.h>
+
+#include <vector>
+
+
+/**
+ * @brief The set of BDM polynomials on tensor product cells
+ *
+ * This class implements the <I>H<SUB>div</SUB></I>-conforming,
+ * vector-valued Brezzi-Douglas-Marini polynomials as described in the
+ * book by Brezzi and Fortin.
+ *
+ * Right now, they are implemented in two dimensions only. There, they
+ * consist of the complete polynomial space of order $k$ plus two
+ * additional vectors.
+ *
+ * @author Guido Kanschat, 2003
+ */
+template <int dim>
+class PolynomialsBDM
+{
+  public:
+				     /**
+				      * Constructor. Creates all basis
+				      * functions for BDM polynomials
+				      * of given degree.
+				      *
+				      * Remark that the degree of a
+				      * BDM space is the degree of the
+				      * largest complete polynomial
+				      * space embedded.
+				      *
+				      * @arg k: the degree of the
+				      * BDM-space
+				      */
+    PolynomialsBDM (const unsigned int k);
+
+				     /**
+				      * Computes the value and the
+				      * first and second derivatives
+				      * of each BDM
+				      * polynomial at @p unit_point.
+				      *
+				      * The size of the vectors must
+				      * either be zero or equal
+				      * <tt>n()</tt>.  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
+				      * @p{compute_value},
+				      * @p{compute_grad} or
+				      * @p{compute_grad_grad}
+				      * functions, see below, in a
+				      * loop over all tensor product
+				      * polynomials.
+				      */
+    void compute (const Point<dim>            &unit_point,
+                  std::vector<Tensor<1,dim> > &values,
+                  std::vector<Tensor<2,dim> > &grads,
+                  std::vector<Tensor<3,dim> > &grad_grads) const;
+    
+				     /**
+				      * Computes the value of the
+				      * @p{i}th BDM
+				      * polynomial at
+				      * @p{unit_point}.
+				      *
+				      * 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 @p{compute} function, see
+				      * above, with
+				      * @p{values.size()==n_tensor_pols}
+				      * to get the point values of all
+				      * tensor polynomials all at once
+				      * and in a much more efficient
+				      * way.
+				      */
+    Tensor<1,dim> compute_value (const unsigned int i,
+				 const Point<dim> &p) const;
+
+				     /**
+				      * Computes the grad of the
+				      * @p{i}th tensor product
+				      * polynomial at
+				      * @p{unit_point}. Here @p{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 @p{compute}
+				      * function, see above, with
+				      * @p{grads.size()==n_tensor_pols}
+				      * to get the point value of all
+				      * tensor polynomials all at once
+				      * and in a much more efficient
+				      * way.
+				      */
+    Tensor<2,dim> compute_grad (const unsigned int i,
+				const Point<dim> &p) const;
+
+				     /**
+				      * Computes the second
+				      * derivative (grad_grad) of the
+				      * @p{i}th tensor product
+				      * polynomial at
+				      * @p{unit_point}. Here @p{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 @p{compute}
+				      * function, see above, with
+				      * @p{grad_grads.size()==n_tensor_pols}
+				      * to get the point value of all
+				      * tensor polynomials all at once
+				      * and in a much more efficient
+				      * way.
+				      */
+    Tensor<3,dim> compute_grad_grad (const unsigned int i,
+                                     const Point<dim> &p) const;
+
+				     /**
+				      * Compute the matrix that has as
+				      * its entry
+				      * <i>a<sub>ij</sub></i> the node
+				      * functional <i>i</i> evaluated
+				      * for basis function
+				      * <i>j</i>. The node functionals
+				      * are the standard BDM
+				      * interpolation operators.
+				      *
+				      * The inverse of this matrix can
+				      * be used to interpolate node
+				      * values to BDM polynomials.
+				      */
+    void compute_node_matrix (Table<2,double>&) const;
+    
+				     /**
+				      * Returns the number of BDM polynomials.
+				      */
+    unsigned int n () const;
+
+				     /**
+				      * Exception.
+				      */
+    DeclException3 (ExcDimensionMismatch2,
+		    int, int, int,
+		    << "Dimension " << arg1 << " not equal to " << arg2 << " nor to " << arg3);
+
+	    
+  private:
+				     /**
+				      * An object representing the
+				      * polynomial space used
+				      * here. The constructor fills
+				      * this with the monomial basis.
+				      */
+    const PolynomialSpace<dim> polynomial_space;
+
+				     /**
+				      * Storage for monomials
+				      */
+    std::vector<Polynomials::Polynomial<double> > monomials;
+    
+				     /**
+				      * Storage for derivatives of monomials
+				      */
+    std::vector<Polynomials::Polynomial<double> > monomial_derivatives;
+    
+				     /**
+				      * Number of BDM
+				      * polynomials.
+				      */
+    unsigned int n_pols;
+
+				     /**
+				      * Auxiliary memory.
+				      */
+    mutable std::vector<double> p_values;
+				     /**
+				      * Auxiliary memory.
+				      */
+    mutable std::vector<Tensor<1,dim> > p_grads;
+				     /**
+				      * Auxiliary memory.
+				      */
+    mutable std::vector<Tensor<2,dim> > p_grad_grads;
+};
+
+
+
+template <int dim>
+inline unsigned int
+PolynomialsBDM<dim>::n() const
+{
+  return n_pols;
+}
+
+#endif
diff --git a/deal.II/base/source/polynomial_space.cc b/deal.II/base/source/polynomial_space.cc
index 0fdfb58753..aefd6e5ff3 100644
--- a/deal.II/base/source/polynomial_space.cc
+++ b/deal.II/base/source/polynomial_space.cc
@@ -287,14 +287,6 @@ PolynomialSpace<dim>::compute (const Point<dim>            &p,
 }
 
 
-template<int dim>
-unsigned int
-PolynomialSpace<dim>::n() const
-{
-  return n_pols;
-}
-
-  
 template class PolynomialSpace<1>;
 template class PolynomialSpace<2>;
 template class PolynomialSpace<3>;
diff --git a/deal.II/base/source/polynomials_bdm.cc b/deal.II/base/source/polynomials_bdm.cc
new file mode 100644
index 0000000000..8c5639df81
--- /dev/null
+++ b/deal.II/base/source/polynomials_bdm.cc
@@ -0,0 +1,166 @@
+//-------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000, 2001, 2002, 2003 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-------------------------------------------------------------------
+
+#include <base/polynomials_bdm.h>
+#include <base/quadrature_lib.h>
+#include <iostream>
+using namespace std;
+using namespace Polynomials;
+
+
+template <int dim>
+PolynomialsBDM<dim>::PolynomialsBDM (const unsigned int k)
+		:
+		polynomial_space (Polynomials::Monomial<double>::generate_complete_basis(k)),
+		monomials(1),
+		monomial_derivatives(1),
+		n_pols(dim * polynomial_space.n()+2),
+		p_values(polynomial_space.n()),
+		p_grads(polynomial_space.n()),
+		p_grad_grads(polynomial_space.n())
+{
+  Assert (dim == 2, ExcNotImplemented());
+  monomials[0] = Monomial<double> (k+1);
+  for (unsigned int i=0;i<monomials.size();++i)
+    monomial_derivatives[i] = monomials[i].derivative();
+}
+
+
+
+template <int dim>
+void
+PolynomialsBDM<dim>::compute (const Point<dim>            &unit_point,
+			      std::vector<Tensor<1,dim> > &values,
+			      std::vector<Tensor<2,dim> > &grads,
+			      std::vector<Tensor<3,dim> > &grad_grads) const
+{
+  Assert(values.size()==n_pols || values.size()==0,
+	 ExcDimensionMismatch2(values.size(), n_pols, 0));
+  Assert(grads.size()==n_pols|| grads.size()==0,
+	 ExcDimensionMismatch2(grads.size(), n_pols, 0));
+  Assert(grad_grads.size()==n_pols|| grad_grads.size()==0,
+	 ExcDimensionMismatch2(grad_grads.size(), n_pols, 0));
+
+  const unsigned int n_sub = polynomial_space.n();
+  p_values.resize((values.size() == 0) ? 0 : n_sub);
+  p_grads.resize((grads.size() == 0) ? 0 : n_sub);
+  p_grad_grads.resize((grad_grads.size() == 0) ? 0 : n_sub);
+
+				   // Compute values of complete space
+				   // and insert into tensors.  Result
+				   // will have first all polynomials
+				   // in the x-component, then y and
+				   // z.
+  polynomial_space.compute (unit_point, p_values, p_grads, p_grad_grads);
+
+  std::fill(values.begin(), values.end(), Tensor<1,dim>());
+  for (unsigned int i=0;i<p_values.size();++i)
+    {
+      for (unsigned int j=0;j<dim;++j)
+	{
+	  values[i+j*n_sub][j] = p_values[i];
+	  std::cerr << i+j*n_sub << ' ' << j << ' ' << p_values[i] << std::endl;
+	}
+      
+    }
+
+				   // Let's hope this is not the transpose
+  std::fill(grads.begin(), grads.end(), Tensor<2,dim>());
+  for (unsigned int i=0;i<p_grads.size();++i)
+    {
+      for (unsigned int j=0;j<dim;++j)
+	grads[i+j*n_sub][j] = p_grads[i];
+    }
+
+				   // Let's hope this is not the transpose
+  std::fill(grad_grads.begin(), grad_grads.end(), Tensor<3,dim>());
+  for (unsigned int i=0;i<p_grad_grads.size();++i)
+    {
+      for (unsigned int j=0;j<dim;++j)
+	grad_grads[i+j*n_sub][j] = p_grad_grads[i];
+    }
+  
+  const unsigned int start = dim*n_sub;
+  if (values.size() != 0)
+    {
+      values[start][0] = monomials[0].value (unit_point(0));
+      values[start][1] = -unit_point(1)
+			 * monomial_derivatives[0].value (unit_point(0));
+      values[start+1][0] = -unit_point(0)
+			 * monomial_derivatives[0].value (unit_point(1));
+      values[start+1][1] = monomials[0].value (unit_point(1));
+    }
+  if (grads.size() != 0)
+    {
+      Assert(false,ExcNotImplemented());
+    }
+  if (grad_grads.size() != 0)
+    {
+      Assert(false,ExcNotImplemented());
+    }
+}
+
+
+template <int dim>
+void
+PolynomialsBDM<dim>::compute_node_matrix (Table<2,double>& A) const
+{
+  std::vector<Polynomial<double> > legendre(2);
+  for (unsigned int i=0;i<legendre.size();++i)
+    legendre[i] = Legendre(i);
+
+  QGauss<1> qface(polynomial_space.degree());
+
+  Table<2,double> integrals (n(), n());
+
+  std::vector<Tensor<1,dim> > values(n());
+  std::vector<Tensor<2,dim> > grads;
+  std::vector<Tensor<3,dim> > grad_grads;
+  values.resize(n());
+
+  for (unsigned int face=0;face<2*dim;++face)
+    for (unsigned int k=0;k<qface.n_quadrature_points;++k)
+      {
+	const double w = qface.weight(k);
+	const double x = qface.point(k)(0);
+	Point<dim> p;
+	switch (face)
+	  {
+	    case 2:
+	      p(1) = 1.;
+	    case 0:
+	      p(0) = x;
+	      break;
+	    case 1:
+	      p(0) = 1.;
+	    case 3:
+	      p(1) = x;
+	      break;	      
+	  }
+	std::cerr << p << std::endl;
+	
+	compute (p, values, grads, grad_grads);
+	for (unsigned int i=0;i<n();++i)
+	  {
+	    for (unsigned int j=0;j<legendre.size();++j)
+	      A(2*face+j,i) += w * values[i][1-face%2] * legendre[j].value(x);
+	  }
+      }
+				   // Volume integrals are missing
+  Assert (polynomial_space.degree() < 2,
+	  ExcNotImplemented());
+}
+
+
+template class PolynomialsBDM<2>;
+