]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
move all boost::any *_nth_derivative() to Tensor<n,dim>
authorjanssen <janssen@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 9 Aug 2012 11:47:14 +0000 (11:47 +0000)
committerjanssen <janssen@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 9 Aug 2012 11:47:14 +0000 (11:47 +0000)
*_nth_derivative() and the boost::any functions become internal

git-svn-id: https://svn.dealii.org/branches/branch_higher_derivatives@25803 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/base/tensor_product_polynomials.h
deal.II/include/deal.II/fe/fe.h
deal.II/include/deal.II/fe/fe_poly.h
deal.II/include/deal.II/fe/fe_poly.templates.h
deal.II/include/deal.II/fe/fe_poly_tensor.h
deal.II/include/deal.II/fe/fe_values.h
deal.II/source/base/tensor_product_polynomials.cc
deal.II/source/fe/fe.cc
deal.II/source/fe/fe_poly.cc
deal.II/source/fe/fe_poly_tensor.cc
deal.II/source/fe/fe_values.cc

index 73c5f21eedf596f981b1de30dd249d3afa9994b6..e0c02f596d8013ae3d56aada257094925142690a 100644 (file)
@@ -129,6 +129,7 @@ class TensorProductPolynomials
                                       * compute_value(),
                                       * compute_grad() or
                                       * compute_hessian()
+                                      * compute_nth_derivative()
                                       * functions, see below, in a
                                       * loop over all tensor product
                                       * polynomials.
@@ -138,6 +139,12 @@ class TensorProductPolynomials
                   std::vector<Tensor<1,dim> > &grads,
                   std::vector<Tensor<2,dim> > &hessians) const;
 
+    void compute (const Point<dim>                      &unit_point,
+                  std::vector<double>                   &values,
+                  std::vector<Tensor<1,dim> >           &grads,
+                  std::vector<Tensor<2,dim> >           &hessians,
+                  std::vector<std::vector<boost::any> > &nth_derivatives) const;
+
                                      /**
                                       * Computes the value of the
                                       * <tt>i</tt>th tensor product
@@ -373,6 +380,12 @@ class AnisotropicPolynomials
                   std::vector<Tensor<1,dim> > &grads,
                   std::vector<Tensor<2,dim> > &hessians) const;
 
+    void compute (const Point<dim>            &unit_point,
+                  std::vector<double>         &values,
+                  std::vector<Tensor<1,dim> > &grads,
+                  std::vector<Tensor<2,dim> > &hessians,
+                  std::vector<std::vector<boost::any> > &nth_derivatives) const;
+
                                      /**
                                       * Computes the value of the
                                       * <tt>i</tt>th tensor product
index 3813fde941e61feeb6e95bcd2756e2d3d35b0506..1f43870c3aa5e549a13a041c1652a3726ecaf21d 100644 (file)
@@ -637,11 +637,6 @@ class FiniteElement : public Subscriptor,
      * @param nth_derivative. The return value
      * @return can only store Tensor<nth_derivative,dim>.
      */
-    virtual
-    boost::any
-    shape_nth_derivative_internal (const unsigned int i,
-                                  const Point<dim> &p,
-                                  const unsigned int nth_derivative) const;
 
     template <int n>
     Tensor<n,dim>
@@ -663,12 +658,17 @@ class FiniteElement : public Subscriptor,
      * @param nth_derivative. The return value
      * @return can only store Tensor<nth_derivative,dim>.
      */
-    virtual
-    boost::any
+
+    template <int n>
+    Tensor<n,dim>
     shape_nth_derivative_component (const unsigned int i,
                                     const Point<dim> &p,
                                     const unsigned int component,
-                                    const unsigned int nth_derivative) const;
+                                    const unsigned int nth_derivative) const
+      {
+       return boost::any_cast<Tensor<n,dim> >
+         (shape_nth_derivative_component_internal (i, p, nth_derivative));
+      }
 
                                      /**
                                       * Check for non-zero values on a
@@ -2664,6 +2664,19 @@ class FiniteElement : public Subscriptor,
                             typename Mapping<dim,spacedim>::InternalDataBase &fe_internal,
                             FEValuesData<dim,spacedim>                    &data) const = 0;
 
+    virtual
+    boost::any
+    shape_nth_derivative_internal (const unsigned int i,
+                                  const Point<dim> &p,
+                                  const unsigned int nth_derivative) const;
+
+    virtual
+    boost::any
+    shape_nth_derivative_component_internal (const unsigned int i,
+                                    const Point<dim> &p,
+                                    const unsigned int component,
+                                    const unsigned int nth_derivative) const;
+
     friend class InternalDataBase;
     friend class FEValuesBase<dim,spacedim>;
     friend class FEValues<dim,spacedim>;
index 2f95eaa759fccf531012eaeb835a3f39ef2fb3a3..238b014e5c6b48ce043ad9a24764a823b1d4b8c9 100644 (file)
@@ -211,11 +211,16 @@ class FE_Poly : public FiniteElement<dim,spacedim>
                                       * for more information about the
                                       * semantics of this function.
                                       */
-    virtual
-    boost::any
+
+    template <int n>
+    Tensor<n,dim>
     shape_nth_derivative (const unsigned int i,
                           const Point<dim> &p,
-                          const unsigned int nth_derivative) const;
+                          const unsigned int nth_derivative) const
+      {
+       return boost::any_cast<Tensor<n,dim> >
+         (shape_nth_derivative_internal (i, p, nth_derivative));
+      }
 
                                      /**
                                       * Return the tensor of the 
@@ -236,12 +241,17 @@ class FE_Poly : public FiniteElement<dim,spacedim>
                                       * were called, provided that the
                                       * specified component is zero.
                                       */
