Adjust some fields for face integration

git-svn-id: https://svn.dealii.org/trunk@31495 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/include/deal.II/matrix_free/fe_evaluation.h b/deal.II/include/deal.II/matrix_free/fe_evaluation.h
index aa3ff08f554c2904e52d378ddb9072baef3717db..7a18e32fd236874d612ad384c532c5ccc27ed321 100644 (file)
--- a/deal.II/include/deal.II/matrix_free/fe_evaluation.h
+++ b/deal.II/include/deal.II/matrix_free/fe_evaluation.h
@@ -3841,9 +3841,9 @@ namespace internal
                     res0 += val0 * in[stride*ind];
                   }
                 if (add == false)
-                  out[stride*col]         = res0;
+                  out[stride*col]  = res0;
                 else
-                  out[stride*col]        += res0;
+                  out[stride*col] += res0;
               }
 
             // increment: in regular case, just go to the next point in
@@ -3874,6 +3874,79 @@ namespace internal
 
 
 
+  // This method applies the tensor product operation to produce face values
+  // out from cell values. As opposed to the apply_tensor_product method, this
+  // method assumes that the directions orthogonal to the face have
+  // fe_degree+1 degrees of freedom per direction and not n_q_points_1d for
+  // those directions lower than the one currently applied
+  template <int dim, int fe_degree, typename Number, int face_direction,
+           bool dof_to_quad, bool add>
+  inline
+  void
+  apply_tensor_product_face (const Number *shape_data,
+                             const Number in [],
+                             Number       out [])
+  {
+    const int n_blocks1 = dim > 1 ? (fe_degree+1) : 1;
+    const int n_blocks2 = dim > 2 ? (fe_degree+1) : 1;
+
+    AssertIndexRange (face_direction, dim);
+    const int mm     = dof_to_quad ? (fe_degree+1) : 1,
+              nn     = dof_to_quad ? 1 : (fe_degree+1);
+
+    const int stride = Utilities::fixed_int_power<fe_degree+1,face_direction>::value;
+
+    for (int i2=0; i2<n_blocks2; ++i2)
+      {
+        for (int i1=0; i1<n_blocks1; ++i1)
+          {
+            if (dof_to_quad == true)
+              {
+                Number res0 = shape_data[0] * in[0];
+                for (int ind=1; ind<mm; ++ind)
+                  res0 += shape_data[ind] * in[stride*ind];
+                if (add == false)
+                  out[0]  = res0;
+                else
+                  out[0] += res0;
+              }
+            else
+              {
+                for (int col=0; col<nn; ++col)
+                  if (add == false)
+                    out[col*stride]  = shape_data[col] * in[0];
+                  else
+                    out[col*stride] += shape_data[col] * in[0];
+              }
+
+            // 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 (face_direction)
+              {
+              case 0:
+                in += mm;
+                out += nn;
+                break;
+              case 1:
+              case 2:
+                ++in;
+                ++out;
+                break;
+              default:
+                Assert (false, ExcNotImplemented());
+              }
+          }
+        if (face_direction == 1)
+          {
+            in += mm*(mm-1);
+            out += nn*(nn-1);
+          }
+      }
+  }
+
+
+
   // 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,

diff --git a/deal.II/include/deal.II/matrix_free/shape_info.h b/deal.II/include/deal.II/matrix_free/shape_info.h
index ea987f2d253a51926684baafb4640e51d7d3d295..e25b8bfb37e4fc371f636f0f7f72871947daa724 100644 (file)
--- a/deal.II/include/deal.II/matrix_free/shape_info.h
+++ b/deal.II/include/deal.II/matrix_free/shape_info.h
@@ -115,21 +115,23 @@ namespace internal
       Table<2,unsigned int>  face_indices;
 
       /**
-       * Stores one-dimensional values of shape
-       * functions on subface. Since there are two
-       * subfaces, store two variants. Not
-       * vectorized.
+       * Stores one-dimensional values of shape functions evaluated in zero
+       * and one, i.e., on the one-dimensional faces. Not vectorized.
        */
       std::vector<Number>    face_value[2];
 
       /**
-       * Stores one-dimensional gradients of shape
-       * functions on subface. Since there are two
-       * subfaces, store two variants. Not
-       * vectorized.
+       * Stores one-dimensional gradients of shape functions evaluated in zero
+       * and one, i.e., on the one-dimensional faces. Not vectorized.
        */
       std::vector<Number>    face_gradient[2];
 
+      /**
+       * Stores one-dimensional values of shape functions on subface. Since
+       * there are two subfaces, store two variants. Not vectorized.
+       */
+      std::vector<Number>    subface_value[2];
+
       /**
        * Non-vectorized version of shape
        * values. Needed when evaluating face info.

diff --git a/deal.II/include/deal.II/matrix_free/shape_info.templates.h b/deal.II/include/deal.II/matrix_free/shape_info.templates.h
index 1f9147f828c242dc7ad932887fb7804d7470ea05..852a6369a0b99e0187e7cd66473d0bb6ab0def15 100644 (file)
--- a/deal.II/include/deal.II/matrix_free/shape_info.templates.h
+++ b/deal.II/include/deal.II/matrix_free/shape_info.templates.h
@@ -91,10 +91,12 @@ namespace internal
       this->shape_values.resize_fast (array_size);
       this->shape_hessians.resize_fast (array_size);
 
+      this->face_value[0].resize(n_dofs_1d);
       this->face_gradient[0].resize(n_dofs_1d);
-      this->face_value[0].resize(array_size);
+      this->subface_value[0].resize(array_size);
+      this->face_value[1].resize(n_dofs_1d);
       this->face_gradient[1].resize(n_dofs_1d);
-      this->face_value[1].resize(array_size);
+      this->subface_value[1].resize(array_size);
       this->shape_values_number.resize (array_size);
       this->shape_gradient_number.resize (array_size);
 
@@ -120,13 +122,15 @@ namespace internal
               shape_hessians[i*n_q_points_1d+q] =
                 fe.shape_grad_grad(my_i,q_point)[0][0];
               q_point[0] *= 0.5;
-              face_value[0][i*n_q_points_1d+q] = fe.shape_value(my_i,q_point);
+              subface_value[0][i*n_q_points_1d+q] = fe.shape_value(my_i,q_point);
               q_point[0] += 0.5;
-              face_value[1][i*n_q_points_1d+q] = fe.shape_value(my_i,q_point);
+              subface_value[1][i*n_q_points_1d+q] = fe.shape_value(my_i,q_point);
             }
           Point<dim> q_point;
+          this->face_value[0][i] = fe.shape_value(my_i,q_point);
           this->face_gradient[0][i] = fe.shape_grad(my_i,q_point)[0];
           q_point[0] = 1;
+          this->face_value[1][i] = fe.shape_value(my_i,q_point);
           this->face_gradient[1][i] = fe.shape_grad(my_i,q_point)[0];
         }