]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
extend DerivativeApproximation to third order derivatives, enable returning the full...
authorleicht <leicht@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 4 Aug 2006 05:02:02 +0000 (05:02 +0000)
committerleicht <leicht@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 4 Aug 2006 05:02:02 +0000 (05:02 +0000)
git-svn-id: https://svn.dealii.org/trunk@13591 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/derivative_approximation.h
deal.II/deal.II/source/numerics/derivative_approximation.cc

index 291db5f75f95a49feb281db8b032139162160db6..62098de3d6a9c7f3cfc93cafd1dfd02834b5a8f4 100644 (file)
@@ -1,3 +1,4 @@
+
 //---------------------------------------------------------------------------
 //    $Id$
 //    Version: $Name$
@@ -265,6 +266,53 @@ class DerivativeApproximation
                                   Vector<float>         &derivative_norm,
                                   const unsigned int     component = 0);
 
+                                    /**
+                                     * This function calculates the
+                                     * <tt>order</tt>-th order approximate
+                                     * derivative and returns the full tensor
+                                     * for a single cell.
+                                     *
+                                     * The last parameter denotes the
+                                     * solution component, for which
+                                     * the gradient is to be
+                                     * computed. It defaults to the
+                                     * first component. For
+                                     * scalar elements, this is the only
+                                     * valid choice; for vector-valued ones,
+                                     * any component between zero and the
+                                     * number of vector components can be
+                                     * given here.
+                                     */
+
+    template <int dim, template <int> class DH, class InputVector, int order>
+    static void
+    approximate_derivative_tensor (const Mapping<dim>                           &mapping,
+                                  const DH<dim>                                &dof,
+                                  const InputVector                            &solution,
+                                  const typename DH<dim>::active_cell_iterator &cell,
+                                  Tensor<order,dim>                            &derivative,
+                                  const unsigned int                            component = 0);
+
+                                    /**
+                                     * Same as above, with
+                                     * <tt>mapping=MappingQ1@<dim@>()</tt>.
+                                     */
+
+    template <int dim, template <int> class DH, class InputVector, int order>
+    static void
+    approximate_derivative_tensor (const DH<dim>                                &dof,
+                                  const InputVector                            &solution,
+                                  const typename DH<dim>::active_cell_iterator &cell,
+                                  Tensor<order,dim>                            &derivative,
+                                  const unsigned int                            component = 0);
+
+                                    /**
+                                     * Return the norm of the derivative.
+                                     */
+    template <int dim, int order>
+    static double
+    derivative_norm(const Tensor<order,dim> &derivative);
+    
                                     /**
                                      * Exception
                                      */
@@ -449,6 +497,130 @@ class DerivativeApproximation
                                          */
        static void symmetrize (Derivative &derivative_tensor);
     };
