]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Do symmetric tensor product evaluation with even-odd formula.
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 22 Oct 2013 09:50:02 +0000 (09:50 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 22 Oct 2013 09:50:02 +0000 (09:50 +0000)
git-svn-id: https://svn.dealii.org/trunk@31380 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/matrix_free/fe_evaluation.h
tests/matrix_free/evaluate_1d_shape.cc [new file with mode: 0644]
tests/matrix_free/evaluate_1d_shape.output [new file with mode: 0644]
tests/matrix_free/evaluate_1d_shape_evenodd.cc [new file with mode: 0644]
tests/matrix_free/evaluate_1d_shape_evenodd.output [new file with mode: 0644]

index edcb3095cfbe06b04af27e36b12c5222c12116d4..aa3ff08f554c2904e52d378ddb9072baef3717db 100644 (file)
@@ -1304,6 +1304,10 @@ protected:
   void apply_hessians (const VectorizedArray<Number> in [],
                        VectorizedArray<Number> out []);
 
+  VectorizedArray<Number> shape_val_evenodd[fe_degree+1][(n_q_points_1d+1)/2];
+  VectorizedArray<Number> shape_gra_evenodd[fe_degree+1][(n_q_points_1d+1)/2];
+  VectorizedArray<Number> shape_hes_evenodd[fe_degree+1][(n_q_points_1d+1)/2];
+
   /**
    * Friend declarations.
    */
@@ -4390,6 +4394,204 @@ namespace internal
 
 
 
+  // This method implements a different approach to the symmetric case for
+  // values, gradients, and Hessians also treated with the above functions: It
+  // is possible to reduce the cost per dimension from N^2 to N^2/2, where N
+  // is the number of 1D dofs (there are only N^2/2 different entries in the
+  // shape matrix, so this is plausible). The approach is based on the idea of
+  // applying the operator on the even and odd part of the input vectors
+  // separately, given that the shape functions evaluated on quadrature points
+  // are symmetric. This method is presented e.g. in the book "Implementing
+  // Spectral Methods for Partial Differential Equations" by David A. Kopriva,
+  // Springer, 2009, section 3.5.3 (Even-Odd-Decomposition). Even though the
+  // experiments in the book say that the method is not efficient for N<20, it
+  // is more efficient in the context where the loop bounds are compile-time
+  // constants (templates).
+  template <int dim, int fe_degree, int n_q_points_1d, typename Number,
+            int direction, bool dof_to_quad, bool add, int type>
+  inline
+  void
+  apply_tensor_product_evenodd (const Number shapes [][(n_q_points_1d+1)/2],
+                                const Number in [],
+                                Number       out [])
+  {
+    AssertIndexRange (type, 3);
+    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 code may look very inefficient at first sight due to the many
+    // different cases with if's at the innermost loop part, but all of the
+    // conditionals can be evaluated at compile time because they are
+    // templates, so the compiler should optimize everything away
+    for (int i2=0; i2<n_blocks2; ++i2)
+      {
+        for (int i1=0; i1<n_blocks1; ++i1)
+          {
+            Number xp[mid], xm[mid];
+            for (int i=0; i<mid; ++i)
+              {
+                if (dof_to_quad == true && type == 1)
+                  {
+                    xp[i] = in[stride*i] - in[stride*(mm-1-i)];
+                    xm[i] = in[stride*i] + in[stride*(mm-1-i)];
+                  }
+                else
+                  {
+                    xp[i] = in[stride*i] + in[stride*(mm-1-i)];
+                    xm[i] = in[stride*i] - in[stride*(mm-1-i)];
+                  }
+              }
+            for (int col=0; col<n_cols; ++col)
+              {
+                Number r0, r1;
+                if (mid > 0)
+                  {
+                    if (dof_to_quad == true)
+                      {
+                        r0 = shapes[0][col]         * xp[0];
+                        r1 = shapes[fe_degree][col] * xm[0];
+                      }
+                    else
+                      {
+                        r0 = shapes[col][0]           * xp[0];
+                        r1 = shapes[fe_degree-col][0] * xm[0];
+                      }
+                    for (int ind=1; ind<mid; ++ind)
+                      {
+                        if (dof_to_quad == true)
+                          {
+                            r0 += shapes[ind][col]           * xp[ind];
+                            r1 += shapes[fe_degree-ind][col] * xm[ind];
+                          }
+                        else
+                          {
+                            r0 += shapes[col][ind]           * xp[ind];
+                            r1 += shapes[fe_degree-col][ind] * xm[ind];
+                          }
+                      }
+                  }
+                else
+                  r0 = r1 = Number();
+                if (mm % 2 == 1 && dof_to_quad == true)
+                  {
+                    if (type == 1)
+                      r1 += shapes[mid][col] * in[stride*mid];
+                    else
+                      r0 += shapes[mid][col] * in[stride*mid];
+                  }
+                else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0))
+                  r0 += shapes[col][mid] * in[stride*mid];
+
+                if (add == false)
+                  {
+                    out[stride*col]         = r0 + r1;
+                    if (type == 1 && dof_to_quad == false)
+                      out[stride*(nn-1-col)]  = r1 - r0;
+                    else
+                      out[stride*(nn-1-col)]  = r0 - r1;
+                  }
+                else
+                  {
+                    out[stride*col]        += r0 + r1;
+                    if (type == 1 && dof_to_quad == false)
+                      out[stride*(nn-1-col)] += r1 - r0;
+                    else
+                      out[stride*(nn-1-col)] += r0 - r1;
+                  }
+              }
+            if ( type == 0 && 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)
+              {
+                Number r0;
+                if (mid > 0)
+                  {
+                    r0  = shapes[0][n_cols] * xp[0];
+                    for (int ind=1; ind<mid; ++ind)
+                      r0 += shapes[ind][n_cols] * xp[ind];
+                  }
+                else
+                  r0 = Number();
+                if (type != 1 && mm % 2 == 1)
+                  r0 += shapes[mid][n_cols] * in[stride*mid];
+
+                if (add == false)
+                  out[stride*n_cols]  = r0;
+                else
+                  out[stride*n_cols] += r0;
+              }
+            else if (dof_to_quad == false && nn%2 == 1)
+              {
+                Number r0;
+                if (mid > 0)
+                  {
+                    if (type == 1)
+                      {
+                        r0 = shapes[n_cols][0] * xm[0];
+                        for (int ind=1; ind<mid; ++ind)
+                          r0 += shapes[n_cols][ind] * xm[ind];
+                      }
+                    else
+                      {
+                        r0 = shapes[n_cols][0] * xp[0];
+                        for (int ind=1; ind<mid; ++ind)
+                          r0 += shapes[n_cols][ind] * xp[ind];
+                      }
+                  }
+                else
+                  r0 = Number();
+
+                if (type == 0 && mm % 2 == 1)
+                  r0 += in[stride*mid];
+                else if (type == 2 && mm % 2 == 1)
+                  r0 += shapes[n_cols][mid] * in[stride*mid];
+
+                if (add == false)
+                  out[stride*n_cols]  = r0;
+                else
+                  out[stride*n_cols] += r0;
+              }
+
+            // 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
@@ -4904,7 +5106,8 @@ namespace internal
     fe_eval.dof_values_initialized = true;
 #endif
   }
-}
+
+} // end of namespace internal
 
 
 
