]> https://gitweb.dealii.org/ - dealii.git/commitdiff
new BDM polynomial space
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Sun, 4 Jan 2004 19:00:08 +0000 (19:00 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Sun, 4 Jan 2004 19:00:08 +0000 (19:00 +0000)
git-svn-id: https://svn.dealii.org/trunk@8275 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/polynomial_space.h
deal.II/base/include/base/polynomials_bdm.h [new file with mode: 0644]
deal.II/base/source/polynomial_space.cc
deal.II/base/source/polynomials_bdm.cc [new file with mode: 0644]

index 2f06766044de1dba0d2615fd552329bf6f99a7ac..a361da56374fdd0e32239ed54187259f7a153dc1 100644 (file)
@@ -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 (file)
index 0000000..1d5381a
--- /dev/null
@@ -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
index 0fdfb58753c56882b503a9ba5ea6ab6921f88998..aefd6e5ff34b83e39f4e852ea0b15464b2119cb8 100644 (file)
@@ -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 (file)
index 0000000..8c5639d
--- /dev/null
@@ -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>;
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.