-    virtual
-    boost::any
+
+    template <int n>
+    Tensor<n,dim>
     shape_nth_derivative_component (const unsigned int i,
                                     const Point<dim> &p,
                                     const unsigned int component,
-                                    const unsigned int nth_derivative) const;
+                                    const unsigned int nth_derivative) const
+      {
+       return boost::any_cast<Tensor<n,dim> >
+         (shape_nth_derivative_component_internal (i, p, nth_derivative));
+      }
   protected:
 
     virtual
@@ -278,6 +288,18 @@ class FE_Poly : public FiniteElement<dim,spacedim>
                             typename Mapping<dim,spacedim>::InternalDataBase      &fe_internal,
                             FEValuesData<dim,spacedim>& data) const ;
 
+    virtual
+    boost::any
+    shape_nth_derivative_internal (const unsigned int i,
+                          const Point<dim> &p,
+                          const unsigned int nth_derivative) const;
+
+    virtual
+    boost::any
+    shape_nth_derivative_component_internal (const unsigned int i,
+                          const Point<dim> &p,
+                          const unsigned int nth_derivative,
+                          const unsigned int component) const;
 
                                      /**
                                       * Determine the values that need
@@ -424,6 +446,8 @@ class FE_Poly : public FiniteElement<dim,spacedim>
         std::vector<std::vector<Tensor<7,dim> > > shape_7th_derivatives;
         std::vector<std::vector<Tensor<8,dim> > > shape_8th_derivatives;
         std::vector<std::vector<Tensor<9,dim> > > shape_9th_derivatives;
+
+        std::vector<std::vector<std::vector<boost::any> > > shape_nth_derivatives;
     };
 
                                      /**
index bac8eaf2311b1a3113aa3e39215b62092ccf1ed4..cefcbcaa534d72b3da2870446c91a9484253de35 100644 (file)
@@ -114,7 +114,7 @@ FE_Poly<POLY,dim,spacedim>::shape_hessian_component (const unsigned int i,
 
 template <class POLY, int dim, int spacedim>
 boost::any
-FE_Poly<POLY,dim,spacedim>::shape_nth_derivative (const unsigned int i,
+FE_Poly<POLY,dim,spacedim>::shape_nth_derivative_internal (const unsigned int i,
                                                   const Point<dim> &p,
                                                   const unsigned int nth_derivative) const
 {
@@ -126,7 +126,7 @@ FE_Poly<POLY,dim,spacedim>::shape_nth_derivative (const unsigned int i,
 
 template <class POLY, int dim, int spacedim>
 boost::any
-FE_Poly<POLY,dim,spacedim>::shape_nth_derivative_component (const unsigned int i,
+FE_Poly<POLY,dim,spacedim>::shape_nth_derivative_component_internal (const unsigned int i,
                                                             const Point<dim> &p,
                                                             const unsigned int component,
                                                             const unsigned int nth_derivative) const
index 14a24361a8523c3708bd9f8e08cc4be709e4fa27..2aed3b68097d0639befc2464be28207103dd754a 100644 (file)
@@ -169,14 +169,22 @@ class FE_PolyTensor : public FiniteElement<dim,spacedim>
                                       * vector valued, an exception is
                                       * thrown.
                                       */
-    virtual boost::any shape_nth_derivative (const unsigned int  i,
+    template <int n> Tensor<n,dim> shape_nth_derivative (const unsigned int  i,
                                              const Point<dim> &p,
-                                             const unsigned int nth_derivative) const;
+                                             const unsigned int nth_derivative) const
+      {
+       return boost::any_cast<Tensor<n,dim> >
+         (shape_nth_derivative_internal (i, p, nth_derivative));
+      }
 
-    virtual boost::any shape_nth_derivative_component (const unsigned int i,
+    template <int n> Tensor<n,dim> shape_nth_derivative_component (const unsigned int i,
                                                 const Point<dim> &p,
                                                 const unsigned int component,
-                                                const unsigned int nth_derivative) const;
+                                                const unsigned int nth_derivative) const
+      {
+       return boost::any_cast<Tensor<n,dim> >
+         (shape_nth_derivative_component_internal (i, p, nth_derivative));
+      }
 
                                      /**
                                       * Given <tt>flags</tt>,
@@ -244,6 +252,17 @@ class FE_PolyTensor : public FiniteElement<dim,spacedim>
                             typename Mapping<dim,spacedim>::InternalDataBase      &fe_internal,
                             FEValuesData<dim,spacedim>& data) const ;
 
+    virtual boost::any shape_nth_derivative_internal (const unsigned int  i,
+                                             const Point<dim> &p,
+                                             const unsigned int nth_derivative) const;
+
+
+    virtual boost::any shape_nth_derivative_component_internal (const unsigned int i,
+                                                const Point<dim> &p,
+                                                const unsigned int component,
+                                                const unsigned int nth_derivative) const;
+
+
                                      /**
                                       * Fields of cell-independent
                                       * data for FE_PolyTensor. Stores
@@ -283,6 +302,16 @@ class FE_PolyTensor : public FiniteElement<dim,spacedim>
                                           * point.
                                           */
         std::vector< std::vector< DerivativeForm<1, dim, spacedim> > > shape_grads;