@@ -5090,6 +5293,43 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
                                                   j-1][0]) < zero_tol,
               ExcMessage(error_message));
 #endif
+
+  // Compute symmetric and skew-symmetric part of shape values for even-odd
+  // decomposition
+  for (unsigned int i=0; i<(fe_degree+1)/2; ++i)
+    for (unsigned int q=0; q<(n_q_points_1d+1)/2; ++q)
+      {
+        shape_val_evenodd[i][q] =
+          0.5 * (this->data.shape_values[i*n_q_points_1d+q] +
+                 this->data.shape_values[i*n_q_points_1d+n_q_points_1d-1-q]);
+        shape_val_evenodd[fe_degree-i][q] =
+          0.5 * (this->data.shape_values[i*n_q_points_1d+q] -
+                 this->data.shape_values[i*n_q_points_1d+n_q_points_1d-1-q]);
+
+        shape_gra_evenodd[i][q] =
+          0.5 * (this->data.shape_gradients[i*n_q_points_1d+q] +
+                 this->data.shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]);
+        shape_gra_evenodd[fe_degree-i][q] =
+          0.5 * (this->data.shape_gradients[i*n_q_points_1d+q] -
+                 this->data.shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]);
+
+        shape_hes_evenodd[i][q] =
+          0.5 * (this->data.shape_hessians[i*n_q_points_1d+q] +
+                 this->data.shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]);
+        shape_hes_evenodd[fe_degree-i][q] =
+          0.5 * (this->data.shape_hessians[i*n_q_points_1d+q] -
+                 this->data.shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]);
+      }
+  if (fe_degree % 2 == 0)
+    for (unsigned int q=0; q<(n_q_points_1d+1)/2; ++q)
+      {
+        shape_val_evenodd[fe_degree/2][q] =
+          this->data.shape_values[(fe_degree/2)*n_q_points_1d+q];
+        shape_gra_evenodd[fe_degree/2][q] =
+          this->data.shape_gradients[(fe_degree/2)*n_q_points_1d+q];
+        shape_hes_evenodd[fe_degree/2][q] =
+          this->data.shape_hessians[(fe_degree/2)*n_q_points_1d+q];
+      }
 }
 
 