+
+    template <int dim>
+    class ThirdDerivative
+    {
+      public:
+                                        /**
+                                         * Declare which data fields have
+                                         * to be updated for the function
+                                         * @p get_projected_derivative
+                                         * to work.
+                                         */
+       static const UpdateFlags update_flags;
+
+                                        /**
+                                         * Declare the data type which
+                                         * holds the derivative described
+                                         * by this class.
+                                         */
+       typedef Tensor<3,dim> Derivative;
+
+                                        /**
+                                         * Likewise declare the data type
+                                         * that holds the derivative
+                                         * projected to a certain
+                                         * directions.
+                                         */
+       typedef Tensor<2,dim> ProjectedDerivative;
+
+                                        /**
+                                         * Given an FEValues object
+                                         * initialized to a cell, and a
+                                         * solution vector, extract the
+                                         * desired derivative at the
+                                         * first quadrature point (which
+                                         * is the only one, as we only
+                                         * evaluate the finite element
+                                         * field at the center of each
+                                         * cell).
+                                         */
+       template <class InputVector>
+       static ProjectedDerivative
+       get_projected_derivative (const FEValues<dim>  &fe_values,
+                                 const InputVector    &solution,
+                                 const unsigned int    component);
+       
+                                        /**
+                                         * Return the norm of the
+                                         * derivative object. Here, for
+                                         * the (symmetric) tensor of
+                                         * second derivatives, we choose
+                                         * the absolute value of the
+                                         * largest eigenvalue, which is
+                                         * the matrix norm associated to
+                                         * the $l_2$ norm of vectors. It
+                                         * is also the largest value of
+                                         * the curvature of the solution.
+                                         */
+       static double derivative_norm (const Derivative &d);
+
+                                        /**
+                                         * If for the present derivative
+                                         * order, symmetrization of the
+                                         * derivative tensor is
+                                         * necessary, then do so on the
+                                         * argument.
+                                         *
+                                         * For the second derivatives,
+                                         * each entry of the tensor is
+                                         * set to the mean of its value
+                                         * and the value of the transpose
+                                         * element.
+                                         *
+                                         * Note that this function
+                                         * actually modifies its
+                                         * argument.
+                                         */
+       static void symmetrize (Derivative &derivative_tensor);
+    };
+    
+    template <int order, int dim>
+    class DerivativeSelector
+    {
+      public:
+                                        /**
+                                         * typedef to select the
+                                         * DerivativeDescription corresponding
+                                         * to the <tt>order</tt>th
+                                         * derivative. In this general template
+                                         * we set an unvalid typedef to void,
+                                         * the real typedefs have to be
+                                         * specialized.
+                                         */
+       typedef void DerivDescr;
+       
+    };
+
+    template <int dim>
+    class DerivativeSelector<1,dim>
+    {
+      public:
+
+       typedef Gradient<dim> DerivDescr;
+    };
+    
+    template <int dim>
+    class DerivativeSelector<2,dim>
+    {
+      public:
+       
+       typedef SecondDerivative<dim> DerivDescr;
+    };
+    
+    template <int dim>
+    class DerivativeSelector<3,dim>
+    {
+      public:
+       
+       typedef ThirdDerivative<dim> DerivDescr;
+    };
+    
+    
+    
+    
+  private:
     
                                     /**
                                      * Convenience typedef denoting
@@ -491,6 +663,9 @@ class DerivativeApproximation
                                      * approximation on the cells in
                                      * the range given by the third
                                      * parameter.
+                                     * Fill the @p derivative_norm vector with
+                                     * the norm of the computed derivative
+                                     * tensors on each cell.
                                      */
     template <class DerivativeDescription, int dim,
               template <int> class DH, class InputVector>
@@ -501,6 +676,21 @@ class DerivativeApproximation
                 const unsigned int     component,
                 const IndexInterval   &index_interval,
                 Vector<float>         &derivative_norm);    
+
+                                    /**
+                                     * Compute the derivative approximation on
+                                     * one cell. This computes the full
+                                     * derivative tensor.
+                                     */
+    template <class DerivativeDescription, int dim,
+              template <int> class DH, class InputVector>
+    static void
+    approximate_cell (const Mapping<dim>                            &mapping,
+                     const DH<dim>                                 &dof,
+                     const InputVector                             &solution,
+                     const unsigned int                             component,
+                     const typename DH<dim>::active_cell_iterator  &cell,
+                     typename DerivativeDescription::Derivative    &derivative);    
 };
 
 
index f33de9df38a5d9767173db905274a435be2ec98f..8bffc3207f33e3d87131fb9026cff50bb2e6b288 100644 (file)
@@ -50,6 +50,9 @@ const UpdateFlags DerivativeApproximation::Gradient<dim>::update_flags = update_
 template <int dim>
 const UpdateFlags DerivativeApproximation::SecondDerivative<dim>::update_flags = update_gradients;
 
+template <int dim>
+const UpdateFlags DerivativeApproximation::ThirdDerivative<dim>::update_flags = update_second_derivatives;
+
 
 
 template <int dim>
@@ -388,6 +391,113 @@ DerivativeApproximation::SecondDerivative<dim>::symmetrize (Derivative &d)
 }
 
 
