// $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
* vector valued polynomials:
*
* <dl>
- * <dt>In 2D:
+ * <dt> In 2D:
* <dd> The 2D-curl of the functions <i>x<sup>k+1</sup>y</i>
* and <i>xy<sup>k+1</sup></i>.
- * <dt>In 3D:
+ * <dt>In 3D:
* <dd> For any <i>i=0,...,k</i> the curls of
* <i>(0,0,xy<sup>i+1</sup>z<sup>k-i</sup>)</i>,
- * <i>(0,x<sup>k-i</sup>yz<sup>i+1</sup>,0)</i> and
- * <i>(x<sup>i+1</sup>y<sup>k-i</sup>z,0,0)</i>
+ * <i>(x<sup>k-i</sup>yz<sup>i+1</sup>,0,0)</i> and
+ * <i>(0,x<sup>i+1</sup>y<sup>k-i</sup>z,0)</i>
* </dl>
*
- * 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 <int dim>
class PolynomialsBDM
const PolynomialSpace<dim> polynomial_space;
/**
- * Storage for monomials
+ * Storage for monomials. In 2D,
+ * this is just the polynomial of
+ * order <i>k</i>. In 3D, we
+ * need all polynomials from
+ * degree zero to <i>k</i>.
*/
std::vector<Polynomials::Polynomial<double> > monomials;
// $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
PolynomialsBDM<dim>::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<double> (k+1);
+ switch(dim)
+ {
+ case 2:
+ monomials[0] = Polynomials::Monomial<double> (k+1);
+ break;
+ case 3:
+ for (unsigned int i=0;i<monomials.size();++i)
+ monomials[i] = Polynomials::Monomial<double> (i);
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
}
// 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)
{
}
}
-
- // 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)
{
// 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<double> monoval0(3);
- std::vector<double> 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<std::vector<double> > monovali(dim, std::vector<double>(4));
+ std::vector<std::vector<double> > monovalk(dim, std::vector<double>(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<dim;++d)
+ monomials[0].value(unit_point(d), monovali[d]);
+ if (values.size() != 0)
+ {
+ values[start][0] = monovali[0][0];
+ values[start][1] = -unit_point(1) * monovali[0][1];
+ values[start+1][0] = -unit_point(0) * monovali[1][1];
+ values[start+1][1] = monovali[1][0];
+ }
+ if (grads.size() != 0)
+ {
+ grads[start][0][0] = monovali[0][1];
+ grads[start][0][1] = 0.;
+ grads[start][1][0] = -unit_point(1) * monovali[0][2];
+ grads[start][1][1] = -monovali[0][1];
+ grads[start+1][0][0] = -monovali[1][1];
+ grads[start+1][0][1] = -unit_point(0) * monovali[1][2];
+ grads[start+1][1][0] = 0.;
+ grads[start+1][1][1] = monovali[1][1];
+ }
+ if (grad_grads.size() != 0)
+ {
+ grad_grads[start][0][0][0] = monovali[0][2];
+ grad_grads[start][0][0][1] = 0.;
+ grad_grads[start][0][1][0] = 0.;
+ grad_grads[start][0][1][1] = 0.;
+ grad_grads[start][1][0][0] = -unit_point(1) * monovali[0][3];
+ grad_grads[start][1][0][1] = -monovali[0][2];
+ grad_grads[start][1][1][0] = -monovali[0][2];
+ grad_grads[start][1][1][1] = 0.;
+ grad_grads[start+1][0][0][0] = 0;
+ grad_grads[start+1][0][0][1] = -monovali[1][2];
+ grad_grads[start+1][0][1][0] = -monovali[1][2];
+ grad_grads[start+1][0][1][1] = -unit_point(0) * monovali[1][3];
+ grad_grads[start+1][1][0][0] = 0.;
+ grad_grads[start+1][1][0][1] = 0.;
+ grad_grads[start+1][1][1][0] = 0.;
+ grad_grads[start+1][1][1][1] = monovali[1][2];
+ }
}
- if (grad_grads.size() != 0)
+ else // dim == 3
{
- Assert(false,ExcNotImplemented());
+ // The number of curls in each
+ // component. Note that the
+ // table in BrezziFortin91 has
+ // a typo, but the text has the
+ // right basis
+
+ // Note that the next basis
+ // function is always obtained
+ // from the previous by cyclic
+ // rotation of the coordinates
+ const unsigned int n_curls = monomials.size() - 1;
+ for (unsigned int i=0;i<n_curls;++i, start+=dim)
+ {
+ for (unsigned int d=0;d<dim;++d)
+ {
+ // p(t) = t^(i+1)
+ monomials[i+1].value(unit_point(d), monovali[d]);
+ // q(t) = t^(k-i)
+ monomials[degree()-i].value(unit_point(d), monovalk[d]);
+ }
+ if (values.size() != 0)
+ {
+ // x p'(y) q(z)
+ values[start][0] = unit_point(0) * monovali[1][1] * monovalk[2][0];
+ // - p(y) q(z)
+ values[start][1] = -monovali[1][0] * monovalk[2][0];
+ values[start][2] = 0.;
+
+ // y p'(z) q(x)
+ values[start+1][1] = unit_point(1) * monovali[2][1] * monovalk[0][0];
+ // - p(z) q(x)
+ values[start+1][2] = -monovali[2][0] * monovalk[0][0];
+ values[start+1][0] = 0.;
+
+ // z p'(x) q(y)
+ values[start+2][2] = unit_point(2) * monovali[0][1] * monovalk[1][0];
+ // -p(x) q(y)
+ values[start+2][0] = -monovali[0][0] * monovalk[1][0];
+ values[start+2][1] = 0.;
+ }
+ if (grads.size() != 0)
+ {
+ grads[start][0][0] = monovali[1][1] * monovalk[2][0];
+ grads[start][0][1] = unit_point(0) * monovali[1][2] * monovalk[2][0];
+ grads[start][0][2] = unit_point(0) * monovali[1][1] * monovalk[2][1];
+ grads[start][1][0] = 0.;
+ grads[start][1][1] = -monovali[1][1] * monovalk[2][0];
+ grads[start][1][2] = -monovali[1][0] * monovalk[2][1];
+ grads[start+2][2][0] = 0.;
+ grads[start+2][2][1] = 0.;
+ grads[start+2][2][2] = 0.;
+
+ grads[start+1][1][1] = monovali[2][1] * monovalk[0][0];
+ grads[start+1][1][2] = unit_point(1) * monovali[2][2] * monovalk[0][0];
+ grads[start+1][1][0] = unit_point(1) * monovali[2][1] * monovalk[0][1];
+ grads[start+1][2][1] = 0.;
+ grads[start+1][2][2] = -monovali[2][1] * monovalk[0][0];
+ grads[start+1][2][0] = -monovali[2][0] * monovalk[0][1];
+ grads[start+1][0][1] = 0.;
+ grads[start+1][0][2] = 0.;
+ grads[start+1][0][0] = 0.;
+
+ grads[start+2][2][2] = monovali[0][1] * monovalk[1][0];
+ grads[start+2][2][0] = unit_point(2) * monovali[0][2] * monovalk[1][0];
+ grads[start+2][2][1] = unit_point(2) * monovali[0][1] * monovalk[1][1];
+ grads[start+2][0][2] = 0.;
+ grads[start+2][0][0] = -monovali[0][1] * monovalk[1][0];
+ grads[start+2][0][1] = -monovali[0][0] * monovalk[1][1];
+ grads[start+2][1][2] = 0.;
+ grads[start+2][1][0] = 0.;
+ grads[start+2][1][1] = 0.;
+ }
+ if (grad_grads.size() != 0)
+ {
+ Assert(false,ExcNotImplemented());
+ }
+ }
+ Assert(start == n_pols, ExcInternalError());
}
}