]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Merged from trunk to branch_manifold_id
authorheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 25 Apr 2013 17:27:25 +0000 (17:27 +0000)
committerheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 25 Apr 2013 17:27:25 +0000 (17:27 +0000)
git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@29387 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/cmake/configure/configure_trilinos.cmake
deal.II/include/deal.II/base/logstream.h
deal.II/include/deal.II/matrix_free/fe_evaluation.h
deal.II/source/base/logstream.cc

index 825d304ab33ec44d90113e27083102f2c97c7c25..346d2165fb3c55b3fefc73d69f872f4fee083232 100644 (file)
@@ -26,6 +26,13 @@ MACRO(FEATURE_TRILINOS_FIND_EXTERNAL var)
     #
     SET(${var} TRUE)
 
+    #
+    # Set TRILINOS_DIR to something meaningful if empty
+    #
+    IF("${TRILINOS_DIR}" STREQUAL "")
+      SET(TRILINOS_DIR "<system location>")
+    ENDIF()
+
     #
     # Check whether all required modules of trilinos are installed:
     #
@@ -34,9 +41,9 @@ MACRO(FEATURE_TRILINOS_FIND_EXTERNAL var)
       )
 
     FOREACH(_module
-      amesos epetra ifpack aztecoo sacado teuchos
+      Amesos Epetra Ifpack AztecOO Sacado Teuchos
       )
-      ITEM_MATCHES(_module_found ${_module}$ ${Trilinos_LIBRARIES})
+      ITEM_MATCHES(_module_found ${_module} ${Trilinos_PACKAGE_LIST})
       IF(_module_found)
         MESSAGE(STATUS "Found ${_module}")
       ELSE()
@@ -51,12 +58,10 @@ MACRO(FEATURE_TRILINOS_FIND_EXTERNAL var)
         "Missing ${_modules_missing}"
         )
       SET(TRILINOS_ADDITIONAL_ERROR_STRING
-        "The Trilinos installation found at\n"
-        "  ${TRILINOS_DIR}\n"
+        "The Trilinos installation (found at \"${TRILINOS_DIR}\")\n"
         "is missing one or more modules necessary for the deal.II Trilinos interfaces:\n"
         "  ${_modules_missing}\n\n"
-        "Please re-install Trilinos with the missing Trilinos subpackages
-        enabled.\n\n"
+        "Please re-install Trilinos with the missing Trilinos subpackages enabled.\n\n"
         )
     ENDIF()
 
@@ -82,8 +87,7 @@ MACRO(FEATURE_TRILINOS_FIND_EXTERNAL var)
         )
       SET(TRILINOS_ADDITIONAL_ERROR_STRING
         ${TRILINOS_ADDITIONAL_ERROR_STRING}
-        "The Trilinos installation found at\n"
-        "  ${TRILINOS_DIR}\n"
+        "The Trilinos installation (found at \"${TRILINOS_DIR}\")\n"
         "with version ${TRILINOS_VERSION_MAJOR}.${TRILINOS_VERSION_MINOR}.${TRILINOS_VERSION_SUBMINOR} has bugs that make\n"
         "it incompatible with deal.II. Please use versions before 10.6 or after\n"
         "10.8.1.\n\n"
@@ -103,8 +107,7 @@ MACRO(FEATURE_TRILINOS_FIND_EXTERNAL var)
         )
       SET(TRILINOS_ADDITIONAL_ERROR_STRING
         ${TRILINOS_ADDITIONAL_ERROR_STRING}
-        "The Trilinos installation found at\n"
-        "  ${TRILINOS_DIR}\n"
+        "The Trilinos installation (found at \"${TRILINOS_DIR}\")\n"
         "has to be configured with the same MPI configuration as deal.II, but found:\n"
         "  DEAL_II_WITH_MPI = ${DEAL_II_WITH_MPI}\n"
         "  TRILINOS_WITH_MPI = ${TRILINOS_WITH_MPI}\n"
@@ -152,8 +155,7 @@ MACRO(FEATURE_TRILINOS_FIND_EXTERNAL var)
             )
           SET(TRILINOS_ADDITIONAL_ERROR_STRING
             ${TRILINOS_ADDITIONAL_ERROR_STRING}
-            "The Trilinos installation found at\n"
-            "  ${TRILINOS_DIR}\n"
+            "The Trilinos installation (found at \"${TRILINOS_DIR}\")\n"
             "is not compatible with the C++ standard selected for\n"
             "this compiler. See the deal.II FAQ page for a solution.\n\n"
             )
index 5a6166acb6f95eb31b21d5e64fd0c3bd713bbcee..98617b69a8342436fa931355f331ec4622723b0a 100644 (file)
@@ -531,27 +531,6 @@ private:
    */
   bool at_newline;
 
-  /**
-   * The flags used to modify how
-   * the stream returned by get_stream() is used.
-   */
-  std::ios_base::fmtflags stream_flags;
-
-  /**
-   * The minimum width of characters
-   * used to represent a number. Used by width() to
-   * change the stream return by get_stream()
-   */
-  std::streamsize stream_width;
-
-  /**
-   * The maximum width of characters
-   * used to represent a floating point number.
-   * Used by precision() to
-   * change the stream return by get_stream()
-   */
-  std::streamsize stream_precision;
-
   /**
    * Print head of line. This prints optional time information and the
    * contents of the prefix stack.
@@ -563,43 +542,49 @@ private:
    * function will return the correct internal ostringstream buffer for
    * operater<<.
    */
-  std::ostringstream &get_stream()
-  {
-    return outstreams.get();
-  }
+  std::ostringstream &get_stream();
 
   /**
    * We use tbb's thread local storage facility to generate a stringstream
    * for every thread that sends log messages.
    */
-  Threads::ThreadLocalStorage<std::ostringstream> outstreams;
-
+  Threads::ThreadLocalStorage<std_cxx1x::shared_ptr<std::ostringstream> > outstreams;
 };
 
 
 /* ----------------------------- Inline functions and templates ---------------- */
 
 
+inline
+std::ostringstream &
+LogStream::get_stream()
+{
+  // see if we have already created this stream. if not, do so and
+  // set the default flags (why we set these flags is lost to
+  // history, but this is what we need to keep several hundred tests
+  // from producing different output)
+  //
+  // note that in all of this we need not worry about thread-safety
+  // because we operate on a thread-local object and by definition
+  // there can only be one access at a time
+  if (outstreams.get() == 0)
+    {
+      outstreams.get().reset (new std::ostringstream);
+      outstreams.get()->setf(std::ios::showpoint | std::ios::left);
+    }
+
+  // then return the stream
+  return *outstreams.get();
+}
+
+
 template <class T>
 inline
 LogStream &
 LogStream::operator<< (const T &t)
 {
-  // save the state of the output stream
-  std::ostringstream &stream = get_stream();
-  
-  const std::ios::fmtflags old_flags = stream.flags(stream_flags);
-  const unsigned int old_precision = stream.precision (stream_precision);
-  const unsigned int old_width = stream.width (stream_width);
-
-  // print to the internal stringstream, using the flags we have just set for
-  // the stream object:
-  stream << t;
-
-  // reset output format to whatever it was before
-  stream.flags (old_flags);
-  stream.precision(old_precision);
-  stream.width(old_width);
+  // print to the internal stringstream
+  get_stream() << t;
 
   return *this;
 }
@@ -610,11 +595,7 @@ inline
 LogStream &
 LogStream::operator<< (const double t)
 {
-  // save the state of out stream
   std::ostringstream &stream = get_stream();
-  const std::ios::fmtflags old_flags = stream.flags(stream_flags);
-  const unsigned int old_precision = stream.precision (stream_precision);
-  const unsigned int old_width = stream.width (stream_width);
 
   // we have to make sure that we don't catch NaN's and +-Inf's with the
   // test, because for these denormals all comparisons are always false.
@@ -627,11 +608,6 @@ LogStream::operator<< (const double t)
   else
     stream << t*(1.+offset);
 
-  // reset output format
-  stream.flags (old_flags);
-  stream.precision(old_precision);
-  stream.width(old_width);
-
   return *this;
 }
 