+template <int dim>
+template <class InputVector>
+inline
+typename DerivativeApproximation::ThirdDerivative<dim>::ProjectedDerivative
+DerivativeApproximation::ThirdDerivative<dim>::
+get_projected_derivative (const FEValues<dim>  &fe_values,
+                         const InputVector    &solution,
+                         const unsigned int    component)
+{
+  if (fe_values.get_fe().n_components() == 1)
+    {
+      std::vector<ProjectedDerivative> values (1);
+      fe_values.get_function_2nd_derivatives (solution, values);
+      return values[0];
+    }
+  else
+    {
+      std::vector<std::vector<ProjectedDerivative> > values
+       (1, std::vector<ProjectedDerivative>(fe_values.get_fe().n_components()));
+      fe_values.get_function_2nd_derivatives (solution, values);
+      return values[0][component];
+    };
+}
+
+
+#if deal_II_dimension == 1
+
+template <>
+inline
+double
+DerivativeApproximation::ThirdDerivative<1>::
+derivative_norm (const Derivative &d)
+{
+  return std::fabs (d[0][0][0]);
+}
+
+#endif
+
+
+template <int dim>
+inline
+double
+DerivativeApproximation::ThirdDerivative<dim>::
+derivative_norm (const Derivative &d)
+{
+                                  // return the Frobenius-norm. this is a
+                                  // member function of Tensor<rank_,dim>
+  return d.norm();
+}
+
+
+template <int dim>
+inline
+void
+DerivativeApproximation::ThirdDerivative<dim>::symmetrize (Derivative &d)
+{
+                                  // symmetrize non-diagonal entries
+
+                                  // first do it in the case, that i,j,k are
+                                  // pairwise different (which can onlky happen
+                                  // in dim >= 3)
+  for (unsigned int i=0; i<dim; ++i)
+    for (unsigned int j=i+1; j<dim; ++j)
+      for (unsigned int k=j+1; k<dim; ++k)
+       {
+         const double s = (d[i][j][k] +
+                           d[i][k][j] +
+                           d[j][i][k] +
+                           d[j][k][i] +
+                           d[k][i][j] +
+                           d[k][j][i]) / 6;
+         d[i][j][k]
+           = d[i][k][j]
+           = d[j][i][k]
+           = d[j][k][i]
+           = d[k][i][j]
+           = d[k][j][i]
+           = s;
+       };
+                                  // now do the case, where two indices are
+                                  // equal
+  for (unsigned int i=0; i<dim; ++i)
+    for (unsigned int j=i+1; j<dim; ++j)
+      {
+                                        // case 1: index i (lower one) is
+                                        // double
+       const double s = (d[i][i][j] +
+                         d[i][j][i] +
+                         d[j][i][i] ) / 3;
+       d[i][i][j]
+         = d[i][j][i]
+         = d[j][i][i]
+         = s;
+
+                                        // case 2: index j (higher one) is
+                                        // double
+       const double t = (d[i][j][j] +
+                         d[j][i][j] +
+                         d[j][j][i] ) / 3;
+       d[i][j][j]
+         = d[j][i][j]
+         = d[j][j][i]
+         = t;
+      };
+  
+}
+
 
 
 template <int dim, template <int> class DH, class InputVector>
@@ -458,6 +568,49 @@ approximate_second_derivative (const DH<dim>         &dof_handler,
 }
 
 
