From a325b28ab255279c95d12a616254b2830ff65de5 Mon Sep 17 00:00:00 2001
From: Maien Hamed <tomaien@hotmail.com>
Date: Mon, 7 Sep 2015 18:03:48 +0200
Subject: [PATCH] removed compute_2nd as hessians are analytically computed

---
 include/deal.II/fe/fe.h |  18 ---
 source/fe/fe.cc         | 242 ----------------------------------------
 2 files changed, 260 deletions(-)

diff --git a/include/deal.II/fe/fe.h b/include/deal.II/fe/fe.h
index d24f96272a..ca8ee1035a 100644
--- a/include/deal.II/fe/fe.h
+++ b/include/deal.II/fe/fe.h
@@ -519,14 +519,6 @@ public:
      */
     virtual ~InternalDataBase ();
 
-    /**
-     * Initialize some pointers used in the computation of second derivatives
-     * by finite differencing of gradients.
-     */
-    void initialize_2nd (const FiniteElement<dim,spacedim> *element,
-                         const Mapping<dim,spacedim>       &mapping,
-                         const Quadrature<dim>    &quadrature);
-
     /**
      * Values updated by the constructor or by reinit.
      */
@@ -2213,16 +2205,6 @@ protected:
   TableIndices<2>
   interface_constraints_size () const;
 
