]> https://gitweb.dealii.org/ - dealii.git/commitdiff
PolynomialsBDM in 3D
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Sun, 10 May 2009 21:34:46 +0000 (21:34 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Sun, 10 May 2009 21:34:46 +0000 (21:34 +0000)
git-svn-id: https://svn.dealii.org/trunk@18829 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/polynomials_bdm.h
deal.II/base/source/polynomials_bdm.cc

index e808736c2560dcb9a748143893c735076b9b049e..89ed5a71243887cbe39c0bfce6fff9b0887308e9 100644 (file)
@@ -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:
  *
  * <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
@@ -137,7 +137,11 @@ 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;
     
index ef7443035bb259d3dc1ac48d2b2ece8408823007..782348ee7d904892401b653a951604c8cb2aaba4 100644 (file)
@@ -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 <int dim>
 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());
+    }
 }
 
 
@@ -60,7 +70,7 @@ PolynomialsBDM<dim>::compute (const Point<dim>            &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<p_values.size();++i)
     {
@@ -70,16 +80,14 @@ PolynomialsBDM<dim>::compute (const Point<dim>            &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<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)
     {
@@ -89,37 +97,136 @@ PolynomialsBDM<dim>::compute (const Point<dim>            &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<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());
     }
 }
 

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.