+template <int dim, template <int> class DH, class InputVector, int order>
+void
+DerivativeApproximation::
+approximate_derivative_tensor (const Mapping<dim>                           &mapping,
+                              const DH<dim>                                &dof,
+                              const InputVector                            &solution,
+                              const typename DH<dim>::active_cell_iterator &cell,
+                              Tensor<order,dim>                            &derivative,
+                              const unsigned int                            component)
+{
+  approximate_cell<typename DerivativeSelector<order,dim>::DerivDescr,dim,DH,InputVector>
+    (mapping,
+     dof,
+     solution,
+     component,
+     cell,
+     derivative);
+}
+
+
+
+template <int dim, template <int> class DH, class InputVector, int order>
+void
+DerivativeApproximation::
+approximate_derivative_tensor (const DH<dim>                                &dof,
+                              const InputVector                            &solution,
+                              const typename DH<dim>::active_cell_iterator &cell,
+                              Tensor<order,dim>                            &derivative,
+                              const unsigned int                            component)
+{
+                                  // just call the respective function with Q1
+                                  // mapping
+  Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+  approximate_derivative_tensor (StaticMappingQ1<dim>::mapping,
+                                dof,
+                                solution,
+                                cell,
+                                derivative,
+                                component);
+}
+
+
+
 namespace WorkAround
 {
 // gcc 2.95 is not happy if we take the address of a template function as in
@@ -538,6 +691,51 @@ DerivativeApproximation::approximate (const Mapping<dim>    &mapping,
                                      const IndexInterval   &index_interval,
                                      Vector<float>         &derivative_norm)
 {
+                                  // iterators over all cells and the
+                                  // respective entries in the output
+                                  // vector:
+  Vector<float>::iterator
+    derivative_norm_on_this_cell
+    = derivative_norm.begin() + index_interval.first;
+  
+  typename DH<dim>::active_cell_iterator cell, endc;
+  cell = endc = dof_handler.begin_active();
+                                  // (static_cast to avoid warnings
+                                  // about unsigned always >=0)
+  std::advance (cell, static_cast<int>(index_interval.first));
+  std::advance (endc, static_cast<int>(index_interval.second));
+
+  for (; cell!=endc; ++cell, ++derivative_norm_on_this_cell)
+    {
+      typename DerivativeDescription::Derivative derivative;
+                                      // call the function doing the actual
+                                      // work on this cell
+      DerivativeApproximation::
+       template approximate_cell<DerivativeDescription,dim,DH,InputVector>
+       (mapping,
+        dof_handler,
+        solution,
+        component,
+        cell,
+        derivative);
+                                      // evaluate the norm and fill the vector
+      *derivative_norm_on_this_cell
+       = DerivativeDescription::derivative_norm (derivative);
+    }
+}
+
+
+template <class DerivativeDescription, int dim,
+          template <int> class DH, class InputVector>
+void 
+DerivativeApproximation::
+approximate_cell (const Mapping<dim>                            &mapping,
+                 const DH<dim>                                 &dof_handler,
+                 const InputVector                             &solution,
+                 const unsigned int                             component,
+                 const typename DH<dim>::active_cell_iterator  &cell,
+                 typename DerivativeDescription::Derivative    &derivative)    
+{ 
   QMidpoint<dim> midpoint_rule;
 
                                   // create collection objects from
@@ -559,19 +757,6 @@ DerivativeApproximation::approximate (const Mapping<dim>    &mapping,
                                   // matrix Y=sum_i y_i y_i^T
   Tensor<2,dim> Y;
   
-                                  // iterators over all cells and the
-                                  // respective entries in the output
-                                  // vector:
-  Vector<float>::iterator
-    derivative_norm_on_this_cell
-    = derivative_norm.begin() + index_interval.first;
-  
-  typename DH<dim>::active_cell_iterator cell, endc;
-  cell = endc = dof_handler.begin_active();
-                                  // (static_cast to avoid warnings
-                                  // about unsigned always >=0)
-  std::advance (cell, static_cast<int>(index_interval.first));
-  std::advance (endc, static_cast<int>(index_interval.second));
 
                                   // vector to hold iterators to all
                                   // active neighbors of a cell
@@ -581,9 +766,6 @@ DerivativeApproximation::approximate (const Mapping<dim>    &mapping,
   active_neighbors.reserve (GeometryInfo<dim>::faces_per_cell *
                            GeometryInfo<dim>::subfaces_per_face);
 
-  for (; cell!=endc; ++cell, ++derivative_norm_on_this_cell)
-    {
-      Y.clear ();
                                       // vector
                                       // g=sum_i y_i (f(x+y_i)-f(x))/|y_i|
                                       // or related type for higher
@@ -665,7 +847,7 @@ DerivativeApproximation::approximate (const Mapping<dim>    &mapping,
                  for (unsigned int c=0; c<neighbor->n_children(); ++c)
                    for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
                      if (neighbor->child(c)->neighbor(f) == cell)
-                       active_neighbors.push_back (neighbor->child(c));
+                         active_neighbors.push_back (neighbor->child(c));
              };
          };
 
@@ -751,60 +933,95 @@ DerivativeApproximation::approximate (const Mapping<dim>    &mapping,
                                        // compute Y^-1 g
       const Tensor<2,dim> Y_inverse = invert(Y);
       
-      typename DerivativeDescription::Derivative derivative;
       contract (derivative, Y_inverse, projected_derivative);
       
                                       // finally symmetrize the derivative
       DerivativeDescription::symmetrize (derivative);
+}
 
-      *derivative_norm_on_this_cell
-       = DerivativeDescription::derivative_norm (derivative);
-    }
+
+template <int dim, int order>
+double
+DerivativeApproximation::
+derivative_norm(const Tensor<order,dim> &derivative)
+{
+  return DerivativeSelector<order,dim>::DerivDescr::derivative_norm(derivative);
 }
 
 
 // --------------------------------------------------------------------
 
 // explicit instantiations
-#define INSTANTIATE(InputVector,DH)                \
-template                                           \
-void                                               \
-DerivativeApproximation::                          \
-approximate_gradient<deal_II_dimension>            \
-(const Mapping<deal_II_dimension> &mapping,        \
- const DH<deal_II_dimension> &dof_handler,         \
- const InputVector  &solution,                     \
- Vector<float>         &derivative_norm,           \
- const unsigned int     component);                \
-                                                   \
-template                                           \
-void                                               \
-DerivativeApproximation::                          \
-approximate_gradient<deal_II_dimension>            \
-(const DH<deal_II_dimension> &dof_handler,         \
- const InputVector     &solution,                  \
- Vector<float>         &derivative_norm,           \
- const unsigned int     component);                \
-                                                   \
-template                                           \
-void                                               \
-DerivativeApproximation::                          \
-approximate_second_derivative<deal_II_dimension>   \
-(const Mapping<deal_II_dimension> &mapping,        \
- const DH<deal_II_dimension> &dof_handler,         \
- const InputVector  &solution,                     \
- Vector<float>         &derivative_norm,           \
- const unsigned int     component);                \
-                                                   \
-template                                           \
-void                                               \
-DerivativeApproximation::                          \
-approximate_second_derivative<deal_II_dimension>   \
-(const DH<deal_II_dimension> &dof_handler,         \
- const InputVector     &solution,                  \
- Vector<float>         &derivative_norm,           \
+#define INSTANTIATE(InputVector,DH)                      \
+template                                                 \
+void                                                     \
+DerivativeApproximation::                                \
+approximate_gradient<deal_II_dimension>                  \
+(const Mapping<deal_II_dimension> &mapping,              \
+ const DH<deal_II_dimension> &dof_handler,               \
+ const InputVector  &solution,                           \
+ Vector<float>         &derivative_norm,                 \
+ const unsigned int     component);                      \
+                                                         \
+template                                                 \
+void                                                     \
+DerivativeApproximation::                                \
+approximate_gradient<deal_II_dimension>                  \
+(const DH<deal_II_dimension> &dof_handler,               \
+ const InputVector     &solution,                        \
+ Vector<float>         &derivative_norm,                 \
+ const unsigned int     component);                      \
+                                                         \
+template                                                 \
+void                                                     \
+DerivativeApproximation::                                \
+approximate_second_derivative<deal_II_dimension>         \
+(const Mapping<deal_II_dimension> &mapping,              \
+ const DH<deal_II_dimension> &dof_handler,               \
+ const InputVector  &solution,                           \
+ Vector<float>         &derivative_norm,                 \
+ const unsigned int     component);                      \
+                                                         \
+template                                                 \
+void                                                     \
+DerivativeApproximation::                                \
+approximate_second_derivative<deal_II_dimension>         \
+(const DH<deal_II_dimension> &dof_handler,               \
+ const InputVector     &solution,                        \
+ Vector<float>         &derivative_norm,                 \
+ const unsigned int     component);                      \
+                                                         \
+template                                                 \
+void                                                     \
+DerivativeApproximation::                                \
+approximate_derivative_tensor<deal_II_dimension>         \
+(const DH<deal_II_dimension> &dof_handler,               \
+ const InputVector     &solution,                        \
+ const DH<deal_II_dimension>::active_cell_iterator &cell,\
+ Tensor<1,deal_II_dimension> &derivative,                \
+ const unsigned int     component);                      \
+                                                         \
+template                                                 \
+void                                                     \
+DerivativeApproximation::                                \
+approximate_derivative_tensor<deal_II_dimension>         \
+(const DH<deal_II_dimension> &dof_handler,               \
+ const InputVector     &solution,                        \
+ const DH<deal_II_dimension>::active_cell_iterator &cell,\
+ Tensor<2,deal_II_dimension> &derivative,                \
+ const unsigned int     component);                      \
+                                                         \
+template                                                 \
+void                                                     \
+DerivativeApproximation::                                \
+approximate_derivative_tensor<deal_II_dimension>         \
+(const DH<deal_II_dimension> &dof_handler,               \
+ const InputVector     &solution,                        \
+ const DH<deal_II_dimension>::active_cell_iterator &cell,\
+ Tensor<3,deal_II_dimension> &derivative,                \
  const unsigned int     component)
 
+
 INSTANTIATE(Vector<double>, DoFHandler);
 INSTANTIATE(Vector<float>, DoFHandler);
 INSTANTIATE(BlockVector<double>, DoFHandler);
@@ -823,6 +1040,22 @@ INSTANTIATE(PETScWrappers::Vector, hp::DoFHandler);
 INSTANTIATE(PETScWrappers::BlockVector, hp::DoFHandler);
 #endif
 
+template
+double
+DerivativeApproximation::
+derivative_norm(const Tensor<1,deal_II_dimension> &derivative);
+
+template
+double
+DerivativeApproximation::
+derivative_norm(const Tensor<2,deal_II_dimension> &derivative);
+
+template
+double
+DerivativeApproximation::
+derivative_norm(const Tensor<3,deal_II_dimension> &derivative);
+
+
 
 // static variables
 // 
@@ -836,3 +1069,7 @@ DerivativeApproximation::Gradient<deal_II_dimension>::update_flags;
 template
 const UpdateFlags
 DerivativeApproximation::SecondDerivative<deal_II_dimension>::update_flags;
+
+template
+const UpdateFlags
+DerivativeApproximation::ThirdDerivative<deal_II_dimension>::update_flags;

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.