@@ -641,11 +617,7 @@ inline
 LogStream &
 LogStream::operator<< (const float t)
 {
-  // save the state of out stream
   std::ostringstream &stream = get_stream();
-  const std::ios::fmtflags old_flags = stream.flags(stream_flags);
-  const unsigned int old_precision = stream.precision (stream_precision);
-  const unsigned int old_width = stream.width (stream_width);
 
   // we have to make sure that we don't catch NaN's and +-Inf's with the
   // test, because for these denormals all comparisons are always false.
@@ -658,11 +630,6 @@ LogStream::operator<< (const float t)
   else
     stream << t*(1.+offset);
 
-  // reset output format
-  stream.flags (old_flags);
-  stream.precision(old_precision);
-  stream.width(old_width);
-
   return *this;
 }
 
index 067488d482ec36b2a2b8dbc60a559c9a08093022..9acb29a336fab92456f3ea4384bf6cd12b25ec1e 100644 (file)
@@ -3789,8 +3789,9 @@ FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
 namespace internal
 {
   // evaluates the given shape data in 1d-3d using the tensor product
-  // form. does not use the tensor product form and corresponds to a usual
-  // matrix-matrix product
+  // form. does not use a particular layout of entries in the matrices
+  // like the functions below and corresponds to a usual matrix-matrix
+  // product
   template <int dim, int fe_degree, int n_q_points_1d, typename Number,
            int direction, bool dof_to_quad, bool add>
   inline
@@ -3861,6 +3862,724 @@ namespace internal
 
 
 
+  // This method specializes the general application of tensor-product based
+  // elements for "symmetric" finite elements, i.e., when the shape functions
+  // are symmetric about 0.5 and the quadrature points are, too. In that case,
+  // the 1D shape values read (sorted lexicographically, rows run over 1D
+  // dofs, columns over quadrature points):
+  // Q2 --> [ 0.687  0 -0.087 ]
+  //        [ 0.4    1  0.4   ]
+  //        [-0.087  0  0.687 ]
+  // Q3 --> [ 0.66   0.003  0.002  0.049 ]
+  //        [ 0.521  1.005 -0.01  -0.230 ]
+  //        [-0.230 -0.01   1.005  0.521 ]
+  //        [ 0.049  0.002  0.003  0.66  ]
+  // Q4 --> [ 0.658  0.022  0 -0.007 -0.032 ]
+  //        [ 0.608  1.059  0  0.039  0.176 ]
+  //        [-0.409 -0.113  1 -0.113 -0.409 ]
+  //        [ 0.176  0.039  0  1.059  0.608 ]
+  //        [-0.032 -0.007  0  0.022  0.658 ]
+  //
+  // In these matrices, we want to use avoid computations involving zeros and
+  // ones and in addition use the symmetry in entries to reduce the number of
+  // read operations.
+  template <int dim, int fe_degree, int n_q_points_1d, typename Number,
+           int direction, bool dof_to_quad, bool add>
+  inline
+  void
+  apply_tensor_product_values (const VectorizedArray<Number> *shape_values,
+                              const VectorizedArray<Number> in [],
+                              VectorizedArray<Number>       out [])
+  {
+    AssertIndexRange (direction, dim);
+    const int mm     = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
+      nn     = dof_to_quad ? n_q_points_1d : (fe_degree+1);
+    const int n_cols = nn / 2;
+    const int mid    = mm / 2;
+
+    const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
+    const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
+    const int stride    = Utilities::fixed_int_power<nn,direction>::value;
+
+    for (int i2=0; i2<n_blocks2; ++i2)
+      {
+        for (int i1=0; i1<n_blocks1; ++i1)
+          {
+            for (int col=0; col<n_cols; ++col)
+              {
+               VectorizedArray<Number> val0, val1, in0, in1, res0, res1;
+               if (dof_to_quad == true)
+                 {
+                   val0 = shape_values[col];
+                   val1 = shape_values[nn-1-col];
+                 }
+               else
+                 {
+                   val0 = shape_values[col*n_q_points_1d];
+                   val1 = shape_values[(col+1)*n_q_points_1d-1];
+                 }
+               if (mid > 0)
+                 {
+                   in0 = in[0];
+                   in1 = in[stride*(mm-1)];
+                   res0 = val0 * in0;
+                   res1 = val1 * in0;
+                   res0 += val1 * in1;
+                   res1 += val0 * in1;
+                   for (int ind=1; ind<mid; ++ind)
+                     {
+                       if (dof_to_quad == true)
+                         {
+                           val0 = shape_values[ind*n_q_points_1d+col];
+                           val1 = shape_values[ind*n_q_points_1d+nn-1-col];
+                         }
+                       else
+                         {
+                           val0 = shape_values[col*n_q_points_1d+ind];
+                           val1 = shape_values[(col+1)*n_q_points_1d-1-ind];
+                         }
+                       in0 = in[stride*ind];
+                       in1 = in[stride*(mm-1-ind)];
+                       res0 += val0 * in0;
+                       res1 += val1 * in0;
+                       res0 += val1 * in1;
+                       res1 += val0 * in1;
+                     }
+                 }
+               else
+                 res0 = res1 = VectorizedArray<Number>();
+               if (dof_to_quad == true)
+                 {
+                   if (mm % 2 == 1)
+                     {
+                       val0 = shape_values[mid*n_q_points_1d+col];
+                       val1 = val0 * in[stride*mid];
+                       res0 += val1;
+                       res1 += val1;
+                     }
+                 }
+               else
+                 {
+                   if (mm % 2 == 1 && nn % 2 == 0)
+                     {
+                       val0 = shape_values[col*n_q_points_1d+mid];
+                       val1 = val0 * in[stride*mid];
+                       res0 += val1;
+                       res1 += val1;
+                     }
+                 }
+               if (add == false)
+                 {
+                   out[stride*col]         = res0;
+                   out[stride*(nn-1-col)]  = res1;
+                 }
+               else
+                 {
+                   out[stride*col]        += res0;
+                   out[stride*(nn-1-col)] += res1;
+                 }
+              }
+            if ( dof_to_quad == true && nn%2==1 && mm%2==1 )
+              {
+               if (add==false)
+                 out[stride*n_cols]  = in[stride*mid];
+               else
+                 out[stride*n_cols] += in[stride*mid];
+              }
+            else if (dof_to_quad == true && nn%2==1)
+              {
+               VectorizedArray<Number> res0;
+               VectorizedArray<Number> val0  = shape_values[n_cols];
+               if (mid > 0)
+                 {
+                   res0  = in[0] + in[stride*(mm-1)];
+                   res0 *= val0;
+                   for (int ind=1; ind<mid; ++ind)
+                     {
+                       val0  = shape_values[ind*n_q_points_1d+n_cols];
+                       VectorizedArray<Number> val1  = in[stride*ind] + in[stride*(mm-1-ind)];
+                       val1 *= val0;
+                       res0 += val1;
+                     }
+                 }
+               else
+                 res0 = VectorizedArray<Number>();
+               if (add == false)
+                 out[stride*n_cols]  = res0;
+               else
+                 out[stride*n_cols] += res0;
+              }
+            else if (dof_to_quad == false && nn%2 == 1)
+              {
+               VectorizedArray<Number> res0;
+               if (mid > 0)
+                 {
+                   VectorizedArray<Number> val0 = shape_values[n_cols*n_q_points_1d];
+                   res0 = in[0] + in[stride*(mm-1)];
+                   res0 *= val0;
+                   for (int ind=1; ind<mid; ++ind)
+                     {
+                       val0  = shape_values[n_cols*n_q_points_1d+ind];
+                       VectorizedArray<Number> val1 = in[stride*ind] + in[stride*(mm-1-ind)];
+                       val1 *= val0;
+                       res0 += val1;
+                     }
+                   if (mm % 2)
+                     res0 += in[stride*mid];
+                 }
+               else
+                 res0 = in[0];
+               if (add == false)
+                 out[stride*n_cols]  = res0;
+               else
+                 out[stride*n_cols] += res0;
+              }
+
+            // increment: in regular case, just go to the next point in
+            // x-direction. If we are at the end of one chunk in x-dir, need to
+            // jump over to the next layer in z-direction
+            switch (direction)
+              {
+              case 0:
+               in += mm;
+               out += nn;
+               break;
+              case 1:
+              case 2:
+               ++in;
+               ++out;
+               break;
+              default:
+               Assert (false, ExcNotImplemented());
+              }
+          }
+        if (direction == 1)
+          {
+            in += nn*(mm-1);
+            out += nn*(nn-1);
+          }
+      }
+  }
+
+
+
+  // evaluates the given shape data in 1d-3d using the tensor product
+  // form assuming the symmetries of unit cell shape gradients for
+  // finite elements in FEEvaluation
+
+  // For the specialized loop used for the gradient computation in
+  // here, the 1D shape values read (sorted lexicographically, rows
+  // run over 1D dofs, columns over quadrature points):
+  // Q2 --> [-2.549 -1  0.549 ]
+  //        [ 3.098  0 -3.098 ]
+  //        [-0.549  1  2.549 ]
+  // Q3 --> [-4.315 -1.03  0.5  -0.44  ]
+  //        [ 6.07  -1.44 -2.97  2.196 ]
+  //        [-2.196  2.97  1.44 -6.07  ]
+  //        [ 0.44  -0.5   1.03  4.315 ]
+  // Q4 --> [-6.316 -1.3    0.333 -0.353  0.413 ]
+  //        [10.111 -2.76  -2.667  2.066 -2.306 ]
+  //        [-5.688  5.773  0     -5.773  5.688 ]
+  //        [ 2.306 -2.066  2.667  2.76 -10.111 ]
+  //        [-0.413  0.353 -0.333 -0.353  0.413 ]
+  //
+  // In these matrices, we want to use avoid computations involving
+  // zeros and ones and in addition use the symmetry in entries to
+  // reduce the number of read operations.
+  template <int dim, int fe_degree, int n_q_points_1d, typename Number,
+           int direction, bool dof_to_quad, bool add>
+  inline
+  void
+  apply_tensor_product_gradients (const VectorizedArray<Number> *shape_gradients,
+                                 const VectorizedArray<Number> in [],
+                                 VectorizedArray<Number>       out [])
+  {
+    AssertIndexRange (direction, dim);
+    const int mm     = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
+      nn     = dof_to_quad ? n_q_points_1d : (fe_degree+1);
+    const int n_cols = nn / 2;
+    const int mid    = mm / 2;
+
+    const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
+    const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
+    const int stride    = Utilities::fixed_int_power<nn,direction>::value;
+
+    for (int i2=0; i2<n_blocks2; ++i2)
+      {
+        for (int i1=0; i1<n_blocks1; ++i1)
+          {
+            for (int col=0; col<n_cols; ++col)
+              {
+                VectorizedArray<Number> val0, val1, in0, in1, res0, res1;
+                if (dof_to_quad == true)
+                  {
+                    val0 = shape_gradients[col];
+                    val1 = shape_gradients[nn-1-col];
+                  }
+                else
+                  {
+                    val0 = shape_gradients[col*n_q_points_1d];
+                    val1 = shape_gradients[(nn-col-1)*n_q_points_1d];
+                  }
+                if (mid > 0)
+                  {
+                    in0 = in[0];
+                    in1 = in[stride*(mm-1)];
+                    res0 = val0 * in0;
+                    res1 = val1 * in0;
+                    res0 -= val1 * in1;
+                    res1 -= val0 * in1;
+                    for (int ind=1; ind<mid; ++ind)
+                      {
+                        if (dof_to_quad == true)
+                          {
+                            val0 = shape_gradients[ind*n_q_points_1d+col];
+                            val1 = shape_gradients[ind*n_q_points_1d+nn-1-col];
+                          }
+                        else
+                          {
+                            val0 = shape_gradients[col*n_q_points_1d+ind];
+                            val1 = shape_gradients[(nn-col-1)*n_q_points_1d+ind];
+                          }
+                        in0 = in[stride*ind];
+                        in1 = in[stride*(mm-1-ind)];
+                        res0 += val0 * in0;
+                        res1 += val1 * in0;
+                        res0 -= val1 * in1;
+                        res1 -= val0 * in1;
+                      }
+                  }
+                else
+                  res0 = res1 = VectorizedArray<Number>();
+                if (mm % 2 == 1)
+                  {
+                    if (dof_to_quad == true)
+                      val0 = shape_gradients[mid*n_q_points_1d+col];
+                    else
+                      val0 = shape_gradients[col*n_q_points_1d+mid];
+                    val1 = val0 * in[stride*mid];
+                    res0 += val1;
+                    res1 -= val1;
+                  }
+                if (add == false)
+                  {
+                    out[stride*col]         = res0;
+                    out[stride*(nn-1-col)]  = res1;
+                  }
+                else
+                  {
+                    out[stride*col]        += res0;
+                    out[stride*(nn-1-col)] += res1;
+                  }
+              }
+            if ( nn%2 == 1 )
+              {
+                VectorizedArray<Number> val0, res0;
+                if (dof_to_quad == true)
+                  val0 = shape_gradients[n_cols];
+                else
+                  val0 = shape_gradients[n_cols*n_q_points_1d];
+                res0  = in[0] - in[stride*(mm-1)];
+                res0 *= val0;
+                for (int ind=1; ind<mid; ++ind)
+                  {
+                    if (dof_to_quad == true)
+                      val0 = shape_gradients[ind*n_q_points_1d+n_cols];
+                    else
+                      val0 = shape_gradients[n_cols*n_q_points_1d+ind];
+                    VectorizedArray<Number> val1  = in[stride*ind] - in[stride*(mm-1-ind)];
+                    val1 *= val0;
+                    res0 += val1;
+                  }
+                if (add == false)
+                  out[stride*n_cols]  = res0;
+                else
+                  out[stride*n_cols] += res0;
+              }
+
+            // increment: in regular case, just go to the next point in
+            // x-direction. for y-part in 3D and if we are at the end of one
+            // chunk in x-dir, need to jump over to the next layer in
+            // z-direction
+            switch (direction)
+              {
+              case 0:
+                in += mm;
+                out += nn;
+                break;
+              case 1:
+              case 2:
+                ++in;
+                ++out;
+                break;
+              default:
+                Assert (false, ExcNotImplemented());
+              }
+          }
+
+        if (direction == 1)
+          {
+            in  += nn * (mm-1);
+            out += nn * (nn-1);
+          }
+      }
+  }
+
+
+
+  // evaluates the given shape data in 1d-3d using the tensor product
+  // form assuming the symmetries of unit cell shape hessians for
+  // finite elements in FEEvaluation
+  template <int dim, int fe_degree, int n_q_points_1d, typename Number,
+           int direction, bool dof_to_quad, bool add>
+  inline
+  void
+  apply_tensor_product_hessians (const VectorizedArray<Number> *shape_hessians,
+                                const VectorizedArray<Number> in [],
+                                VectorizedArray<Number>       out [])
+  {
+    AssertIndexRange (direction, dim);
+    const int mm     = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
+      nn     = dof_to_quad ? n_q_points_1d : (fe_degree+1);
+    const int n_cols = nn / 2;
+    const int mid    = mm / 2;
+
+    const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
+    const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
+    const int stride    = Utilities::fixed_int_power<nn,direction>::value;
+
+    for (int i2=0; i2<n_blocks2; ++i2)
+      {
+        for (int i1=0; i1<n_blocks1; ++i1)
+          {
+            for (int col=0; col<n_cols; ++col)
+              {
+                VectorizedArray<Number> val0, val1, in0, in1, res0, res1;
+                if (dof_to_quad == true)
+                  {
+                    val0 = shape_hessians[col];
+                    val1 = shape_hessians[nn-1-col];
+                  }
+                else
+                  {
+                    val0 = shape_hessians[col*n_q_points_1d];
+                    val1 = shape_hessians[(col+1)*n_q_points_1d-1];
+                  }
+                if (mid > 0)
+                  {
+                    in0 = in[0];
+                    in1 = in[stride*(mm-1)];
+                    res0 = val0 * in0;
+                    res1 = val1 * in0;
+                    res0 += val1 * in1;
+                    res1 += val0 * in1;
+                    for (int ind=1; ind<mid; ++ind)
+                      {
+                        if (dof_to_quad == true)
+                          {
+                            val0 = shape_hessians[ind*n_q_points_1d+col];
+                            val1 = shape_hessians[ind*n_q_points_1d+nn-1-col];
+                          }
+                        else
+                          {
+                            val0 = shape_hessians[col*n_q_points_1d+ind];
+                            val1 = shape_hessians[(col+1)*n_q_points_1d-1-ind];
+                          }
+                        in0 = in[stride*ind];
+                        in1 = in[stride*(mm-1-ind)];
+                        res0 += val0 * in0;
+                        res1 += val1 * in0;
+                        res0 += val1 * in1;
+                        res1 += val0 * in1;
+                      }
+                  }
+                else
+                  res0 = res1 = VectorizedArray<Number>();
+                if (mm % 2 == 1)
+                  {
+                    if (dof_to_quad == true)
+                      val0 = shape_hessians[mid*n_q_points_1d+col];
+                    else
+                      val0 = shape_hessians[col*n_q_points_1d+mid];
+                    val1 = val0 * in[stride*mid];
+                    res0 += val1;
+                    res1 += val1;
+                  }
+                if (add == false)
+                  {
+                    out[stride*col]         = res0;
+                    out[stride*(nn-1-col)]  = res1;
+                  }
+                else
+                  {
+                    out[stride*col]        += res0;
+                    out[stride*(nn-1-col)] += res1;
+                  }
+              }
+            if ( nn%2 == 1 )
+              {
+                VectorizedArray<Number> val0, res0;
+                if (dof_to_quad == true)
+                  val0 = shape_hessians[n_cols];
+                else
+                  val0 = shape_hessians[n_cols*n_q_points_1d];
+                if (mid > 0)
+                  {
+                    res0  = in[0] + in[stride*(mm-1)];
+                    res0 *= val0;
+                    for (int ind=1; ind<mid; ++ind)
+                      {
+                        if (dof_to_quad == true)
+                          val0 = shape_hessians[ind*n_q_points_1d+n_cols];
+                        else
+                          val0 = shape_hessians[n_cols*n_q_points_1d+ind];
+                        VectorizedArray<Number> val1  = in[stride*ind] + in[stride*(mm-1-ind)];
+                        val1 *= val0;
+                        res0 += val1;
+                      }
+                  }
+                else
+                  res0 = VectorizedArray<Number>();
+                if (mm % 2 == 1)
+                  {
+                    if (dof_to_quad == true)
+                      val0 = shape_hessians[mid*n_q_points_1d+n_cols];
+                    else
+                      val0 = shape_hessians[n_cols*n_q_points_1d+mid];
+                    res0 += val0 * in[stride*mid];
+                  }
+                if (add == false)
+                  out[stride*n_cols]  = res0;
+                else
+                  out[stride*n_cols] += res0;
+              }
+
+            // increment: in regular case, just go to the next point in
+            // x-direction. If we are at the end of one chunk in x-dir, need to
+            // jump over to the next layer in z-direction
+            switch (direction)
+              {
+              case 0:
+                in += mm;
+                out += nn;
+                break;
+              case 1:
+              case 2:
+                ++in;
+                ++out;
+                break;
+              default:
+                Assert (false, ExcNotImplemented());
+              }
+          }
+        if (direction == 1)
+          {
+            in += nn*(mm-1);
+            out += nn*(nn-1);
+          }
+      }
+  }
+
+
+
+  // evaluates the given shape data in 1d-3d using the tensor product
+  // form assuming the symmetries of unit cell shape gradients for
+  // finite elements in FEEvaluationGL
+
+  // This function specializes the application of the tensor product loop for
+  // Gauss-Lobatto elements which are symmetric about 0.5 just as the general
+  // class of elements treated by FEEvaluation, have diagonal shape matrices
+  // for the values and have the following gradient matrices (notice the zeros
+  // on the diagonal in the interior points, which is due to the construction
+  // of Legendre polynomials):
+  // Q2 --> [-3 -1  1 ]
+  //        [ 4  0 -4 ]
+  //        [-1  1  3 ]
+  // Q3 --> [-6    -1.618  0.618 -1    ]
+  //        [ 8.09  0     -2.236  3.09 ]
+  //        [-3.09  2.236  0     -8.09 ]
+  //        [ 1    -0.618  1.618  6    ]
+  // Q4 --> [-10    -2.482  0.75  -0.518  1     ]
+  //        [ 13.51  0     -2.673  1.528 -2.82  ]
+  //        [-5.333  3.491  0     -3.491  5.333 ]
+  //        [ 2.82  -1.528  2.673  0    -13.51  ]
+  //        [-1      0.518 -0.75   2.482 10     ]
+  template <int dim, int fe_degree, typename Number,
+           int direction, bool dof_to_quad, bool add>
+  inline
+  void
+  apply_tensor_product_gradients_gl (const VectorizedArray<Number> *shape_gradients,
+                                    const VectorizedArray<Number> in [],
+                                    VectorizedArray<Number>       out [])
+  {
+    AssertIndexRange (direction, dim);
+    const int mm     = fe_degree+1;
+    const int nn     = fe_degree+1;
+    const int n_cols = nn / 2;
+    const int mid    = mm / 2;
+
+    const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
+    const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
+    const int stride    = Utilities::fixed_int_power<nn,direction>::value;
+
+    for (int i2=0; i2<n_blocks2; ++i2)
+      {
+        for (int i1=0; i1<n_blocks1; ++i1)
+          {
+            for (int col=0; col<n_cols; ++col)
+              {
+                VectorizedArray<Number> val0, val1, in0, in1, res0, res1;
+                if (mid > 0)
+                  {
+                    if (dof_to_quad == true)
+                      {
+                        val0 = shape_gradients[col];
+                        val1 = shape_gradients[nn-1-col];
+                      }
+                    else
+                      {
+                        val0 = shape_gradients[col*mm];
+                        val1 = shape_gradients[(nn-col-1)*mm];
+                      }
+                    in0 = in[0];
+                    in1 = in[stride*(mm-1)];
+                    if (col == 0)
+                      {
+                        if ((mm+dof_to_quad)%2 == 1)
+                          {
+                            res0 = val0 * in0;
+                            res1 = -in0;
+                            res0 += in1;
+                            res1 -= val0 * in1;
+                          }
+                        else
+                          {
+                            res0 = val0 * in0;
+                            res0 -= in1;
+                            res1 = in0;
+                            res1 -= val0 * in1;
+                          }
+                      }
+                    else
+                      {
+                        res0 = val0 * in0;
+                        res1 = val1 * in0;
+                        res0 -= val1 * in1;
+                        res1 -= val0 * in1;
+                      }
+                    for (int ind=1; ind<mid; ++ind)
+                      {
+                        if (dof_to_quad == true)
+                          {
+                            val0 = shape_gradients[ind*mm+col];
+                            val1 = shape_gradients[ind*mm+nn-1-col];
+                          }
+                        else
+                          {
+                            val0 = shape_gradients[col*mm+ind];
+                            val1 = shape_gradients[(nn-col-1)*mm+ind];
+                          }
+
+                        // at inner points, the gradient is zero for ind==col
+                        in0 = in[stride*ind];
+                        in1 = in[stride*(mm-1-ind)];
+                        if (ind == col)
+                          {
+                            res1 += val1 * in0;
+                            res0 -= val1 * in1;
+                          }
+                        else
+                          {
+                            res0 += val0 * in0;
+                            res1 += val1 * in0;
+                            res0 -= val1 * in1;
+                            res1 -= val0 * in1;
+                          }
+                      }
+                  }
+                else
+                  res0 = res1 = VectorizedArray<Number>();
+                if (mm % 2 == 1)
+                  {
+                    if (dof_to_quad == true)
+                      val0 = shape_gradients[mid*mm+col];
+                    else
+                      val0 = shape_gradients[col*mm+mid];
+                    val1 = val0 * in[stride*mid];
+                    res0 += val1;
+                    res1 -= val1;
+                  }
+                if (add == false)
+                  {
+                    out[stride*col]         = res0;
+                    out[stride*(nn-1-col)]  = res1;
+                  }
+                else
+                  {
+                    out[stride*col]        += res0;
+                    out[stride*(nn-1-col)] += res1;
+                  }
+              }
+            if ( nn%2 == 1 )
+              {
+                VectorizedArray<Number> val0, res0;
+                if (dof_to_quad == true)
+                  val0 = shape_gradients[n_cols];
+                else
+                  val0 = shape_gradients[n_cols*mm];
+                if (mid > 0)
+                  {
+                    res0  = in[0] - in[stride*(mm-1)];
+                    res0 *= val0;
+                    for (int ind=1; ind<mid; ++ind)
+                      {
+                        if (dof_to_quad == true)
+                          val0 = shape_gradients[ind*mm+n_cols];
+                        else
+                          val0 = shape_gradients[n_cols*mm+ind];
+                        VectorizedArray<Number> val1  = in[stride*ind] - in[stride*(mm-1-ind)];
+                        val1 *= val0;
+                        res0 += val1;
+                      }
+                  }
+                else
+                  res0 = VectorizedArray<Number>();
+                if (add == false)
+                  out[stride*n_cols]  = res0;
+                else
+                  out[stride*n_cols] += res0;
+              }
+
+            // increment: in regular case, just go to the next point in
+            // x-direction. for y-part in 3D and if we are at the end of one
+            // chunk in x-dir, need to jump over to the next layer in
+            // z-direction
+            switch (direction)
+              {
+              case 0:
+                in += mm;
+                out += nn;
+                break;
+              case 1:
+              case 2:
+                ++in;
+                ++out;
+                break;
+              default:
+                Assert (false, ExcNotImplemented());
+              }
+          }
+
+        if (direction == 1)
+          {
+            in  += nn * (mm-1);
+            out += nn * (nn-1);
+          }
+      }
+  }
+
+
+
   // This performs the evaluation of function values, gradients and Hessians
   // for tensor-product finite elements. The operation is used for both
   // FEEvaluationGeneral and FEEvaluation, which provide different functions
@@ -4393,8 +5112,6 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
 
 
 
-/*----------------- optimized implementation tensor product symmetric case --*/
-
 template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
          typename Number>
 template <int direction, bool dof_to_quad, bool add>
@@ -4404,196 +5121,9 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
 ::apply_values (const VectorizedArray<Number> in [],
                 VectorizedArray<Number>       out [])
 {
-  AssertIndexRange (direction, dim);
-  const int mm     = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
-            nn     = dof_to_quad ? n_q_points_1d : (fe_degree+1);
-  const int n_cols = nn / 2;
-  const int mid    = mm / 2;
-
-  const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
-  const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
-  const int stride    = Utilities::fixed_int_power<nn,direction>::value;
-
-  // This loop specializes the general application of tensor-product based
-  // elements for "symmetric" finite elements, i.e., when the shape functions
-  // are symmetric about 0.5 and the quadrature points are, too. In that case,
-  // the 1D shape values read (sorted lexicographically, rows run over 1D
-  // dofs, columns over quadrature points):
-  // Q2 --> [ 0.687  0 -0.087 ]
-  //        [ 0.4    1  0.4   ]
-  //        [-0.087  0  0.687 ]
-  // Q3 --> [ 0.66   0.003  0.002  0.049 ]
-  //        [ 0.521  1.005 -0.01  -0.230 ]
-  //        [-0.230 -0.01   1.005  0.521 ]
-  //        [ 0.049  0.002  0.003  0.66  ]
-  // Q4 --> [ 0.658  0.022  0 -0.007 -0.032 ]
-  //        [ 0.608  1.059  0  0.039  0.176 ]
-  //        [-0.409 -0.113  1 -0.113 -0.409 ]
-  //        [ 0.176  0.039  0  1.059  0.608 ]
-  //        [-0.032 -0.007  0  0.022  0.658 ]
-  //
-  // In these matrices, we want to use avoid computations involving zeros and
-  // ones and in addition use the symmetry in entries to reduce the number of
-  // read operations.
-  const VectorizedArray<Number> *shape_values = this->data.shape_values.begin();
-  for (int i2=0; i2<n_blocks2; ++i2)
-    {
-      for (int i1=0; i1<n_blocks1; ++i1)
-        {
-          for (int col=0; col<n_cols; ++col)
-            {
-              VectorizedArray<Number> val0, val1, in0, in1, res0, res1;
-              if (dof_to_quad == true)
-                {
-                  val0 = shape_values[col];
-                  val1 = shape_values[nn-1-col];
-                }
-              else
-                {
-                  val0 = shape_values[col*n_q_points_1d];
-                  val1 = shape_values[(col+1)*n_q_points_1d-1];
-                }
-              if (mid > 0)
-                {
-                  in0 = in[0];
-                  in1 = in[stride*(mm-1)];
-                  res0 = val0 * in0;
-                  res1 = val1 * in0;
-                  res0 += val1 * in1;
-                  res1 += val0 * in1;
-                  for (int ind=1; ind<mid; ++ind)
-                    {
-                      if (dof_to_quad == true)
-                        {
-                          val0 = shape_values[ind*n_q_points_1d+col];
-                          val1 = shape_values[ind*n_q_points_1d+nn-1-col];
-                        }
-                      else
-                        {
-                          val0 = shape_values[col*n_q_points_1d+ind];
-                          val1 = shape_values[(col+1)*n_q_points_1d-1-ind];
-                        }
-                      in0 = in[stride*ind];
-                      in1 = in[stride*(mm-1-ind)];
-                      res0 += val0 * in0;
-                      res1 += val1 * in0;
-                      res0 += val1 * in1;
-                      res1 += val0 * in1;
-                    }
-                }
-              else
-                res0 = res1 = VectorizedArray<Number>();
-              if (dof_to_quad == true)
-                {
-                  if (mm % 2 == 1)
-                    {
-                      val0 = shape_values[mid*n_q_points_1d+col];
-                      val1 = val0 * in[stride*mid];
-                      res0 += val1;
-                      res1 += val1;
-                    }
-                }
-              else
-                {
-                  if (mm % 2 == 1 && nn % 2 == 0)
-                    {
-                      val0 = shape_values[col*n_q_points_1d+mid];
-                      val1 = val0 * in[stride*mid];
-                      res0 += val1;
-                      res1 += val1;
-                    }
-                }
-              if (add == false)
-                {
-                  out[stride*col]         = res0;
-                  out[stride*(nn-1-col)]  = res1;
-                }
-              else
-                {
-                  out[stride*col]        += res0;
-                  out[stride*(nn-1-col)] += res1;
-                }
-            }
-          if ( dof_to_quad == true && nn%2==1 && mm%2==1 )
-            {
-              if (add==false)
-                out[stride*n_cols]  = in[stride*mid];
-              else
-                out[stride*n_cols] += in[stride*mid];
-            }
-          else if (dof_to_quad == true && nn%2==1)
-            {
-              VectorizedArray<Number> res0;
-              VectorizedArray<Number> val0  = shape_values[n_cols];
-              if (mid > 0)
-                {
-                  res0  = in[0] + in[stride*(mm-1)];
-                  res0 *= val0;
-                  for (int ind=1; ind<mid; ++ind)
-                    {
-                      val0  = shape_values[ind*n_q_points_1d+n_cols];
-                      VectorizedArray<Number> val1  = in[stride*ind] + in[stride*(mm-1-ind)];
-                      val1 *= val0;
-                      res0 += val1;
-                    }
-                }
-              else
-                res0 = VectorizedArray<Number>();
-              if (add == false)
-                out[stride*n_cols]  = res0;
-              else
-                out[stride*n_cols] += res0;
-            }
-          else if (dof_to_quad == false && nn%2 == 1)
-            {
-              VectorizedArray<Number> res0;
-              if (mid > 0)
-                {
-                  VectorizedArray<Number> val0 = shape_values[n_cols*n_q_points_1d];
-                  res0 = in[0] + in[stride*(mm-1)];
-                  res0 *= val0;
-                  for (int ind=1; ind<mid; ++ind)
-                    {
-                      val0  = shape_values[n_cols*n_q_points_1d+ind];
-                      VectorizedArray<Number> val1 = in[stride*ind] + in[stride*(mm-1-ind)];
-                      val1 *= val0;
-                      res0 += val1;
-                    }
-                  if (mm % 2)
-                    res0 += in[stride*mid];
-                }
-              else
-                res0 = in[0];
-              if (add == false)
-                out[stride*n_cols]  = res0;
-              else
-                out[stride*n_cols] += res0;
-            }
-
-          // increment: in regular case, just go to the next point in
-          // x-direction. If we are at the end of one chunk in x-dir, need to
-          // jump over to the next layer in z-direction
-          switch (direction)
-            {
-            case 0:
-              in += mm;
-              out += nn;
-              break;
-            case 1:
-            case 2:
-              ++in;
-              ++out;
-              break;
-            default:
-              Assert (false, ExcNotImplemented());
-            }
-        }
-      if (direction == 1)
-        {
-          in += nn*(mm-1);
-          out += nn*(nn-1);
-        }
-    }
+  internal::apply_tensor_product_values<dim,fe_degree,n_q_points_1d,Number,
+           direction, dof_to_quad, add>
+           (this->data.shape_values.begin(), in, out);
 }
 
 
@@ -4607,155 +5137,9 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
 ::apply_gradients (const VectorizedArray<Number> in [],
                    VectorizedArray<Number>       out [])
 {
-  AssertIndexRange (direction, dim);
-  const int mm     = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
-            nn     = dof_to_quad ? n_q_points_1d : (fe_degree+1);
-  const int n_cols = nn / 2;
-  const int mid    = mm / 2;
-
-  const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
-  const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
-  const int stride    = Utilities::fixed_int_power<nn,direction>::value;
-
-  const VectorizedArray<Number> *shape_gradients = this->data.shape_gradients.begin();
-  for (int i2=0; i2<n_blocks2; ++i2)
-    {
-      for (int i1=0; i1<n_blocks1; ++i1)
-        {
-          // For the specialized loop used for the gradient computation in
-          // here, the 1D shape values read (sorted lexicographically, rows
-          // run over 1D dofs, columns over quadrature points):
-          // Q2 --> [-2.549 -1  0.549 ]
-          //        [ 3.098  0 -3.098 ]
-          //        [-0.549  1  2.549 ]
-          // Q3 --> [-4.315 -1.03  0.5  -0.44  ]
-          //        [ 6.07  -1.44 -2.97  2.196 ]
-          //        [-2.196  2.97  1.44 -6.07  ]
-          //        [ 0.44  -0.5   1.03  4.315 ]
-          // Q4 --> [-6.316 -1.3    0.333 -0.353  0.413 ]
-          //        [10.111 -2.76  -2.667  2.066 -2.306 ]
-          //        [-5.688  5.773  0     -5.773  5.688 ]
-          //        [ 2.306 -2.066  2.667  2.76 -10.111 ]
-          //        [-0.413  0.353 -0.333 -0.353  0.413 ]
-          //
-          // In these matrices, we want to use avoid computations involving
-          // zeros and ones and in addition use the symmetry in entries to
-          // reduce the number of read operations.
-          for (int col=0; col<n_cols; ++col)
-            {
-              VectorizedArray<Number> val0, val1, in0, in1, res0, res1;
-              if (dof_to_quad == true)
-                {
-                  val0 = shape_gradients[col];
-                  val1 = shape_gradients[nn-1-col];
-                }
-              else
-                {
-                  val0 = shape_gradients[col*n_q_points_1d];
-                  val1 = shape_gradients[(nn-col-1)*n_q_points_1d];
-                }
-              if (mid > 0)
-                {
-                  in0 = in[0];
-                  in1 = in[stride*(mm-1)];
-                  res0 = val0 * in0;
-                  res1 = val1 * in0;
-                  res0 -= val1 * in1;
-                  res1 -= val0 * in1;
-                  for (int ind=1; ind<mid; ++ind)
-                    {
-                      if (dof_to_quad == true)
-                        {
-                          val0 = shape_gradients[ind*n_q_points_1d+col];
-                          val1 = shape_gradients[ind*n_q_points_1d+nn-1-col];
-                        }
-                      else
-                        {
-                          val0 = shape_gradients[col*n_q_points_1d+ind];
-                          val1 = shape_gradients[(nn-col-1)*n_q_points_1d+ind];
-                        }
-                      in0 = in[stride*ind];
-                      in1 = in[stride*(mm-1-ind)];
-                      res0 += val0 * in0;
-                      res1 += val1 * in0;
-                      res0 -= val1 * in1;
-                      res1 -= val0 * in1;
-                    }
-                }
-              else
-                res0 = res1 = VectorizedArray<Number>();
-              if (mm % 2 == 1)
-                {
-                  if (dof_to_quad == true)
-                    val0 = shape_gradients[mid*n_q_points_1d+col];
-                  else
-                    val0 = shape_gradients[col*n_q_points_1d+mid];
-                  val1 = val0 * in[stride*mid];
-                  res0 += val1;
-                  res1 -= val1;
-                }
-              if (add == false)
-                {
-                  out[stride*col]         = res0;
-                  out[stride*(nn-1-col)]  = res1;
-                }
-              else
-                {
-                  out[stride*col]        += res0;
-                  out[stride*(nn-1-col)] += res1;
-                }
-            }
-          if ( nn%2 == 1 )
-            {
-              VectorizedArray<Number> val0, res0;
-              if (dof_to_quad == true)
-                val0 = shape_gradients[n_cols];
-              else
-                val0 = shape_gradients[n_cols*n_q_points_1d];
-              res0  = in[0] - in[stride*(mm-1)];
-              res0 *= val0;
-              for (int ind=1; ind<mid; ++ind)
-                {
-                  if (dof_to_quad == true)
-                    val0 = shape_gradients[ind*n_q_points_1d+n_cols];
-                  else
-                    val0 = shape_gradients[n_cols*n_q_points_1d+ind];
-                  VectorizedArray<Number> val1  = in[stride*ind] - in[stride*(mm-1-ind)];
-                  val1 *= val0;
-                  res0 += val1;
-                }
-              if (add == false)
-                out[stride*n_cols]  = res0;
-              else
-                out[stride*n_cols] += res0;
-            }
-
-          // increment: in regular case, just go to the next point in
-          // x-direction. for y-part in 3D and if we are at the end of one
-          // chunk in x-dir, need to jump over to the next layer in
-          // z-direction
-          switch (direction)
-            {
-            case 0:
-              in += mm;
-              out += nn;
-              break;
-            case 1:
-            case 2:
-              ++in;
-              ++out;
-              break;
-            default:
-              Assert (false, ExcNotImplemented());
-            }
-        }
-
-      if (direction == 1)
-        {
-          in  += nn * (mm-1);
-          out += nn * (nn-1);
-        }
-    }
+  internal::apply_tensor_product_gradients<dim,fe_degree,n_q_points_1d,Number,
+           direction, dof_to_quad, add>
+           (this->data.shape_gradients.begin(), in, out);
 }
 
 
@@ -4772,146 +5156,9 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
 ::apply_hessians (const VectorizedArray<Number> in [],
                   VectorizedArray<Number>       out [])
 {
-  AssertIndexRange (direction, dim);
-  const int mm     = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
-            nn     = dof_to_quad ? n_q_points_1d : (fe_degree+1);
-  const int n_cols = nn / 2;
-  const int mid    = mm / 2;
-
-  const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
-  const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
-  const int stride    = Utilities::fixed_int_power<nn,direction>::value;
-
-  for (int i2=0; i2<n_blocks2; ++i2)
-    {
-      for (int i1=0; i1<n_blocks1; ++i1)
-        {
-          for (int col=0; col<n_cols; ++col)
-            {
-              VectorizedArray<Number> val0, val1, in0, in1, res0, res1;
-              if (dof_to_quad == true)
-                {
-                  val0 = this->data.shape_hessians[col];
-                  val1 = this->data.shape_hessians[nn-1-col];
-                }
-              else
-                {
-                  val0 = this->data.shape_hessians[col*n_q_points_1d];
-                  val1 = this->data.shape_hessians[(col+1)*n_q_points_1d-1];
-                }
-              if (mid > 0)
-                {
-                  in0 = in[0];
-                  in1 = in[stride*(mm-1)];
-                  res0 = val0 * in0;
-                  res1 = val1 * in0;
-                  res0 += val1 * in1;
-                  res1 += val0 * in1;
-                  for (int ind=1; ind<mid; ++ind)
-                    {
-                      if (dof_to_quad == true)
-                        {
-                          val0 = this->data.shape_hessians[ind*n_q_points_1d+col];
-                          val1 = this->data.shape_hessians[ind*n_q_points_1d+nn-1-col];
-                        }
-                      else
-                        {
-                          val0 = this->data.shape_hessians[col*n_q_points_1d+ind];
-                          val1 = this->data.shape_hessians[(col+1)*n_q_points_1d-1-ind];
-                        }
-                      in0 = in[stride*ind];
-                      in1 = in[stride*(mm-1-ind)];
-                      res0 += val0 * in0;
-                      res1 += val1 * in0;
-                      res0 += val1 * in1;
-                      res1 += val0 * in1;
-                    }
-                }
-              else
-                res0 = res1 = VectorizedArray<Number>();
-              if (mm % 2 == 1)
-                {
-                  if (dof_to_quad == true)
-                    val0 = this->data.shape_hessians[mid*n_q_points_1d+col];
-                  else
-                    val0 = this->data.shape_hessians[col*n_q_points_1d+mid];
-                  val1 = val0 * in[stride*mid];
-                  res0 += val1;
-                  res1 += val1;
-                }
-              if (add == false)
-                {
-                  out[stride*col]         = res0;
-                  out[stride*(nn-1-col)]  = res1;
-                }
-              else
-                {
-                  out[stride*col]        += res0;
-                  out[stride*(nn-1-col)] += res1;
-                }
-            }
-          if ( nn%2 == 1 )
-            {
-              VectorizedArray<Number> val0, res0;
-              if (dof_to_quad == true)
-                val0 = this->data.shape_hessians[n_cols];
-              else
-                val0 = this->data.shape_hessians[n_cols*n_q_points_1d];
-              if (mid > 0)
-                {
-                  res0  = in[0] + in[stride*(mm-1)];
-                  res0 *= val0;
-                  for (int ind=1; ind<mid; ++ind)
-                    {
-                      if (dof_to_quad == true)
-                        val0 = this->data.shape_hessians[ind*n_q_points_1d+n_cols];
-                      else
-                        val0 = this->data.shape_hessians[n_cols*n_q_points_1d+ind];
-                      VectorizedArray<Number> val1  = in[stride*ind] + in[stride*(mm-1-ind)];
-                      val1 *= val0;
-                      res0 += val1;
-                    }
-                }
-              else
-                res0 = VectorizedArray<Number>();
-              if (mm % 2 == 1)
-                {
-                  if (dof_to_quad == true)
-                    val0 = this->data.shape_hessians[mid*n_q_points_1d+n_cols];
-                  else
-                    val0 = this->data.shape_hessians[n_cols*n_q_points_1d+mid];
-                  res0 += val0 * in[stride*mid];
-                }
-              if (add == false)
-                out[stride*n_cols]  = res0;
-              else
-                out[stride*n_cols] += res0;
-            }
-
-          // increment: in regular case, just go to the next point in
-          // x-direction. If we are at the end of one chunk in x-dir, need to
-          // jump over to the next layer in z-direction
-          switch (direction)
-            {
-            case 0:
-              in += mm;
-              out += nn;
-              break;
-            case 1:
-            case 2:
-              ++in;
-              ++out;
-              break;
-            default:
-              Assert (false, ExcNotImplemented());
-            }
-        }
-      if (direction == 1)
-        {
-          in += nn*(mm-1);
-          out += nn*(nn-1);
-        }
-    }
+  internal::apply_tensor_product_hessians<dim,fe_degree,n_q_points_1d,Number,
+           direction, dof_to_quad, add>
+           (this->data.shape_hessians.begin(), in, out);
 }
 
 
@@ -5152,193 +5399,14 @@ FEEvaluationGL<dim,fe_degree,n_components_,Number>
 ::apply_gradients (const VectorizedArray<Number> in [],
                    VectorizedArray<Number>       out [])
 {
-  AssertIndexRange (direction, dim);
-  const int mm     = fe_degree+1;
-  const int nn     = fe_degree+1;
-  const int n_cols = nn / 2;
-  const int mid    = mm / 2;
-
-  const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
-  const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
-  const int stride    = Utilities::fixed_int_power<nn,direction>::value;
-
-  // This loop specializes the application of the tensor product loop for
-  // Gauss-Lobatto elements which are symmetric about 0.5 just as the general
-  // class of elements treated by FEEvaluation, have diagonal shape matrices
-  // for the values and have the following gradient matrices (notice the zeros
-  // on the diagonal in the interior points, which is due to the construction
-  // of Legendre polynomials):
-  // Q2 --> [-3 -1  1 ]
-  //        [ 4  0 -4 ]
-  //        [-1  1  3 ]
-  // Q3 --> [-6    -1.618  0.618 -1    ]
-  //        [ 8.09  0     -2.236  3.09 ]
-  //        [-3.09  2.236  0     -8.09 ]
-  //        [ 1    -0.618  1.618  6    ]
-  // Q4 --> [-10    -2.482  0.75  -0.518  1     ]
-  //        [ 13.51  0     -2.673  1.528 -2.82  ]
-  //        [-5.333  3.491  0     -3.491  5.333 ]
-  //        [ 2.82  -1.528  2.673  0    -13.51  ]
-  //        [-1      0.518 -0.75   2.482 10     ]
-  for (int i2=0; i2<n_blocks2; ++i2)
-    {
-      for (int i1=0; i1<n_blocks1; ++i1)
-        {
-          for (int col=0; col<n_cols; ++col)
-            {
-              VectorizedArray<Number> val0, val1, in0, in1, res0, res1;
-              if (mid > 0)
-                {
-                  if (dof_to_quad == true)
-                    {
-                      val0 = this->data.shape_gradients[col];
-                      val1 = this->data.shape_gradients[nn-1-col];
-                    }
-                  else
-                    {
-                      val0 = this->data.shape_gradients[col*mm];
-                      val1 = this->data.shape_gradients[(nn-col-1)*mm];
-                    }
-                  in0 = in[0];
-                  in1 = in[stride*(mm-1)];
-                  if (col == 0)
-                    {
-                      if ((mm+dof_to_quad)%2 == 1)
-                        {
-                          res0 = val0 * in0;
-                          res1 = -in0;
-                          res0 += in1;
-                          res1 -= val0 * in1;
-                        }
-                      else
-                        {
-                          res0 = val0 * in0;
-                          res0 -= in1;
-                          res1 = in0;
-                          res1 -= val0 * in1;
-                        }
-                    }
-                  else
-                    {
-                      res0 = val0 * in0;
-                      res1 = val1 * in0;
-                      res0 -= val1 * in1;
-                      res1 -= val0 * in1;
-                    }
-                  for (int ind=1; ind<mid; ++ind)
-                    {
-                      if (dof_to_quad == true)
-                        {
-                          val0 = this->data.shape_gradients[ind*mm+col];
-                          val1 = this->data.shape_gradients[ind*mm+nn-1-col];
-                        }
-                      else
-                        {
-                          val0 = this->data.shape_gradients[col*mm+ind];
-                          val1 = this->data.shape_gradients[(nn-col-1)*mm+ind];
-                        }
-
-                      // at inner points, the gradient is zero for ind==col
-                      in0 = in[stride*ind];
-                      in1 = in[stride*(mm-1-ind)];
-                      if (ind == col)
-                        {
-                          res1 += val1 * in0;
-                          res0 -= val1 * in1;
-                        }
-                      else
-                        {
-                          res0 += val0 * in0;
-                          res1 += val1 * in0;
-                          res0 -= val1 * in1;
-                          res1 -= val0 * in1;
-                        }
-                    }
-                }
-              else
-                res0 = res1 = VectorizedArray<Number>();
-              if (mm % 2 == 1)
-                {
-                  if (dof_to_quad == true)
-                    val0 = this->data.shape_gradients[mid*mm+col];
-                  else
-                    val0 = this->data.shape_gradients[col*mm+mid];
-                  val1 = val0 * in[stride*mid];
-                  res0 += val1;
-                  res1 -= val1;
-                }
-              if (add == false)
-                {
-                  out[stride*col]         = res0;
-                  out[stride*(nn-1-col)]  = res1;
-                }
-              else
-                {
-                  out[stride*col]        += res0;
-                  out[stride*(nn-1-col)] += res1;
-                }
-            }
-          if ( nn%2 == 1 )
-            {
-              VectorizedArray<Number> val0, res0;
-              if (dof_to_quad == true)
-                val0 = this->data.shape_gradients[n_cols];
-              else
-                val0 = this->data.shape_gradients[n_cols*mm];
-              if (mid > 0)
-                {
-                  res0  = in[0] - in[stride*(mm-1)];
-                  res0 *= val0;
-                  for (int ind=1; ind<mid; ++ind)
-                    {
-                      if (dof_to_quad == true)
-                        val0 = this->data.shape_gradients[ind*mm+n_cols];
-                      else
-                        val0 = this->data.shape_gradients[n_cols*mm+ind];
-                      VectorizedArray<Number> val1  = in[stride*ind] - in[stride*(mm-1-ind)];
-                      val1 *= val0;
-                      res0 += val1;
-                    }
-                }
-              else
-                res0 = VectorizedArray<Number>();
-              if (add == false)
-                out[stride*n_cols]  = res0;
-              else
-                out[stride*n_cols] += res0;
-            }
-
-          // increment: in regular case, just go to the next point in
-          // x-direction. for y-part in 3D and if we are at the end of one
-          // chunk in x-dir, need to jump over to the next layer in
-          // z-direction
-          switch (direction)
-            {
-            case 0:
-              in += mm;
-              out += nn;
-              break;
-            case 1:
-            case 2:
-              ++in;
-              ++out;
-              break;
-            default:
-              Assert (false, ExcNotImplemented());
-            }
-        }
+  internal::apply_tensor_product_gradients_gl<dim,fe_degree,Number,
+           direction, dof_to_quad, add>
+           (this->data.shape_gradients.begin(), in, out);
+}
 
-      if (direction == 1)
-        {
-          in  += nn * (mm-1);
-          out += nn * (nn-1);
-        }
-    }
 
 #endif  // ifndef DOXYGEN
 
-}
-
 
 DEAL_II_NAMESPACE_CLOSE
 
