]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
missing functions in PolynomialSpace coded
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 8 Jul 2002 13:39:43 +0000 (13:39 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 8 Jul 2002 13:39:43 +0000 (13:39 +0000)
git-svn-id: https://svn.dealii.org/trunk@6223 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 2cd7a37a1f56a1cd5ab8e39ca4dfe276ca66ddb2..586642a612c59d609d7d3f22aa86b16e43b11184 100644 (file)
@@ -152,6 +152,20 @@ class PolynomialSpace
                                      */
     const unsigned int n_pols;
 
+                                    /**
+                                     * Compute numbers in x, y and z
+                                     * direction. Given an index
+                                     * @p{n} in the d-dimensional
+                                     * polynomial space, compute the
+                                     * indices i,j,k such that
+                                     * @p{p_n(x,y,z) =
+                                     * p_i(x)p_j(y)p_k(z)}.
+                                     */
+    void compute_index(unsigned int n,
+                      unsigned int& nx,
+                      unsigned int& ny,
+                      unsigned int& nz) const;
+    
                                     /**
                                      * Static function used in the
                                      * constructor to compute the
index 460a9f9d9da19200fbbe97aaa0eb5dbf724b3fcf..1c01c060b1d073adf4df70babfc716c1e3eec6a7 100644 (file)
@@ -30,34 +30,103 @@ PolynomialSpace<dim>::compute_n_pols (const unsigned int n)
 };
 
 
+template <int dim>
+void
+PolynomialSpace<dim>::compute_index(unsigned int n,
+                                   unsigned int& nx,
+                                   unsigned int& ny,
+                                   unsigned int& nz) const
+{
+  const unsigned int n_1d=polynomials.size();
+  unsigned int k=0;
+  for (unsigned int iz=0;iz<((dim>2) ? n_1d : 1);++iz)
+    for (unsigned int iy=0;iy<((dim>1) ? n_1d-iz : 1);++iy)
+      for (unsigned int ix=0; ix<n_1d-iy-iz; ++ix)
+       if (k++ == n)
+         {
+           nz = iz;
+           ny = iy;
+           nx = ix;
+           return;
+         }
+}
+
 
 template <int dim>
 double
-PolynomialSpace<dim>::compute_value(const unsigned int /*i*/,
-                                   const Point<dim> & /*p*/) const
+PolynomialSpace<dim>::compute_value(const unsigned int i,
+                                   const Point<dim> & p) const
 {
-  Assert(false, ExcNotImplemented());
-  return 0.;
+  unsigned int ix = 0;
+  unsigned int iy = 0;
+  unsigned int iz = 0;
+  compute_index(i,ix,iy,iz);
+  
+  double result = polynomials[ix].value(p(0));
+  if (dim>1)
+    result *= polynomials[iy].value(p(1));
+  if (dim>2)
+    result *= polynomials[iz].value(p(2));
+  return result;
 }
 
   
 template <int dim>
 Tensor<1,dim>
-PolynomialSpace<dim>::compute_grad(const unsigned int /*i*/,
-                                  const Point<dim> &/*p*/) const
+PolynomialSpace<dim>::compute_grad(const unsigned int i,
+                                  const Point<dim> &p) const
 {
-  Assert(false, ExcNotImplemented());
-  return Tensor<1,dim>();
+  unsigned int ix[3];
+  compute_index(i,ix[0],ix[1],ix[2]);
+  
+  Tensor<1,dim> result;
+  for (unsigned int d=0;d<dim;++d)
+    result[d] = 1.;
+  
+  std::vector<double> v(2);
+  for (unsigned int d=0;d<dim;++d)
+    {
+      polynomials[ix[d]].value(p(d), v);
+      result[d] *= v[1];
+      for (unsigned int d1=0;d1<dim;++d1)
+       if (d1 != d)
+         result[d1] *= v[0];
+    }
+  return result;
 }
 
 
 template <int dim>
 Tensor<2,dim>
-PolynomialSpace<dim>::compute_grad_grad(const unsigned int /*i*/,
-                                       const Point<dim> &/*p*/) const
+PolynomialSpace<dim>::compute_grad_grad(const unsigned int i,
+                                       const Point<dim> &p) const
 {
-  Assert(false, ExcNotImplemented());
-  return Tensor<2,dim>();
+  unsigned int ix[3];
+  compute_index(i,ix[0],ix[1],ix[2]);
+  
+  Tensor<2,dim> result;
+  for (unsigned int d=0;d<dim;++d)
+    for (unsigned int d1=0;d1<dim;++d1)
+      result[d][d1] = 1.;
+  
+  std::vector<double> v(3);
+  for (unsigned int d=0;d<dim;++d)
+    {
+      polynomials[ix[d]].value(p(d), v);
+      result[d][d] *= v[2];
+      for (unsigned int d1=0;d1<dim;++d1)
+       {
+         if (d1 != d)
+           {
+             result[d][d1] *= v[1];
+             result[d1][d] *= v[1];
+             for (unsigned int d2=0;d2<dim;++d2)
+               if (d2 != d)
+                 result[d1][d2] *= v[0];
+           }
+       }
+    }
+  return result;
 }
 
 

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.