From d451278082486a68aad7949d6cdec9cff6c1d62f Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Sun, 10 May 2009 21:34:46 +0000 Subject: [PATCH] PolynomialsBDM in 3D git-svn-id: https://svn.dealii.org/trunk@18829 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/polynomials_bdm.h | 20 ++- deal.II/base/source/polynomials_bdm.cc | 175 ++++++++++++++++---- 2 files changed, 153 insertions(+), 42 deletions(-) diff --git a/deal.II/base/include/base/polynomials_bdm.h b/deal.II/base/include/base/polynomials_bdm.h index e808736c25..89ed5a7124 100644 --- a/deal.II/base/include/base/polynomials_bdm.h +++ b/deal.II/base/include/base/polynomials_bdm.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2004, 2005, 2006 by the deal.II authors +// Copyright (C) 2004, 2005, 2006, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -43,19 +43,19 @@ DEAL_II_NAMESPACE_OPEN * vector valued polynomials: * *
- *
In 2D: + *
In 2D: *
The 2D-curl of the functions xk+1y * and xyk+1. - *
In 3D: + *
In 3D: *
For any i=0,...,k the curls of * (0,0,xyi+1zk-i), - * (0,xk-iyzi+1,0) and - * (xi+1yk-iz,0,0) + * (xk-iyzi+1,0,0) and + * (0,xi+1yk-iz,0) *
* - * Right now, they are implemented in two dimensions only. + * @todo Second derivatives in 3D are missing. * - * @author Guido Kanschat, 2003, 2005 + * @author Guido Kanschat, 2003, 2005, 2009 */ template class PolynomialsBDM @@ -137,7 +137,11 @@ class PolynomialsBDM const PolynomialSpace polynomial_space; /** - * Storage for monomials + * Storage for monomials. In 2D, + * this is just the polynomial of + * order k. In 3D, we + * need all polynomials from + * degree zero to k. */ std::vector > monomials; diff --git a/deal.II/base/source/polynomials_bdm.cc b/deal.II/base/source/polynomials_bdm.cc index ef7443035b..782348ee7d 100644 --- a/deal.II/base/source/polynomials_bdm.cc +++ b/deal.II/base/source/polynomials_bdm.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2004, 2005, 2006, 2007, 2008 by the deal.II authors +// Copyright (C) 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -23,14 +23,24 @@ template PolynomialsBDM::PolynomialsBDM (const unsigned int k) : polynomial_space (Polynomials::Legendre::generate_complete_basis(k)), - monomials(1), + monomials((dim==2) ? (1) : (k+2)), n_pols(compute_n_pols(k)), p_values(polynomial_space.n()), p_grads(polynomial_space.n()), p_grad_grads(polynomial_space.n()) { - Assert (dim == 2, ExcNotImplemented()); - monomials[0] = Polynomials::Monomial (k+1); + switch(dim) + { + case 2: + monomials[0] = Polynomials::Monomial (k+1); + break; + case 3: + for (unsigned int i=0;i (i); + break; + default: + Assert(false, ExcNotImplemented()); + } } @@ -60,7 +70,7 @@ PolynomialsBDM::compute (const Point &unit_point, // 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::compute (const Point &unit_point, } } - - // Let's hope this is not the transpose + std::fill(grads.begin(), grads.end(), Tensor<2,dim>()); for (unsigned int i=0;i()); for (unsigned int i=0;i::compute (const Point &unit_point, // This is the first polynomial not // covered by the P_k subspace - const unsigned int start = dim*n_sub; + unsigned int start = dim*n_sub; // Store values of auxiliary - // polynomials and their + // polynomials and their three // derivatives - std::vector monoval0(3); - std::vector monoval1(3); - - monomials[0].value(unit_point(0), monoval0); - monomials[0].value(unit_point(1), monoval1); - if (values.size() != 0) - { - values[start][0] = monoval0[0]; - values[start][1] = -unit_point(1) * monoval0[1]; - values[start+1][0] = -unit_point(0) * monoval1[1]; - values[start+1][1] = monoval1[0]; - } - if (grads.size() != 0) + std::vector > monovali(dim, std::vector(4)); + std::vector > monovalk(dim, std::vector(4)); + + if (dim == 2) { - grads[start][0][0] = monoval0[1]; - grads[start][0][1] = 0.; - grads[start][1][0] = -unit_point(1) * monoval0[2]; - grads[start][1][1] = -monoval0[1]; - grads[start+1][0][0] = -monoval1[1]; - grads[start+1][0][1] = -unit_point(0) * monoval1[2]; - grads[start+1][1][0] = 0.; - grads[start+1][1][1] = monoval1[1]; + for (unsigned int d=0;d