index 7594dc2750a4cf6b44e32db262f9b4fbb06079eb..588f1af187cc150c9a0128b5e9ce7e8f3599a933 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2008, 2009, 2010, 2011, 2012 by the deal.II authors
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2008, 2009, 2010, 2011, 2012, 2013 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -57,10 +57,7 @@ LogStream::LogStream()
   float_threshold(0.),
   offset(0),
   old_cerr(0),
-  at_newline(true),
-  stream_flags(std::ios::showpoint | std::ios::left),
-  stream_width(std::cout.width()),
-  stream_precision(std::cout.precision())
+  at_newline(true)
 {
   get_prefixes().push("DEAL:");
 
@@ -134,10 +131,6 @@ LogStream::operator<< (std::ostream& (*p) (std::ostream &))
 {
 
   std::ostringstream &stream = get_stream();
-  // save the state of out stream
-  std::ios::fmtflags old_flags = stream.flags(stream_flags);
-  unsigned int old_precision = stream.precision (stream_precision);
-  unsigned int old_width = stream.width (stream_width);
 
   // Print to the internal stringstream:
   stream << p;
@@ -170,11 +163,6 @@ LogStream::operator<< (std::ostream& (*p) (std::ostream &))
       stream.str("");
     }
 
-  // reset output format
-  stream.flags (old_flags);
-  stream.precision(old_precision);
-  stream.width(old_width);
-
   return *this;
 }
 
@@ -260,33 +248,21 @@ void LogStream::pop ()
 std::ios::fmtflags
 LogStream::flags(const std::ios::fmtflags f)
 {
-  Threads::Mutex::ScopedLock lock(log_lock);
-
-  std::ios::fmtflags tmp = stream_flags;
-  stream_flags = f;
-  return tmp;
+  return get_stream().flags (f);
 }
 
 
 std::streamsize
 LogStream::precision (const std::streamsize prec)
 {
-  Threads::Mutex::ScopedLock lock(log_lock);
-
-  std::streamsize tmp = stream_precision;
-  stream_precision = prec;
-  return tmp;
+  return get_stream().precision (prec);
 }
 
 
 std::streamsize
 LogStream::width (const std::streamsize wide)
 {
-  Threads::Mutex::ScopedLock lock(log_lock);
-
-  std::streamsize tmp = stream_width;
-  stream_width = wide;
-  return tmp;
+  return get_stream().width (wide);
 }
 
 

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.