@@ -5129,9 +5369,9 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
 ::apply_values (const VectorizedArray<Number> in [],
                 VectorizedArray<Number>       out [])
 {
-  internal::apply_tensor_product_values<dim,fe_degree,n_q_points_1d,
-           VectorizedArray<Number>, direction, dof_to_quad, add>
-           (this->data.shape_values.begin(), in, out);
+  internal::apply_tensor_product_evenodd<dim,fe_degree,n_q_points_1d,
+           VectorizedArray<Number>, direction, dof_to_quad, add, 0>
+           (shape_val_evenodd, in, out);
 }
 
 
@@ -5145,9 +5385,9 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
 ::apply_gradients (const VectorizedArray<Number> in [],
                    VectorizedArray<Number>       out [])
 {
-  internal::apply_tensor_product_gradients<dim,fe_degree,n_q_points_1d,
-           VectorizedArray<Number>, direction, dof_to_quad, add>
-           (this->data.shape_gradients.begin(), in, out);
+  internal::apply_tensor_product_evenodd<dim,fe_degree,n_q_points_1d,
+           VectorizedArray<Number>, direction, dof_to_quad, add, 1>
+           (shape_gra_evenodd, in, out);
 }
 
 
@@ -5164,9 +5404,9 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
 ::apply_hessians (const VectorizedArray<Number> in [],
                   VectorizedArray<Number>       out [])
 {
-  internal::apply_tensor_product_hessians<dim,fe_degree,n_q_points_1d,
-           VectorizedArray<Number>, direction, dof_to_quad, add>
-           (this->data.shape_hessians.begin(), in, out);
+  internal::apply_tensor_product_evenodd<dim,fe_degree,n_q_points_1d,
+           VectorizedArray<Number>, direction, dof_to_quad, add, 2>
+           (shape_hes_evenodd, in, out);
 }
 
 
