From 52af1de61432bd8243ea3e67aff3ea6be5914555 Mon Sep 17 00:00:00 2001 From: guido Date: Sun, 4 Jan 2004 19:00:08 +0000 Subject: [PATCH] new BDM polynomial space git-svn-id: https://svn.dealii.org/trunk@8275 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/polynomial_space.h | 31 +++ deal.II/base/include/base/polynomials_bdm.h | 246 +++++++++++++++++++ deal.II/base/source/polynomial_space.cc | 8 - deal.II/base/source/polynomials_bdm.cc | 166 +++++++++++++ 4 files changed, 443 insertions(+), 8 deletions(-) create mode 100644 deal.II/base/include/base/polynomials_bdm.h create mode 100644 deal.II/base/source/polynomials_bdm.cc 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 &pols) {} +template +inline +unsigned int +PolynomialSpace::n() const +{ + return n_pols; +} + + + +template +inline +unsigned int +PolynomialSpace::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 +#include +#include +#include +#include +#include +#include + +#include + + +/** + * @brief The set of BDM polynomials on tensor product cells + * + * This class implements the Hdiv-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 +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 + * 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 + * @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 &unit_point, + std::vector > &values, + std::vector > &grads, + std::vector > &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 &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 &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 &p) const; + + /** + * Compute the matrix that has as + * its entry + * aij the node + * functional i evaluated + * for basis function + * j. 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 polynomial_space; + + /** + * Storage for monomials + */ + std::vector > monomials; + + /** + * Storage for derivatives of monomials + */ + std::vector > monomial_derivatives; + + /** + * Number of BDM + * polynomials. + */ + unsigned int n_pols; + + /** + * Auxiliary memory. + */ + mutable std::vector p_values; + /** + * Auxiliary memory. + */ + mutable std::vector > p_grads; + /** + * Auxiliary memory. + */ + mutable std::vector > p_grad_grads; +}; + + + +template +inline unsigned int +PolynomialsBDM::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::compute (const Point &p, } -template -unsigned int -PolynomialSpace::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 +#include +#include +using namespace std; +using namespace Polynomials; + + +template +PolynomialsBDM::PolynomialsBDM (const unsigned int k) + : + polynomial_space (Polynomials::Monomial::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 (k+1); + for (unsigned int i=0;i +void +PolynomialsBDM::compute (const Point &unit_point, + std::vector > &values, + std::vector > &grads, + std::vector > &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()); + for (unsigned int i=0;i()); + for (unsigned int i=0;i +void +PolynomialsBDM::compute_node_matrix (Table<2,double>& A) const +{ + std::vector > legendre(2); + for (unsigned int i=0;i qface(polynomial_space.degree()); + + Table<2,double> integrals (n(), n()); + + std::vector > values(n()); + std::vector > grads; + std::vector > grad_grads; + values.resize(n()); + + for (unsigned int face=0;face<2*dim;++face) + for (unsigned int k=0;k 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; + -- 2.39.5