-  /**
-   * Compute second derivatives by finite differences of gradients.
-   */
-  void compute_2nd (const Mapping<dim,spacedim>                                         &mapping,
-                    const typename Triangulation<dim,spacedim>::cell_iterator           &cell,
-                    const unsigned int                                                   offset,
-                    const typename Mapping<dim,spacedim>::InternalDataBase              &mapping_internal,
-                    const InternalDataBase                                              &fe_internal,
-                    dealii::internal::FEValues::FiniteElementRelatedData<dim, spacedim> &data) const;
-
   /**
    * Given the pattern of nonzero components for each shape function, compute
    * for each entry how many components are non-zero for each shape function.
diff --git a/source/fe/fe.cc b/source/fe/fe.cc
index 28d06a68f6..291eabc240 100644
--- a/source/fe/fe.cc
+++ b/source/fe/fe.cc
@@ -73,58 +73,6 @@ FiniteElement<dim, spacedim>::InternalDataBase::memory_consumption () const
 
 
 
-
-template <int dim, int spacedim>
-void
-FiniteElement<dim,spacedim>::
-InternalDataBase::initialize_2nd (const FiniteElement<dim,spacedim> *element,
-                                  const Mapping<dim,spacedim>       &mapping,
-                                  const Quadrature<dim>    &quadrature)
-{
-  // if we shall compute second
-  // derivatives, then we do so by
-  // finite differencing the
-  // gradients. that we do by
-  // evaluating the gradients of
-  // shape values at points shifted
-  // star-like a little in each
-  // coordinate direction around each
-  // quadrature point.
-  //
-  // therefore generate 2*dim (the
-  // number of evaluation points)
-  // FEValues objects with slightly
-  // shifted positions
-  std::vector<Point<dim> > diff_points (quadrature.size());
-
-  differences.resize(2*dim);
-  for (unsigned int d=0; d<dim; ++d)
-    {
-      Tensor<1,dim> shift;
-      shift[d] = fd_step_length;
-
-      // generate points and FEValues
-      // objects shifted in
-      // plus-direction. note that
-      // they only need to compute
-      // gradients, not more
-      for (unsigned int i=0; i<diff_points.size(); ++i)
-        diff_points[i] = quadrature.point(i) + shift;
-      const Quadrature<dim> plus_quad (diff_points, quadrature.get_weights());
-      differences[d] = new FEValues<dim,spacedim> (mapping, *element,
-                                                   plus_quad, update_gradients);
-
-      // now same in minus-direction
-      for (unsigned int i=0; i<diff_points.size(); ++i)
-        diff_points[i] = quadrature.point(i) - shift;
-      const Quadrature<dim> minus_quad (diff_points, quadrature.get_weights());
-      differences[d+dim] = new FEValues<dim,spacedim> (mapping, *element,
-                                                       minus_quad, update_gradients);
-    }
-}
-
-
-
 template <int dim, int spacedim>
 UpdateFlags
 FiniteElement<dim,spacedim>::InternalDataBase::current_update_flags () const
@@ -1283,196 +1231,6 @@ FiniteElement<dim,spacedim>::memory_consumption () const
 
 
 
-template<>
-void
-FiniteElement<1,2>::compute_2nd (
-  const Mapping<1,2> &,
-  const Triangulation<1,2>::cell_iterator &,
-  const unsigned int,
-  const Mapping<1,2>::InternalDataBase &,
-  const InternalDataBase &,
-  internal::FEValues::FiniteElementRelatedData<1,2> &) const
-{
-  Assert(false, ExcNotImplemented());
-}
-
-
-template<>
-void
-FiniteElement<1,3>::compute_2nd (
-  const Mapping<1,3> &,
-  const Triangulation<1,3>::cell_iterator &,
-  const unsigned int,
-  const Mapping<1,3>::InternalDataBase &,
-  const InternalDataBase &,
-  internal::FEValues::FiniteElementRelatedData<1,3> &) const
-{
-  Assert(false, ExcNotImplemented());
-}
-
-
-
-template<>
-void
-FiniteElement<2,3>::compute_2nd (
-  const Mapping<2,3> &,
-  const Triangulation<2,3>::cell_iterator &,
-  const unsigned int,
-  const Mapping<2,3>::InternalDataBase &,
-  const InternalDataBase &,
-  internal::FEValues::FiniteElementRelatedData<2,3> &) const
-{
-  Assert(false, ExcNotImplemented());
-}
-
-
-
-template <int dim, int spacedim>
-void
-FiniteElement<dim,spacedim>::compute_2nd (
-  const Mapping<dim,spacedim>                   &mapping,
-  const typename Triangulation<dim,spacedim>::cell_iterator &cell,
-  const unsigned int offset,
-  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_internal,
-  const InternalDataBase                     &fe_internal,
-  internal::FEValues::FiniteElementRelatedData<dim,spacedim>                    &data) const
-{
-  Assert ((fe_internal.update_each | fe_internal.update_once)
-          & update_hessians,
-          ExcInternalError());
-
-  // there is nothing to do if there are no degrees of freedom (e.g., in an
-  // FE_Nothing)
-  if (this->dofs_per_cell == 0)
-    return;
-
-// make sure we have as many entries as there are nonzero components
-//  Assert (data.shape_hessians.size() ==
-//          std::accumulate (n_nonzero_components_table.begin(),
-//                        n_nonzero_components_table.end(),
-//                        0U),
-//        ExcInternalError());
-  // Number of quadrature points
-  const unsigned int n_q_points = data.shape_hessians[0].size();
-
-  // first reinit the fe_values
-  // objects used for the finite
-  // differencing stuff
-  for (unsigned int d=0; d<dim; ++d)
-    {
-      fe_internal.differences[d]->reinit(cell);
-      fe_internal.differences[d+dim]->reinit(cell);
-      Assert(offset <= fe_internal.differences[d]->n_quadrature_points - n_q_points,
-             ExcIndexRange(offset, 0, fe_internal.differences[d]->n_quadrature_points
-                           - n_q_points));
-    }
-
-  // collection of difference
-  // quotients of gradients in each
-  // direction (first index) and at
-  // all q-points (second index)
-  std::vector<std::vector<Tensor<1,dim> > >
-  diff_quot (spacedim, std::vector<Tensor<1,dim> > (n_q_points));
-  std::vector<Tensor<1,spacedim> > diff_quot2 (n_q_points);
-
-  // for all nonzero components of
-  // all shape functions at all
-  // quadrature points and difference
-  // quotients in all directions:
-  unsigned int total_index = 0;
-  for (unsigned int shape_index=0; shape_index<this->dofs_per_cell; ++shape_index)
-    for (unsigned int n=0; n<n_nonzero_components(shape_index); ++n, ++total_index)
-      {
-        for (unsigned int d1=0; d1<dim; ++d1)
-          for (unsigned int q=0; q<n_q_points; ++q)
-            {
-              // get gradient at points
-              // shifted slightly to
-              // the right and to the
-              // left in the present
-              // coordinate direction
-              //
-              // note that things
-              // might be more
-              // difficult if the
-              // shape function has
-              // more than one
-              // non-zero component,
-              // so find out about
-              // the actual component
-              // if necessary
-              Tensor<1,spacedim> right, left;
-              if (is_primitive(shape_index))
-                {
-                  right = fe_internal.differences[d1]->shape_grad(shape_index, q+offset);
-                  left  = fe_internal.differences[d1+dim]->shape_grad(shape_index, q+offset);
-                }
-              else
-                {
-                  // get the
-                  // component index
-                  // of the n-th
-                  // nonzero
-                  // component
-                  unsigned int component=0;
-                  for (unsigned int nonzero_comp=0; component<this->n_components();
-                       ++component)
-                    if (nonzero_components[shape_index][component] == true)
-                      {
-                        ++nonzero_comp;
-                        // check
-                        // whether we
-                        // have found
-                        // the
-                        // component
-                        // we are
-                        // looking
-                        // for. note
-                        // that
-                        // nonzero_comp
-                        // is 1-based
-                        // by the way
-                        // we compute
-                        // it
-                        if (nonzero_comp == n+1)
-                          break;
-                      }
-                  Assert (component < this->n_components(),
-                          ExcInternalError());
-
-                  right = fe_internal.differences[d1]
-                          ->shape_grad_component(shape_index, q+offset, component);
-                  left  = fe_internal.differences[d1+dim]
-                          ->shape_grad_component(shape_index, q+offset, component);
-                };
-
-              // compute the second
-              // derivative from a
-              // symmetric difference
-              // approximation
-              for (unsigned int d=0; d<spacedim; ++d)
-                diff_quot[d][q][d1] = 1./(2*fd_step_length) * (right[d]-left[d]);
-            }
-
-        // up to now we still have
-        // difference quotients on the
-        // unit cell, so transform it
-        // to something on the real
-        // cell
-        for (unsigned int d=0; d<spacedim; ++d)
-          {
-            mapping.transform (diff_quot[d], mapping_covariant, mapping_internal, diff_quot2);
-
-            for (unsigned int q=0; q<n_q_points; ++q)
-              for (unsigned int d1=0; d1<spacedim; ++d1)
-                data.shape_hessians[total_index][q][d][d1]
-                  = diff_quot2[q][d1];
-          }
-      }
-}
-
-
-
 template <int dim, int spacedim>
 std::vector<unsigned int>
 FiniteElement<dim,spacedim>::compute_n_nonzero_components (
-- 
2.39.5