+        std::vector<std::vector<DerivativeForm<2,dim, spacedim> > > shape_hessians;
+        std::vector<std::vector<DerivativeForm<3,dim, spacedim> > > shape_3rd_derivatives;
+        std::vector<std::vector<DerivativeForm<4,dim, spacedim> > > shape_4th_derivatives;
+        std::vector<std::vector<DerivativeForm<5,dim, spacedim> > > shape_5th_derivatives;
+        std::vector<std::vector<DerivativeForm<6,dim, spacedim> > > shape_6th_derivatives;
+        std::vector<std::vector<DerivativeForm<7,dim, spacedim> > > shape_7th_derivatives;
+        std::vector<std::vector<DerivativeForm<8,dim, spacedim> > > shape_8th_derivatives;
+        std::vector<std::vector<DerivativeForm<9,dim, spacedim> > > shape_9th_derivatives;
+
+        std::vector<std::vector<std::vector<boost::any> > > shape_nth_derivatives;
     };
 
                                      /**
index 8ab63b280ebd0d2bf3038761d0ef627bf8d39bea..69e9bde9d67fe0c8f2433b9e3faa6e8129a247a0 100644 (file)
@@ -1474,6 +1474,15 @@ class FEValuesData
                                       * derivatives.
                                       */
     typedef std::vector<std::vector<Tensor<2,spacedim> > > HessianVector;
+    typedef std::vector<std::vector<Tensor<3,spacedim> > > ThirdDerivativeVector;
+    typedef std::vector<std::vector<Tensor<4,spacedim> > > FourthDerivativeVector;
+    typedef std::vector<std::vector<Tensor<5,spacedim> > > FifthDerivativeVector;
+    typedef std::vector<std::vector<Tensor<6,spacedim> > > SixthDerivativeVector;
+    typedef std::vector<std::vector<Tensor<7,spacedim> > > SeventhDerivativeVector;
+    typedef std::vector<std::vector<Tensor<8,spacedim> > > EighthDerivativeVector;
+    typedef std::vector<std::vector<Tensor<9,spacedim> > > NinthDerivativeVector;
+
+    typedef std::vector<std::vector<std::vector<boost::any> > > NthDerivativeVector;
 
                                      /**
                                       * Store the values of the shape
@@ -1503,6 +1512,14 @@ class FEValuesData
                                       * this field.
                                       */
     HessianVector shape_hessians;
+    ThirdDerivativeVector shape_3rd_derivatives;
+    FourthDerivativeVector shape_4th_derivatives;
+    FifthDerivativeVector shape_5th_derivatives;
+    SixthDerivativeVector shape_6th_derivatives;
+    SeventhDerivativeVector shape_7th_derivatives;
+    EighthDerivativeVector shape_8th_derivatives;
+    NinthDerivativeVector shape_9th_derivatives;
+    NthDerivativeVector shape_nth_derivatives;
 
                                      /**
                                       * Store an array of weights
@@ -2019,6 +2036,81 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
                                     const unsigned int point_no,
                                     const unsigned int component) const;
 
+                                     /**
+                                      * <tt>nth_derivative</tt> derivatives of
+                                      * the <tt>function_no</tt>th shape function at
+                                      * the <tt>point_no</tt>th quadrature point
+                                      * with respect to real cell
+                                      * coordinates. If you want to
+                                      * get the derivatives in one of
+                                      * the coordinate directions, use
+                                      * the appropriate function of
+                                      * the Tensor class to
+                                      * extract one component. Since
+                                      * only a reference to the
+                                      * derivative values is returned,
+                                      * there should be no major
+                                      * performance drawback.
+                                      *
+                                      * If the shape function is
+                                      * vector-valued, then this
+                                      * returns the only non-zero
+                                      * component. If the shape
+                                      * function has more than one
+                                      * non-zero component (i.e. it is
+                                      * not primitive), then throw an
+                                      * exception of type
+                                      * ExcShapeFunctionNotPrimitive. In
+                                      * that case, use the
+                                      * shape_nth_derivative_component()
+                                      * function.
+                                      *
+                                      * The same holds for the arguments
+                                      * of this function as for the
+                                      * shape_value() function.
+                                      */
+    const boost::any &
+    shape_nth_derivative (const unsigned int function_no,
+                   const unsigned int point_no,
+                   const unsigned int nth_derivative) const;
+
+
+                                     /**
+                                      * Return one vector component of
+                                      * the gradient of a shape
+                                      * function at a quadrature
+                                      * point. If the finite element
+                                      * is scalar, then only component
+                                      * zero is allowed and the return
+                                      * value equals that of the
+                                      * shape_nth_derivative()
+                                      * function. If the finite
+                                      * element is vector valued but
+                                      * all shape functions are
+                                      * primitive (i.e. they are
+                                      * non-zero in only one
+                                      * component), then the value
+                                      * returned by
+                                      * shape_nth_derivative()
+                                      * equals that of this function
+                                      * for exactly one
+                                      * component. This function is
+                                      * therefore only of greater
+                                      * interest if the shape function
+                                      * is not primitive, but then it
+                                      * is necessary since the other
+                                      * function cannot be used.
+                                      *
+                                      * The same holds for the arguments
+                                      * of this function as for the
+                                      * shape_value_component() function.
+                                      */
+    boost::any
+    shape_nth_derivative_component (const unsigned int function_no,
+                             const unsigned int point_no,
+                             const unsigned int component,
+                             const unsigned int nth_derivative) const;
+
 
                                      //@}
                                      /// @name Access to values of global finite element fields