diff --git a/tests/matrix_free/evaluate_1d_shape.cc b/tests/matrix_free/evaluate_1d_shape.cc
new file mode 100644 (file)
index 0000000..b19a94c
--- /dev/null
@@ -0,0 +1,174 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// this function tests the correctness of the 1d evaluation functions used in
+// FEEvaluation. These functions are marked 'internal' but it is much easier
+// to check their correctness directly rather than from the results in
+// dependent functions
+
+#include "../tests.h"
+#include <iostream>
+#include <fstream>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+
+template <int M, int N, int type, bool add>
+void test()
+{
+  deallog << "Test " << M << " x " << N << std::endl;
+  double shape[M][N];
+  for (unsigned int i=0; i<(M+1)/2; ++i)
+    for (unsigned int j=0; j<N; ++j)
+      {
+        shape[i][j] = -1. + 2. * (double)rand()/RAND_MAX;
+        if (type == 1)
+          shape[M-1-i][N-1-j] = -shape[i][j];
+        else
+          shape[M-1-i][N-1-j] = shape[i][j];
+      }
+  if (type == 0 && M%2 == 1 && N%2 == 1)
+    {
+      for (unsigned int i=0; i<M; ++i)
+        shape[i][N/2] = 0.;
+      shape[M/2][N/2] = 1;
+    }
+  if (type == 1 && M%2 == 1 && N%2 == 1)
+    shape[M/2][N/2] = 0.;
+
+  double x[N], x_ref[N], y[M], y_ref[M];
+  for (unsigned int i=0; i<N; ++i)
+    x[i] = (double)rand()/RAND_MAX;
+
+  // compute reference
+  for (unsigned int i=0; i<M; ++i)
+    {
+      y[i] = 1.;
+      y_ref[i] = add ? y[i] : 0.;
+      for (unsigned int j=0; j<N; ++j)
+        y_ref[i] += shape[i][j] * x[j];
+    }
+
+  // apply function for tensor product
+  if (type == 0)
+    internal::apply_tensor_product_values<1,M-1,N,double,0,false,add>
+      (&shape[0][0],x,y);
+  if (type == 1)
+    internal::apply_tensor_product_gradients<1,M-1,N,double,0,false,add>
+      (&shape[0][0],x,y);
+  if (type == 2)
+    internal::apply_tensor_product_hessians<1,M-1,N,double,0,false,add>
+      (&shape[0][0],x,y);
+
+  deallog << "Errors no transpose: ";
+  for (unsigned int i=0; i<M; ++i)
+    deallog << y[i] - y_ref[i] << " ";
+  deallog << std::endl;
+
+
+  for (unsigned int i=0; i<M; ++i)
+    y[i] = (double)rand()/RAND_MAX;
+
+  // compute reference
+  for (unsigned int i=0; i<N; ++i)
+    {
+      x[i] = 2.;
+      x_ref[i] = add ? x[i] : 0.;
+      for (unsigned int j=0; j<M; ++j)
+        x_ref[i] += shape[j][i] * y[j];
+    }
+
+  // apply function for tensor product
+  if (type == 0)
+    internal::apply_tensor_product_values<1,M-1,N,double,0,true,add>
+      (&shape[0][0],y,x);
+  if (type == 1)
+    internal::apply_tensor_product_gradients<1,M-1,N,double,0,true,add>
+      (&shape[0][0],y,x);
+  if (type == 2)
+    internal::apply_tensor_product_hessians<1,M-1,N,double,0,true,add>
+      (&shape[0][0],y,x);
+
+  deallog << "Errors transpose:    ";
+  for (unsigned int i=0; i<N; ++i)
+    deallog << x[i] - x_ref[i] << " ";
+  deallog << std::endl;
+}
+
+int main ()
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1e-14);
+
+  deallog.push("values");
+  test<4,4,0,false>();
+  test<3,3,0,false>();
+  test<4,3,0,false>();
+  test<3,4,0,false>();
+  test<3,5,0,false>();
+  deallog.pop();
+
+  deallog.push("gradients");
+  test<4,4,1,false>();
+  test<3,3,1,false>();
+  test<4,3,1,false>();
+  test<3,4,1,false>();
+  test<3,5,1,false>();
+  deallog.pop();
+
+  deallog.push("hessians");
+  test<4,4,2,false>();
+  test<3,3,2,false>();
+  test<4,3,2,false>();
+  test<3,4,2,false>();
+  test<3,5,2,false>();
+  deallog.pop();
+
+  deallog.push("add");
+
+  deallog.push("values");
+  test<4,4,0,true>();
+  test<3,3,0,true>();
+  test<4,3,0,true>();
+  test<3,4,0,true>();
+  test<3,5,0,true>();
+  deallog.pop();
+
+  deallog.push("gradients");
+  test<4,4,1,true>();
+  test<3,3,1,true>();
+  test<4,3,1,true>();
+  test<3,4,1,true>();
+  test<3,5,1,true>();
+  deallog.pop();
+
+  deallog.push("hessians");
+  test<4,4,2,true>();
+  test<3,3,2,true>();
+  test<4,3,2,true>();
+  test<3,4,2,true>();
+  test<3,5,2,true>();
+  deallog.pop();
+
+  deallog.pop();
+
+  return 0;
+}
+    
diff --git a/tests/matrix_free/evaluate_1d_shape.output b/tests/matrix_free/evaluate_1d_shape.output
new file mode 100644 (file)
index 0000000..2bd9530
--- /dev/null
@@ -0,0 +1,91 @@
+
+DEAL:values::Test 4 x 4
+DEAL:values::Errors no transpose: 0 0 0 0 
+DEAL:values::Errors transpose:    0 0 0 0 
+DEAL:values::Test 3 x 3
+DEAL:values::Errors no transpose: 0 0 0 
+DEAL:values::Errors transpose:    0 0 0 
+DEAL:values::Test 4 x 3
+DEAL:values::Errors no transpose: 0 0 0 0 
+DEAL:values::Errors transpose:    0 0 0 
+DEAL:values::Test 3 x 4
+DEAL:values::Errors no transpose: 0 0 0 
+DEAL:values::Errors transpose:    0 0 0 0 
+DEAL:values::Test 3 x 5
+DEAL:values::Errors no transpose: 0 0 0 
+DEAL:values::Errors transpose:    0 0 0 0 0 
+DEAL:gradients::Test 4 x 4
+DEAL:gradients::Errors no transpose: 0 0 0 0 
+DEAL:gradients::Errors transpose:    0 0 0 0 
+DEAL:gradients::Test 3 x 3
+DEAL:gradients::Errors no transpose: 0 0 0 
+DEAL:gradients::Errors transpose:    0 0 0 
+DEAL:gradients::Test 4 x 3
+DEAL:gradients::Errors no transpose: 0 0 0 0 
+DEAL:gradients::Errors transpose:    0 0 0 
+DEAL:gradients::Test 3 x 4
+DEAL:gradients::Errors no transpose: 0 0 0 
+DEAL:gradients::Errors transpose:    0 0 0 0 
+DEAL:gradients::Test 3 x 5
+DEAL:gradients::Errors no transpose: 0 0 0 
+DEAL:gradients::Errors transpose:    0 0 0 0 0 
+DEAL:hessians::Test 4 x 4
+DEAL:hessians::Errors no transpose: 0 0 0 0 
+DEAL:hessians::Errors transpose:    0 0 0 0 
+DEAL:hessians::Test 3 x 3
+DEAL:hessians::Errors no transpose: 0 0 0 
+DEAL:hessians::Errors transpose:    0 0 0 
+DEAL:hessians::Test 4 x 3
+DEAL:hessians::Errors no transpose: 0 0 0 0 
+DEAL:hessians::Errors transpose:    0 0 0 
+DEAL:hessians::Test 3 x 4
+DEAL:hessians::Errors no transpose: 0 0 0 
+DEAL:hessians::Errors transpose:    0 0 0 0 
+DEAL:hessians::Test 3 x 5
+DEAL:hessians::Errors no transpose: 0 0 0 
+DEAL:hessians::Errors transpose:    0 0 0 0 0 
+DEAL:add:values::Test 4 x 4
+DEAL:add:values::Errors no transpose: 0 0 0 0 
+DEAL:add:values::Errors transpose:    0 0 0 0 
+DEAL:add:values::Test 3 x 3
+DEAL:add:values::Errors no transpose: 0 0 0 
+DEAL:add:values::Errors transpose:    0 0 0 
+DEAL:add:values::Test 4 x 3
+DEAL:add:values::Errors no transpose: 0 0 0 0 
+DEAL:add:values::Errors transpose:    0 0 0 
+DEAL:add:values::Test 3 x 4
+DEAL:add:values::Errors no transpose: 0 0 0 
+DEAL:add:values::Errors transpose:    0 0 0 0 
+DEAL:add:values::Test 3 x 5
+DEAL:add:values::Errors no transpose: 0 0 0 
+DEAL:add:values::Errors transpose:    0 0 0 0 0 
+DEAL:add:gradients::Test 4 x 4
+DEAL:add:gradients::Errors no transpose: 0 0 0 0 
+DEAL:add:gradients::Errors transpose:    0 0 0 0 
+DEAL:add:gradients::Test 3 x 3
+DEAL:add:gradients::Errors no transpose: 0 0 0 
+DEAL:add:gradients::Errors transpose:    0 0 0 
+DEAL:add:gradients::Test 4 x 3
+DEAL:add:gradients::Errors no transpose: 0 0 0 0 
+DEAL:add:gradients::Errors transpose:    0 0 0 
+DEAL:add:gradients::Test 3 x 4
+DEAL:add:gradients::Errors no transpose: 0 0 0 
+DEAL:add:gradients::Errors transpose:    0 0 0 0 
+DEAL:add:gradients::Test 3 x 5
+DEAL:add:gradients::Errors no transpose: 0 0 0 
+DEAL:add:gradients::Errors transpose:    0 0 0 0 0 
+DEAL:add:hessians::Test 4 x 4
+DEAL:add:hessians::Errors no transpose: 0 0 0 0 
+DEAL:add:hessians::Errors transpose:    0 0 0 0 
+DEAL:add:hessians::Test 3 x 3
+DEAL:add:hessians::Errors no transpose: 0 0 0 
+DEAL:add:hessians::Errors transpose:    0 0 0 
+DEAL:add:hessians::Test 4 x 3
+DEAL:add:hessians::Errors no transpose: 0 0 0 0 
+DEAL:add:hessians::Errors transpose:    0 0 0 
+DEAL:add:hessians::Test 3 x 4
+DEAL:add:hessians::Errors no transpose: 0 0 0 
+DEAL:add:hessians::Errors transpose:    0 0 0 0 
+DEAL:add:hessians::Test 3 x 5
+DEAL:add:hessians::Errors no transpose: 0 0 0 
+DEAL:add:hessians::Errors transpose:    0 0 0 0 0 
diff --git a/tests/matrix_free/evaluate_1d_shape_evenodd.cc b/tests/matrix_free/evaluate_1d_shape_evenodd.cc
new file mode 100644 (file)
index 0000000..2eab664
--- /dev/null
@@ -0,0 +1,172 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// this function tests the correctness of the 1d evaluation functions used in
+// FEEvaluation. These functions are marked 'internal' but it is much easier
+// to check their correctness directly rather than from the results in
+// dependent functions. this function tests the even-odd path of the
+// evaluation functions
+
+#include "../tests.h"
+#include <iostream>
+#include <fstream>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+
+template <int M, int N, int type, bool add>
+void test()
+{
+  deallog << "Test " << M << " x " << N << std::endl;
+  double shape[M][N];
+  for (unsigned int i=0; i<(M+1)/2; ++i)
+    for (unsigned int j=0; j<N; ++j)
+      {
+        shape[i][j] = -1. + 2. * (double)rand()/RAND_MAX;
+        if (type == 1)
+          shape[M-1-i][N-1-j] = -shape[i][j];
+        else
+          shape[M-1-i][N-1-j] = shape[i][j];
+      }
+  if (type == 0 && M%2 == 1 && N%2 == 1)
+    {
+      for (unsigned int i=0; i<M; ++i)
+        shape[i][N/2] = 0.;
+      shape[M/2][N/2] = 1;
+    }
+  if (type == 1 && M%2 == 1 && N%2 == 1)
+    shape[M/2][N/2] = 0.;
+
+  // create symmetrized shape array exactly as expected by the evenodd
+  // function
+  double shape_sym[M][(N+1)/2];
+  for (unsigned int i=0; i<M/2; ++i)
+    for (unsigned int q=0; q<(N+1)/2; ++q)
+      {
+        shape_sym[i][q] = 0.5 * (shape[i][q] + shape[i][N-1-q]);
+        shape_sym[M-1-i][q] = 0.5 * (shape[i][q] - shape[i][N-1-q]);
+      }
+  if (M % 2 == 1)
+    for (unsigned int q=0; q<(N+1)/2; ++q)
+      shape_sym[(M-1)/2][q] = shape[(M-1)/2][q];
+
+  double x[N], x_ref[N], y[M], y_ref[M];
+  for (unsigned int i=0; i<N; ++i)
+    x[i] = (double)rand()/RAND_MAX;
+
+  // compute reference
+  for (unsigned int i=0; i<M; ++i)
+    {
+      y[i] = 1.;
+      y_ref[i] = add ? y[i] : 0.;
+      for (unsigned int j=0; j<N; ++j)
+        y_ref[i] += shape[i][j] * x[j];
+    }
+
+  // apply function for tensor product
+  internal::apply_tensor_product_evenodd<1,M-1,N,double,0,false,add,type>(shape_sym,x,y);
+
+  deallog << "Errors no transpose: ";
+  for (unsigned int i=0; i<M; ++i)
+    deallog << y[i] - y_ref[i] << " ";
+  deallog << std::endl;
+
+
+  for (unsigned int i=0; i<M; ++i)
+    y[i] = (double)rand()/RAND_MAX;
+
+  // compute reference
+  for (unsigned int i=0; i<N; ++i)
+    {
+      x[i] = 2.;
+      x_ref[i] = add ? x[i] : 0.;
+      for (unsigned int j=0; j<M; ++j)
+        x_ref[i] += shape[j][i] * y[j];
+    }
+
+  // apply function for tensor product
+  internal::apply_tensor_product_evenodd<1,M-1,N,double,0,true,add,type>(shape_sym,y,x);
+
+  deallog << "Errors transpose:    ";
+  for (unsigned int i=0; i<N; ++i)
+    deallog << x[i] - x_ref[i] << " ";
+  deallog << std::endl;
+}
+
+int main ()
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1e-14);
+
+  deallog.push("values");
+  test<4,4,0,false>();
+  test<3,3,0,false>();
+  test<4,3,0,false>();
+  test<3,4,0,false>();
+  test<3,5,0,false>();
+  deallog.pop();
+
+  deallog.push("gradients");
+  test<4,4,1,false>();
+  test<3,3,1,false>();
+  test<4,3,1,false>();
+  test<3,4,1,false>();
+  test<3,5,1,false>();
+  deallog.pop();
+
+  deallog.push("hessians");
+  test<4,4,2,false>();
+  test<3,3,2,false>();
+  test<4,3,2,false>();
+  test<3,4,2,false>();
+  test<3,5,2,false>();
+  deallog.pop();
+
+  deallog.push("add");
+
+  deallog.push("values");
+  test<4,4,0,true>();
+  test<3,3,0,true>();
+  test<4,3,0,true>();
+  test<3,4,0,true>();
+  test<3,5,0,true>();
+  deallog.pop();
+
+  deallog.push("gradients");
+  test<4,4,1,true>();
+  test<3,3,1,true>();
+  test<4,3,1,true>();
+  test<3,4,1,true>();
+  test<3,5,1,true>();
+  deallog.pop();
+
+  deallog.push("hessians");
+  test<4,4,2,true>();
+  test<3,3,2,true>();
+  test<4,3,2,true>();
+  test<3,4,2,true>();
+  test<3,5,2,true>();
+  deallog.pop();
+
+  deallog.pop();
+
+  return 0;
+}
+    
diff --git a/tests/matrix_free/evaluate_1d_shape_evenodd.output b/tests/matrix_free/evaluate_1d_shape_evenodd.output
new file mode 100644 (file)
index 0000000..2bd9530
--- /dev/null
@@ -0,0 +1,91 @@
+
+DEAL:values::Test 4 x 4
+DEAL:values::Errors no transpose: 0 0 0 0 
+DEAL:values::Errors transpose:    0 0 0 0 
+DEAL:values::Test 3 x 3
+DEAL:values::Errors no transpose: 0 0 0 
+DEAL:values::Errors transpose:    0 0 0 
+DEAL:values::Test 4 x 3
+DEAL:values::Errors no transpose: 0 0 0 0 
+DEAL:values::Errors transpose:    0 0 0 
+DEAL:values::Test 3 x 4
+DEAL:values::Errors no transpose: 0 0 0 
+DEAL:values::Errors transpose:    0 0 0 0 
+DEAL:values::Test 3 x 5
+DEAL:values::Errors no transpose: 0 0 0 
+DEAL:values::Errors transpose:    0 0 0 0 0 
+DEAL:gradients::Test 4 x 4
+DEAL:gradients::Errors no transpose: 0 0 0 0 
+DEAL:gradients::Errors transpose:    0 0 0 0 
+DEAL:gradients::Test 3 x 3
+DEAL:gradients::Errors no transpose: 0 0 0 
+DEAL:gradients::Errors transpose:    0 0 0 
+DEAL:gradients::Test 4 x 3
+DEAL:gradients::Errors no transpose: 0 0 0 0 
+DEAL:gradients::Errors transpose:    0 0 0 
+DEAL:gradients::Test 3 x 4
+DEAL:gradients::Errors no transpose: 0 0 0 
+DEAL:gradients::Errors transpose:    0 0 0 0 
+DEAL:gradients::Test 3 x 5
+DEAL:gradients::Errors no transpose: 0 0 0 
+DEAL:gradients::Errors transpose:    0 0 0 0 0 
+DEAL:hessians::Test 4 x 4
+DEAL:hessians::Errors no transpose: 0 0 0 0 
+DEAL:hessians::Errors transpose:    0 0 0 0 
+DEAL:hessians::Test 3 x 3
+DEAL:hessians::Errors no transpose: 0 0 0 
+DEAL:hessians::Errors transpose:    0 0 0 
+DEAL:hessians::Test 4 x 3
+DEAL:hessians::Errors no transpose: 0 0 0 0 
+DEAL:hessians::Errors transpose:    0 0 0 
+DEAL:hessians::Test 3 x 4
+DEAL:hessians::Errors no transpose: 0 0 0 
+DEAL:hessians::Errors transpose:    0 0 0 0 
+DEAL:hessians::Test 3 x 5
+DEAL:hessians::Errors no transpose: 0 0 0 
+DEAL:hessians::Errors transpose:    0 0 0 0 0 
+DEAL:add:values::Test 4 x 4
+DEAL:add:values::Errors no transpose: 0 0 0 0 
+DEAL:add:values::Errors transpose:    0 0 0 0 
+DEAL:add:values::Test 3 x 3
+DEAL:add:values::Errors no transpose: 0 0 0 
+DEAL:add:values::Errors transpose:    0 0 0 
+DEAL:add:values::Test 4 x 3
+DEAL:add:values::Errors no transpose: 0 0 0 0 
+DEAL:add:values::Errors transpose:    0 0 0 
+DEAL:add:values::Test 3 x 4
+DEAL:add:values::Errors no transpose: 0 0 0 
+DEAL:add:values::Errors transpose:    0 0 0 0 
+DEAL:add:values::Test 3 x 5
+DEAL:add:values::Errors no transpose: 0 0 0 
+DEAL:add:values::Errors transpose:    0 0 0 0 0 
+DEAL:add:gradients::Test 4 x 4
+DEAL:add:gradients::Errors no transpose: 0 0 0 0 
+DEAL:add:gradients::Errors transpose:    0 0 0 0 
+DEAL:add:gradients::Test 3 x 3
+DEAL:add:gradients::Errors no transpose: 0 0 0 
+DEAL:add:gradients::Errors transpose:    0 0 0 
+DEAL:add:gradients::Test 4 x 3
+DEAL:add:gradients::Errors no transpose: 0 0 0 0 
+DEAL:add:gradients::Errors transpose:    0 0 0 
+DEAL:add:gradients::Test 3 x 4
+DEAL:add:gradients::Errors no transpose: 0 0 0 
+DEAL:add:gradients::Errors transpose:    0 0 0 0 
+DEAL:add:gradients::Test 3 x 5
+DEAL:add:gradients::Errors no transpose: 0 0 0 
+DEAL:add:gradients::Errors transpose:    0 0 0 0 0 
+DEAL:add:hessians::Test 4 x 4
+DEAL:add:hessians::Errors no transpose: 0 0 0 0 
+DEAL:add:hessians::Errors transpose:    0 0 0 0 
+DEAL:add:hessians::Test 3 x 3
+DEAL:add:hessians::Errors no transpose: 0 0 0 
+DEAL:add:hessians::Errors transpose:    0 0 0 
+DEAL:add:hessians::Test 4 x 3
+DEAL:add:hessians::Errors no transpose: 0 0 0 0 
+DEAL:add:hessians::Errors transpose:    0 0 0 
+DEAL:add:hessians::Test 3 x 4
+DEAL:add:hessians::Errors no transpose: 0 0 0 
+DEAL:add:hessians::Errors transpose:    0 0 0 0 
+DEAL:add:hessians::Test 3 x 5
+DEAL:add:hessians::Errors no transpose: 0 0 0 
+DEAL:add:hessians::Errors transpose:    0 0 0 0 0 

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.