@@ -4931,6 +5023,112 @@ FEValuesBase<dim,spacedim>::shape_2nd_derivative_component (const unsigned int i
 
 
 
+template <int dim, int spacedim>
+inline
+const boost::any &
+FEValuesBase<dim,spacedim>::shape_nth_derivative (const unsigned int i,
+                                           const unsigned int j,
+                                           const unsigned int nth_derivative) const
+{
+  Assert (i < fe->dofs_per_cell,
+          ExcIndexRange (i, 0, fe->dofs_per_cell));
+  Assert (this->update_flags & update_nth_derivatives(nth_derivative),
+          ExcAccessToUninitializedField());
+  Assert (fe->is_primitive (i),
+          ExcShapeFunctionNotPrimitive(i));
+  Assert (nth_derivative<this->shape_nth_derivatives.size(),
+          ExcIndexRange (nth_derivative, 0, this->shape_nth_derivatives.size()));
+  Assert (i<this->shape_nth_derivatives[nth_derivative].size(),
+          ExcIndexRange (i, 0, this->shape_nth_derivatives[nth_derivative].size()));
+  Assert (j<this->shape_nth_derivatives[nth_derivative][0].size(),
+          ExcIndexRange (j, 0, this->shape_nth_derivatives[nth_derivative][0].size()));
+
+                                   // if the entire FE is primitive,
+                                   // then we can take a short-cut:
+  if (fe->is_primitive())
+    return this->shape_nth_derivatives[nth_derivative][i][j];
+  else
+                                     // otherwise, use the mapping
+                                     // between shape function numbers
+                                     // and rows. note that by the
+                                     // assertions above, we know that
+                                     // this particular shape function
+                                     // is primitive, so there is no
+                                     // question to which vector
+                                     // component the call of this
+                                     // function refers
+    return this->shape_nth_derivatives[nth_derivative][this->shape_function_to_row_table[i]][j];
+}
+
+
+
+
+template <int dim, int spacedim>
+inline
+boost::any
+FEValuesBase<dim,spacedim>::shape_nth_derivative_component (const unsigned int i,
+                                                     const unsigned int j,
+                                                     const unsigned int component,
+                                                     const unsigned int nth_derivative) const
+{
+  Assert (i < fe->dofs_per_cell,
+          ExcIndexRange (i, 0, fe->dofs_per_cell));
+  Assert (this->update_flags & update_nth_derivatives(nth_derivative),
+          ExcAccessToUninitializedField());
+  Assert (component < fe->n_components(),
+          ExcIndexRange(component, 0, fe->n_components()));
+
+                                   // if this particular shape
+                                   // function is primitive, then we
+                                   // can take a short-cut by checking
+                                   // whether the requested component
+                                   // is the only non-zero one (note
+                                   // that calling
+                                   // system_to_component_table only
+                                   // works if the shape function is
+                                   // primitive):
+  if (fe->is_primitive(i))
+    {
+      if (component == fe->system_to_component_index(i).first)
+        return this->shape_nth_derivatives[nth_derivative][this->shape_function_to_row_table[i]][j];
+      else
+        return boost::any ();
+    }
+  else
+    {
+                                       // no, this shape function is
+                                       // not primitive. then we have
+                                       // to loop over its components
+                                       // to find the corresponding
+                                       // row in the arrays of this
+                                       // object. before that check
+                                       // whether the shape function
+                                       // is non-zero at all within
+                                       // this component:
+      if (fe->get_nonzero_components(i)[component] == false)
+        return boost::any ();
+
+                                       // count how many non-zero
+                                       // component the shape function
+                                       // has before the one we are
+                                       // looking for, and add this to
+                                       // the offset of the first
+                                       // non-zero component of this
+                                       // shape function in the arrays
+                                       // we index presently:
+      const unsigned int
+        row = (this->shape_function_to_row_table[i]
+               +
+               std::count (fe->get_nonzero_components(i).begin(),
+                           fe->get_nonzero_components(i).begin()+component,
+                           true));
+      return this->shape_nth_derivatives[nth_derivative][row][j];
+    }
+}
+
+
+
+
 template <int dim, int spacedim>
 inline
 const FiniteElement<dim,spacedim> &
index 7b052f0f6730a9db278ba60bd101beb1c815894e..e0144ccbbba85e55f1f559781876379e4d4110f0 100644 (file)
@@ -354,6 +354,149 @@ compute (const Point<dim>            &p,
 
 
 
+
+template <int dim>
+void
+TensorProductPolynomials<dim>::
+compute (const Point<dim>            &p,
+         std::vector<double>         &values,
+         std::vector<Tensor<1,dim> > &grads,
+         std::vector<Tensor<2,dim> > &hessians,
+         std::vector<std::vector<boost::any> > &nth_derivatives) const
+{
+  Assert (values.size()==n_tensor_pols    || values.size()==0,
+          ExcDimensionMismatch2(values.size(), n_tensor_pols, 0));
+  Assert (grads.size()==n_tensor_pols     || grads.size()==0,
+          ExcDimensionMismatch2(grads.size(), n_tensor_pols, 0));
+  Assert (hessians.size()==n_tensor_pols|| hessians.size()==0,
+          ExcDimensionMismatch2(hessians.size(), n_tensor_pols, 0));
+
+  for(unsigned int i=0; i<nth_derivatives.size(); ++i)
+  {
+    //TODO: check all the other derivatives BJ
+    //Assert (hessians.size()==n_tensor_pols|| hessians.size()==0,
+    //    ExcDimensionMismatch2(hessians.size(), n_tensor_pols, 0));
+  }
+
+  const bool update_values     = (values.size() == n_tensor_pols),
+             update_grads      = (grads.size()==n_tensor_pols),
+             update_hessians = (hessians.size()==n_tensor_pols),
+             update_3rd_derivatives = (nth_derivatives[3].size()==n_tensor_pols),
+             update_4th_derivatives = (nth_derivatives[4].size()==n_tensor_pols),
+             update_5th_derivatives = (nth_derivatives[5].size()==n_tensor_pols),
+             update_6th_derivatives = (nth_derivatives[6].size()==n_tensor_pols),
+             update_7th_derivatives = (nth_derivatives[7].size()==n_tensor_pols),
+             update_8th_derivatives = (nth_derivatives[8].size()==n_tensor_pols),
+             update_9th_derivatives = (nth_derivatives[9].size()==n_tensor_pols);
+
+                                   // check how many
+                                   // values/derivatives we have to
+                                   // compute
+  unsigned int n_values_and_derivatives = 0;
+  if (update_values)
+    n_values_and_derivatives = 1;
+  if (update_grads)
+    n_values_and_derivatives = 2;
+  if (update_hessians)
+    n_values_and_derivatives = 3;
+  if (update_3rd_derivatives)
+    n_values_and_derivatives = 4;
+  if (update_4th_derivatives)
+    n_values_and_derivatives = 5;
+  if (update_5th_derivatives)
+    n_values_and_derivatives = 6;
+  if (update_6th_derivatives)
+    n_values_and_derivatives = 7;
+  if (update_7th_derivatives)
+    n_values_and_derivatives = 8;
+  if (update_8th_derivatives)
+    n_values_and_derivatives = 9;
+  if (update_9th_derivatives)
+    n_values_and_derivatives = 10;
+
+
+                                   // compute the values (and derivatives, if
+                                   // necessary) of all polynomials at this
+                                   // evaluation point. to avoid many
+                                   // reallocation, use one std::vector for
+                                   // polynomial evaluation and store the
+                                   // result as Tensor<1,3> (that has enough
+                                   // fields for any evaluation of values and
+                                   // derivatives)
+  Table<2,Tensor<1,3> > v(dim, polynomials.size());
+  {
+    std::vector<double> tmp (n_values_and_derivatives);
+    for (unsigned int d=0; d<dim; ++d)
+      for (unsigned int i=0; i<polynomials.size(); ++i)
+        {
+          polynomials[i].value(p(d), tmp);
+          for (unsigned int e=0; e<n_values_and_derivatives; ++e)
+            v(d,i)[e] = tmp[e];
+        };
+  }
+
+  for (unsigned int i=0; i<n_tensor_pols; ++i)
+    {
+                                       // first get the
+                                       // one-dimensional indices of
+                                       // this particular tensor
+                                       // product polynomial
+      unsigned int indices[dim];
+      compute_index (i, indices);
+
+      if (update_values)
+        {
+          values[i] = 1;
+          for (unsigned int x=0; x<dim; ++x)
+            values[i] *= v(x,indices[x])[0];
+        }
+
+      if (update_grads)
+        for (unsigned int d=0; d<dim; ++d)
+          {
+            grads[i][d] = 1.;
+            for (unsigned int x=0; x<dim; ++x)
+              grads[i][d] *= v(x,indices[x])[d==x];
+          }
+
+      if (update_hessians)
+        for (unsigned int d1=0; d1<dim; ++d1)
+          for (unsigned int d2=0; d2<dim; ++d2)
+            {
+              hessians[i][d1][d2] = 1.;
+              for (unsigned int x=0; x<dim; ++x)
+                {
+                  unsigned int derivative=0;
+                  if (d1==x || d2==x)
+                    {
+                      if (d1==d2)
+                        derivative=2;
+                      else
+                        derivative=1;
+                    }
+                  hessians[i][d1][d2]
+                    *= v(x,indices[x])[derivative];
+                }
+            }
+      if (update_3rd_derivatives)
+      {/*do something clever here*/}
+      if (update_4th_derivatives)
+      {/*do something clever here*/}
+      if (update_5th_derivatives)
+      {/*do something clever here*/}
+      if (update_6th_derivatives)
+      {/*do something clever here*/}
+      if (update_7th_derivatives)
+      {/*do something clever here*/}
+      if (update_8th_derivatives)
+      {/*do something clever here*/}
+      if (update_9th_derivatives)
+      {/*do something clever here*/}
+    }
+}
+
+
+
 template <int dim>
 unsigned int
 TensorProductPolynomials<dim>::n() const
index e8b465211e5d556a039448370c7fa31eb4597966..f2724c0b6dd4905dc7dc54e35a4ea5ee4428ba83 100644 (file)
@@ -60,7 +60,9 @@ InternalDataBase::initialize_2nd (const FiniteElement<dim,spacedim> *element,
                                    // shifted positions
   std::vector<Point<dim> > diff_points (quadrature.size());
 
+  deallog << "before resize" << std::endl;
   differences.resize(2*dim);
+  deallog << "after resize  " << differences.size() << std::endl;
   for (unsigned int d=0; d<dim; ++d)
     {
       Point<dim> shift;
@@ -306,7 +308,7 @@ FiniteElement<dim,spacedim>::shape_hessian_component (const unsigned int,
 
 template <int dim, int spacedim>
 boost::any
-FiniteElement<dim,spacedim>::shape_nth_derivative (const unsigned int,
+FiniteElement<dim,spacedim>::shape_nth_derivative_internal (const unsigned int,
                                                      const Point<dim> &,
                                                      const unsigned int) const
 {
@@ -318,7 +320,7 @@ FiniteElement<dim,spacedim>::shape_nth_derivative (const unsigned int,
 
 template <int dim, int spacedim>
 boost::any
-FiniteElement<dim,spacedim>::shape_nth_derivative_component(const unsigned int,
+FiniteElement<dim,spacedim>::shape_nth_derivative_component_internal(const unsigned int,
                                                             const Point<dim> &,
                                                             const unsigned int,
                                                             const unsigned int) const
index 0601ef0c55d77078993d79304500f3d53c077a30..1f50582f01b3212414beba9a6fc9d37ff9891a5a 100644 (file)
@@ -51,6 +51,30 @@ FE_Poly<TensorProductPolynomials<1>,1,2>::fill_fe_values
       if (flags & update_gradients && cell_similarity != CellSimilarity::translation)
         mapping.transform(fe_data.shape_gradients[k], data.shape_gradients[k],
                           mapping_data, mapping_covariant);
+
+      if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_3rd_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_4th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_5th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_6th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_7th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_8th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_9th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
     }
 
   if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
@@ -88,6 +112,30 @@ FE_Poly<TensorProductPolynomials<2>,2,3>::fill_fe_values
       if (flags & update_gradients && cell_similarity != CellSimilarity::translation)
         mapping.transform(fe_data.shape_gradients[k], data.shape_gradients[k],
                           mapping_data, mapping_covariant);
+
+      if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_3rd_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_4th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_5th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_6th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_7th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_8th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_9th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
     }
 
   if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
@@ -126,6 +174,30 @@ FE_Poly<PolynomialSpace<1>,1,2>::fill_fe_values (
       if (flags & update_gradients && cell_similarity != CellSimilarity::translation)
         mapping.transform(fe_data.shape_gradients[k], data.shape_gradients[k],
                           mapping_data, mapping_covariant);
+
+      if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_3rd_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_4th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_5th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_6th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_7th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_8th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_9th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
     }
 
   if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
@@ -160,6 +232,30 @@ FE_Poly<PolynomialSpace<2>,2,3>::fill_fe_values
       if (flags & update_gradients && cell_similarity != CellSimilarity::translation)
         mapping.transform(fe_data.shape_gradients[k], data.shape_gradients[k],
                           mapping_data, mapping_covariant);
+
+      if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_3rd_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_4th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_5th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_6th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_7th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_8th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
+
+      if (flags & update_9th_derivatives && cell_similarity != CellSimilarity::translation)
+      { /*do something clever*/ }
     }
 
   if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
index 63de560b0d85febeeafeb118e2fa8e86c54610fa..0d428a0b77be1b03ab053294b190c7c3d8a6ffed 100644 (file)
@@ -238,7 +238,7 @@ FE_PolyTensor<POLY,dim,spacedim>::shape_hessian_component (const unsigned int i,
 
 template <class POLY, int dim, int spacedim>
 boost::any
-FE_PolyTensor<POLY,dim,spacedim>::shape_nth_derivative (const unsigned int, const Point<dim> &, const unsigned int) const
+FE_PolyTensor<POLY,dim,spacedim>::shape_nth_derivative_internal (const unsigned int, const Point<dim> &, const unsigned int) const
 {
   typedef    FiniteElement<dim,spacedim> FEE;
   Assert(false, typename FEE::ExcFENotPrimitive());
@@ -249,7 +249,7 @@ FE_PolyTensor<POLY,dim,spacedim>::shape_nth_derivative (const unsigned int, cons
 
 template <class POLY, int dim, int spacedim>
 boost::any
-FE_PolyTensor<POLY,dim,spacedim>::shape_nth_derivative_component (const unsigned int i,
+FE_PolyTensor<POLY,dim,spacedim>::shape_nth_derivative_component_internal (const unsigned int i,
                                                     const Point<dim> &p,
                                                     const unsigned int component,
                                                     const unsigned int nth_derivative) const
@@ -307,6 +307,19 @@ FE_PolyTensor<POLY,dim,spacedim>::get_data (
   std::vector<Tensor<2,dim> > grads(0);
   std::vector<Tensor<3,dim> > hessians(0);
 
+
+  //lieber jeder vector einzeln??
+  std::vector<Tensor<4,dim> > third_derivatives(0);
+  std::vector<Tensor<5,dim> > fourth_derivatives(0);
+  std::vector<Tensor<6,dim> > fifth_derivatives(0);
+  std::vector<Tensor<7,dim> > sixth_derivatives(0);
+  std::vector<Tensor<8,dim> > seventh_derivatives(0);
+  std::vector<Tensor<9,dim> > eighth_derivatives(0);
+  std::vector<Tensor<10,dim> > ninth_derivatives(0);
+
+  //oder einer für alles??
+  std::vector<std::vector<boost::any> > nth_derivatives(10);
+
                                    // initialize fields only if really
                                    // necessary. otherwise, don't
                                    // allocate memory
@@ -332,8 +345,86 @@ FE_PolyTensor<POLY,dim,spacedim>::get_data (
                                    // that
   if (flags & update_hessians)
     {
-//      hessians.resize (this->dofs_per_cell);
+      hessians.resize (this->dofs_per_cell);
+      //old version BJ
       data->initialize_2nd (this, mapping, quadrature);
+
+      //new version BJ
+      //data->shape_hessians.resize (this->dofs_per_cell);
+      //for (unsigned int i=0;i<this->dofs_per_cell;++i)
+      //  data->shape_hessians[i].resize (n_q_points);
+    }
+
+  if (flags & update_3rd_derivatives)
+    {
+      nth_derivatives[3].resize (this->dofs_per_cell);
+      data->shape_3rd_derivatives.resize (this->dofs_per_cell);
+      for (unsigned int i=0;i<this->dofs_per_cell;++i)
+        data->shape_3rd_derivatives[i].resize (n_q_points);
+    }
+
+  if (flags & update_4th_derivatives)
+    {
+      nth_derivatives[4].resize (this->dofs_per_cell);
+      data->shape_4th_derivatives.resize (this->dofs_per_cell);
+      for (unsigned int i=0;i<this->dofs_per_cell;++i)
+        data->shape_4th_derivatives[i].resize (n_q_points);
+    }
+
+  if (flags & update_5th_derivatives)
+    {
+      nth_derivatives[5].resize (this->dofs_per_cell);
+      data->shape_5th_derivatives.resize (this->dofs_per_cell);
+      for (unsigned int i=0;i<this->dofs_per_cell;++i)
+        data->shape_5th_derivatives[i].resize (n_q_points);
+    }
+
+  if (flags & update_4th_derivatives)
+    {
+      nth_derivatives[4].resize (this->dofs_per_cell);
+      data->shape_4th_derivatives.resize (this->dofs_per_cell);
+      for (unsigned int i=0;i<this->dofs_per_cell;++i)
+        data->shape_4th_derivatives[i].resize (n_q_points);
+    }
+
+  if (flags & update_5th_derivatives)
+    {
+      nth_derivatives[5].resize (this->dofs_per_cell);
+      data->shape_5th_derivatives.resize (this->dofs_per_cell);
+      for (unsigned int i=0;i<this->dofs_per_cell;++i)
+        data->shape_5th_derivatives[i].resize (n_q_points);
+    }
+
+  if (flags & update_6th_derivatives)
+    {
+      nth_derivatives[6].resize (this->dofs_per_cell);
+      data->shape_6th_derivatives.resize (this->dofs_per_cell);
+      for (unsigned int i=0;i<this->dofs_per_cell;++i)
+        data->shape_6th_derivatives[i].resize (n_q_points);
+    }
+
+  if (flags & update_7th_derivatives)
+    {
+      nth_derivatives[7].resize (this->dofs_per_cell);
+      data->shape_7th_derivatives.resize (this->dofs_per_cell);
+      for (unsigned int i=0;i<this->dofs_per_cell;++i)
+        data->shape_7th_derivatives[i].resize (n_q_points);
+    }
+
+  if (flags & update_8th_derivatives)
+    {
+      nth_derivatives[8].resize (this->dofs_per_cell);
+      data->shape_8th_derivatives.resize (this->dofs_per_cell);
+      for (unsigned int i=0;i<this->dofs_per_cell;++i)
+        data->shape_8th_derivatives[i].resize (n_q_points);
+    }
+
+  if (flags & update_9th_derivatives)
+    {
+      nth_derivatives[9].resize (this->dofs_per_cell);
+      data->shape_9th_derivatives.resize (this->dofs_per_cell);
+      for (unsigned int i=0;i<this->dofs_per_cell;++i)
+        data->shape_9th_derivatives[i].resize (n_q_points);
     }
 
                                    // Compute shape function values
@@ -344,41 +435,61 @@ FE_PolyTensor<POLY,dim,spacedim>::get_data (
                                    // functions v_j
   if (flags & (update_values | update_gradients))
     for (unsigned int k=0; k<n_q_points; ++k)
-      {
-        poly_space.compute(quadrature.point(k),
-                           values, grads, hessians);
+    {
+      poly_space.compute(quadrature.point(k),
+          values, grads, hessians);
 
-        if (flags & update_values)
+      if (flags & update_values)
+      {
+        if (inverse_node_matrix.n_cols() == 0)
+          for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+            data->shape_values[i][k] = values[i];
+        else
+          for (unsigned int i=0; i<this->dofs_per_cell; ++i)
           {
-            if (inverse_node_matrix.n_cols() == 0)
-              for (unsigned int i=0; i<this->dofs_per_cell; ++i)
-                data->shape_values[i][k] = values[i];
-            else
-              for (unsigned int i=0; i<this->dofs_per_cell; ++i)
-                {
-                  Tensor<1,dim> add_values;
-                  for (unsigned int j=0; j<this->dofs_per_cell; ++j)
-                    add_values += inverse_node_matrix(j,i) * values[j];
-                  data->shape_values[i][k] = add_values;
-                }
+            Tensor<1,dim> add_values;
+            for (unsigned int j=0; j<this->dofs_per_cell; ++j)
+              add_values += inverse_node_matrix(j,i) * values[j];
+            data->shape_values[i][k] = add_values;
           }
+      }
 
-        if (flags & update_gradients)
+      if (flags & update_gradients)
+      {
+        if (inverse_node_matrix.n_cols() == 0)
+          for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+            data->shape_grads[i][k] = grads[i];
+        else
+          for (unsigned int i=0; i<this->dofs_per_cell; ++i)
           {
-            if (inverse_node_matrix.n_cols() == 0)
-              for (unsigned int i=0; i<this->dofs_per_cell; ++i)
-                data->shape_grads[i][k] = grads[i];
-            else
-              for (unsigned int i=0; i<this->dofs_per_cell; ++i)
-                {
-                  Tensor<2,dim> add_grads;
-                  for (unsigned int j=0; j<this->dofs_per_cell; ++j)
-                    add_grads += inverse_node_matrix(j,i) * grads[j];
-                  data->shape_grads[i][k] = add_grads;
-                }
+            Tensor<2,dim> add_grads;
+            for (unsigned int j=0; j<this->dofs_per_cell; ++j)
+              add_grads += inverse_node_matrix(j,i) * grads[j];
+            data->shape_grads[i][k] = add_grads;
           }
-
       }
+    }
+
+  //new version BJ
+  if (flags & update_hessians)
+    for (unsigned int k=0; k<n_q_points; ++k)
+    {
+      poly_space.compute(quadrature.point(k),
+          values, grads, hessians);
+
+      if (inverse_node_matrix.n_cols() == 0)
+        for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+          data->shape_values[i][k] = values[i];
+      else
+        for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+        {
+          Tensor<1,dim> add_values;
+          for (unsigned int j=0; j<this->dofs_per_cell; ++j)
+            add_values += inverse_node_matrix(j,i) * values[j];
+          data->shape_values[i][k] = add_values;
+        }
+    }
+
   return data;
 }
 
@@ -435,7 +546,6 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values (
   for (unsigned int i=0; i<this->dofs_per_cell; ++i)
     {
       const unsigned int first = data.shape_function_to_row_table[i];
-
       if (flags & update_values && cell_similarity != CellSimilarity::translation)
         switch (mapping_type)
           {
index a839d72159e169ad77faea960d4ae8ca0bd976d0..281af89ddb200a5dc2123b21ef4bdd4dc10fb909 100644 (file)
@@ -1737,16 +1737,90 @@ FEValuesData<dim,spacedim>::initialize (const unsigned int        n_quadrature_p
                                    // that we will need to their
                                    // correct size
   if (flags & update_values)
+  {
     this->shape_values.reinit(n_nonzero_shape_components,
                               n_quadrature_points);
+    this->shape_nth_derivatives.resize(1);
+    this->shape_nth_derivatives[0].resize(n_nonzero_shape_components,
+                              std::vector<boost::any> (n_quadrature_points));
+  }
 
   if (flags & update_gradients)
+  {
     this->shape_gradients.resize (n_nonzero_shape_components,
                                   std::vector<Tensor<1,spacedim> > (n_quadrature_points));
+    this->shape_nth_derivatives.resize(2);
+    this->shape_nth_derivatives[1].resize (n_nonzero_shape_components,
+                                  std::vector<boost::any> (n_quadrature_points));
+  }
 
   if (flags & update_hessians)
+  {
     this->shape_hessians.resize (n_nonzero_shape_components,
                                  std::vector<Tensor<2,spacedim> > (n_quadrature_points));
+    this->shape_nth_derivatives.resize(3);
+    this->shape_nth_derivatives[2].resize (n_nonzero_shape_components,
+                                  std::vector<boost::any> (n_quadrature_points));
+  }
+
+  if (flags & update_derivatives(3,9))
+    this->shape_nth_derivatives.resize(10);
+
+  if (flags & update_3rd_derivatives)
+  {
+    this->shape_3rd_derivatives.resize (n_nonzero_shape_components,
+                                 std::vector<Tensor<3,spacedim> > (n_quadrature_points));
+    this->shape_nth_derivatives[3].resize(n_nonzero_shape_components, 
+                                 std::vector<boost::any> (n_quadrature_points));
+  }
+
+  if (flags & update_4th_derivatives)
+  {
+    this->shape_4th_derivatives.resize (n_nonzero_shape_components,
+                                 std::vector<Tensor<4,spacedim> > (n_quadrature_points));
+    this->shape_nth_derivatives[4].resize(n_nonzero_shape_components, 
+                                 std::vector<boost::any> (n_quadrature_points));
+  }
+
+  if (flags & update_5th_derivatives)
+  {
+    this->shape_5th_derivatives.resize (n_nonzero_shape_components,
+                                 std::vector<Tensor<5,spacedim> > (n_quadrature_points));
+    this->shape_nth_derivatives[5].resize(n_nonzero_shape_components, 
+                                 std::vector<boost::any> (n_quadrature_points));
+  }
+
+  if (flags & update_6th_derivatives)
+  {
+    this->shape_6th_derivatives.resize (n_nonzero_shape_components,
+                                 std::vector<Tensor<6,spacedim> > (n_quadrature_points));
+    this->shape_nth_derivatives[6].resize(n_nonzero_shape_components, 
+                                 std::vector<boost::any> (n_quadrature_points));
+  }
+
+  if (flags & update_7th_derivatives)
+  {
+    this->shape_7th_derivatives.resize (n_nonzero_shape_components,
+                                 std::vector<Tensor<7,spacedim> > (n_quadrature_points));
+    this->shape_nth_derivatives[7].resize(n_nonzero_shape_components, 
+                                 std::vector<boost::any> (n_quadrature_points));
+  }
+
+  if (flags & update_8th_derivatives)
+  {
+    this->shape_8th_derivatives.resize (n_nonzero_shape_components,
+                                 std::vector<Tensor<8,spacedim> > (n_quadrature_points));
+    this->shape_nth_derivatives[8].resize(n_nonzero_shape_components, 
+                                 std::vector<boost::any> (n_quadrature_points));
+  }
+
+  if (flags & update_9th_derivatives)
+  {
+    this->shape_9th_derivatives.resize (n_nonzero_shape_components,
+                                 std::vector<Tensor<9,spacedim> > (n_quadrature_points));
+    this->shape_nth_derivatives[9].resize(n_nonzero_shape_components, 
+                                 std::vector<boost::any> (n_quadrature_points));
+  }
 
   if (flags & update_quadrature_points)
     this->quadrature_points.resize(n_quadrature_points);

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.