--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii__matrix_free_evaluation_kernels_h
+#define dealii__matrix_free_evaluation_kernels_h
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/vectorization.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/matrix_free/tensor_product_kernels.h>
+#include <deal.II/matrix_free/shape_info.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+
+namespace internal
+{
+ // Select evaluator type from element shape function type
+ template <MatrixFreeFunctions::ElementType element, bool is_long>
+ struct EvaluatorSelector {};
+
+ template <bool is_long>
+ struct EvaluatorSelector<MatrixFreeFunctions::tensor_general,is_long>
+ {
+ static const EvaluatorVariant variant = evaluate_general;
+ };
+
+ template <>
+ struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric,false>
+ {
+ static const EvaluatorVariant variant = evaluate_symmetric;
+ };
+
+ template <> struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric,true>
+ {
+ static const EvaluatorVariant variant = evaluate_evenodd;
+ };
+
+ template <bool is_long>
+ struct EvaluatorSelector<MatrixFreeFunctions::truncated_tensor,is_long>
+ {
+ static const EvaluatorVariant variant = evaluate_general;
+ };
+
+ template <>
+ struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0,false>
+ {
+ static const EvaluatorVariant variant = evaluate_general;
+ };
+
+ template <>
+ struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0,true>
+ {
+ static const EvaluatorVariant variant = evaluate_evenodd;
+ };
+
+ template <bool is_long>
+ struct EvaluatorSelector<MatrixFreeFunctions::tensor_gausslobatto,is_long>
+ {
+ static const EvaluatorVariant variant = evaluate_evenodd;
+ };
+
+
+
+ // This struct performs the evaluation of function values, gradients and
+ // Hessians for tensor-product finite elements. The operation is used for
+ // both the symmetric and non-symmetric case, which use different apply
+ // functions 'values', 'gradients' in the individual coordinate
+ // directions. The apply functions for values are provided through one of
+ // the template classes EvaluatorTensorProduct which in turn are selected
+ // from the MatrixFreeFunctions::ElementType template argument.
+ //
+ // There is a specialization made for Gauss-Lobatto elements further down
+ // where the 'values' operation is identity, which allows us to write
+ // shorter code.
+ template <MatrixFreeFunctions::ElementType type, int dim, int fe_degree,
+ int n_q_points_1d, int n_components, typename Number>
+ struct FEEvaluationImpl
+ {
+ static
+ void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ VectorizedArray<Number> *values_dofs_actual[],
+ VectorizedArray<Number> *values_quad[],
+ VectorizedArray<Number> *gradients_quad[][dim],
+ VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
+ VectorizedArray<Number> *scratch_data,
+ const bool evaluate_val,
+ const bool evaluate_grad,
+ const bool evaluate_lapl);
+
+ static
+ void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ VectorizedArray<Number> *values_dofs_actual[],
+ VectorizedArray<Number> *values_quad[],
+ VectorizedArray<Number> *gradients_quad[][dim],
+ VectorizedArray<Number> *scratch_data,
+ const bool evaluate_val,
+ const bool evaluate_grad);
+ };
+
+
+ template <MatrixFreeFunctions::ElementType type, int dim, int fe_degree,
+ int n_q_points_1d, int n_components, typename Number>
+ inline
+ void
+ FEEvaluationImpl<type,dim,fe_degree,n_q_points_1d,n_components,Number>
+ ::evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ VectorizedArray<Number> *values_dofs_actual[],
+ VectorizedArray<Number> *values_quad[],
+ VectorizedArray<Number> *gradients_quad[][dim],
+ VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
+ VectorizedArray<Number> *scratch_data,
+ const bool evaluate_val,
+ const bool evaluate_grad,
+ const bool evaluate_lapl)
+ {
+ if (evaluate_val == false && evaluate_grad == false && evaluate_lapl == false)
+ return;
+
+ const EvaluatorVariant variant =
+ EvaluatorSelector<type,(fe_degree+n_q_points_1d>4)>::variant;
+ typedef EvaluatorTensorProduct<variant, dim, fe_degree, n_q_points_1d,
+ VectorizedArray<Number> > Eval;
+ Eval eval (variant == evaluate_evenodd ? shape_info.shape_val_evenodd :
+ shape_info.shape_values,
+ variant == evaluate_evenodd ? shape_info.shape_gra_evenodd :
+ shape_info.shape_gradients,
+ variant == evaluate_evenodd ? shape_info.shape_hes_evenodd :
+ shape_info.shape_hessians,
+ shape_info.fe_degree,
+ shape_info.n_q_points_1d);
+
+ const unsigned int temp_size = Eval::dofs_per_cell == numbers::invalid_unsigned_int ? 0
+ : (Eval::dofs_per_cell > Eval::n_q_points ?
+ Eval::dofs_per_cell : Eval::n_q_points);
+ VectorizedArray<Number> temp_data[(temp_size > 0 && temp_size < 100) ? 2*temp_size : 1];
+ VectorizedArray<Number> *temp1;
+ VectorizedArray<Number> *temp2;
+ if (temp_size == 0)
+ {
+ temp1 = scratch_data;
+ temp2 = temp1 + std::max(Utilities::fixed_power<dim>(shape_info.fe_degree+1),
+ Utilities::fixed_power<dim>(shape_info.n_q_points_1d));
+ }
+ else if (temp_size > 100)
+ {
+ temp1 = scratch_data;
+ temp2 = temp1 + temp_size;
+ }
+ else
+ {
+ temp1 = &temp_data[0];
+ temp2 = temp1 + temp_size;
+ }
+
+ VectorizedArray<Number> **values_dofs = values_dofs_actual;
+ VectorizedArray<Number> *expanded_dof_values[n_components];
+ if (type == MatrixFreeFunctions::truncated_tensor)
+ {
+ values_dofs = expanded_dof_values;
+ for (unsigned int c=0; c<n_components; ++c)
+ expanded_dof_values[c] = scratch_data+2*(std::max(shape_info.dofs_per_cell,
+ shape_info.n_q_points)) +
+ c*Utilities::fixed_power<dim>(shape_info.fe_degree+1);
+ const int degree = fe_degree != -1 ? fe_degree : shape_info.fe_degree;
+ unsigned int count_p = 0, count_q = 0;
+ for (int i=0; i<(dim>2?degree+1:1); ++i)
+ {
+ for (int j=0; j<(dim>1?degree+1-i:1); ++j)
+ {
+ for (int k=0; k<degree+1-j-i; ++k, ++count_p, ++count_q)
+ for (unsigned int c=0; c<n_components; ++c)
+ expanded_dof_values[c][count_q] = values_dofs_actual[c][count_p];
+ for (int k=degree+1-j-i; k<degree+1; ++k, ++count_q)
+ for (unsigned int c=0; c<n_components; ++c)
+ expanded_dof_values[c][count_q] = VectorizedArray<Number>();
+ }
+ for (int j=degree+1-i; j<degree+1; ++j)
+ for (int k=0; k<degree+1; ++k, ++count_q)
+ for (unsigned int c=0; c<n_components; ++c)
+ expanded_dof_values[c][count_q] = VectorizedArray<Number>();
+ }
+ AssertDimension(count_q, Utilities::fixed_power<dim>(shape_info.fe_degree+1));
+ }
+
+ // These avoid compiler errors; they are only used in sensible context but
+ // compilers typically cannot detect when we access something like
+ // gradients_quad[2] only for dim==3.
+ const unsigned int d1 = dim>1?1:0;
+ const unsigned int d2 = dim>2?2:0;
+ const unsigned int d3 = dim>2?3:0;
+ const unsigned int d4 = dim>2?4:0;
+ const unsigned int d5 = dim>2?5:0;
+
+ switch (dim)
+ {
+ case 1:
+ for (unsigned int c=0; c<n_components; c++)
+ {
+ if (evaluate_val == true)
+ eval.template values<0,true,false> (values_dofs[c], values_quad[c]);
+ if (evaluate_grad == true)
+ eval.template gradients<0,true,false>(values_dofs[c], gradients_quad[c][0]);
+ if (evaluate_lapl == true)
+ eval.template hessians<0,true,false> (values_dofs[c], hessians_quad[c][0]);
+ }
+ break;
+
+ case 2:
+ for (unsigned int c=0; c<n_components; c++)
+ {
+ // grad x
+ if (evaluate_grad == true)
+ {
+ eval.template gradients<0,true,false> (values_dofs[c], temp1);
+ eval.template values<1,true,false> (temp1, gradients_quad[c][0]);
+ }
+ if (evaluate_lapl == true)
+ {
+ // grad xy
+ if (evaluate_grad == false)
+ eval.template gradients<0,true,false>(values_dofs[c], temp1);
+ eval.template gradients<1,true,false> (temp1, hessians_quad[c][d1+d1]);
+
+ // grad xx
+ eval.template hessians<0,true,false>(values_dofs[c], temp1);
+ eval.template values<1,true,false> (temp1, hessians_quad[c][0]);
+ }
+
+ // grad y
+ eval.template values<0,true,false> (values_dofs[c], temp1);
+ if (evaluate_grad == true)
+ eval.template gradients<1,true,false> (temp1, gradients_quad[c][d1]);
+
+ // grad yy
+ if (evaluate_lapl == true)
+ eval.template hessians<1,true,false> (temp1, hessians_quad[c][d1]);
+
+ // val: can use values applied in x
+ if (evaluate_val == true)
+ eval.template values<1,true,false> (temp1, values_quad[c]);
+ }
+ break;
+
+ case 3:
+ for (unsigned int c=0; c<n_components; c++)
+ {
+ if (evaluate_grad == true)
+ {
+ // grad x
+ eval.template gradients<0,true,false> (values_dofs[c], temp1);
+ eval.template values<1,true,false> (temp1, temp2);
+ eval.template values<2,true,false> (temp2, gradients_quad[c][0]);
+ }
+
+ if (evaluate_lapl == true)
+ {
+ // grad xz
+ if (evaluate_grad == false)
+ {
+ eval.template gradients<0,true,false> (values_dofs[c], temp1);
+ eval.template values<1,true,false> (temp1, temp2);
+ }
+ eval.template gradients<2,true,false> (temp2, hessians_quad[c][d4]);
+
+ // grad xy
+ eval.template gradients<1,true,false> (temp1, temp2);
+ eval.template values<2,true,false> (temp2, hessians_quad[c][d3]);
+
+ // grad xx
+ eval.template hessians<0,true,false>(values_dofs[c], temp1);
+ eval.template values<1,true,false> (temp1, temp2);
+ eval.template values<2,true,false> (temp2, hessians_quad[c][0]);
+ }
+
+ // grad y
+ eval.template values<0,true,false> (values_dofs[c], temp1);
+ if (evaluate_grad == true)
+ {
+ eval.template gradients<1,true,false>(temp1, temp2);
+ eval.template values<2,true,false> (temp2, gradients_quad[c][d1]);
+ }
+
+ if (evaluate_lapl == true)
+ {
+ // grad yz
+ if (evaluate_grad == false)
+ eval.template gradients<1,true,false>(temp1, temp2);
+ eval.template gradients<2,true,false> (temp2, hessians_quad[c][d5]);
+
+ // grad yy
+ eval.template hessians<1,true,false> (temp1, temp2);
+ eval.template values<2,true,false> (temp2, hessians_quad[c][d1]);
+ }
+
+ // grad z: can use the values applied in x direction stored in temp1
+ eval.template values<1,true,false> (temp1, temp2);
+ if (evaluate_grad == true)
+ eval.template gradients<2,true,false> (temp2, gradients_quad[c][d2]);
+
+ // grad zz: can use the values applied in x and y direction stored
+ // in temp2
+ if (evaluate_lapl == true)
+ eval.template hessians<2,true,false>(temp2, hessians_quad[c][d2]);
+
+ // val: can use the values applied in x & y direction stored in temp2
+ if (evaluate_val == true)
+ eval.template values<2,true,false> (temp2, values_quad[c]);
+ }
+ break;
+
+ default:
+ AssertThrow(false, ExcNotImplemented());
+ }
+
+ // case additional dof for FE_Q_DG0: add values; gradients and second
+ // derivatives evaluate to zero
+ if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0 && evaluate_val)
+ for (unsigned int c=0; c<n_components; ++c)
+ for (unsigned int q=0; q<shape_info.n_q_points; ++q)
+ values_quad[c][q] += values_dofs[c][shape_info.dofs_per_cell-1];
+ }
+
+
+
+ template <MatrixFreeFunctions::ElementType type, int dim, int fe_degree,
+ int n_q_points_1d, int n_components, typename Number>
+ inline
+ void
+ FEEvaluationImpl<type,dim,fe_degree,n_q_points_1d,n_components,Number>
+ ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ VectorizedArray<Number> *values_dofs_actual[],
+ VectorizedArray<Number> *values_quad[],
+ VectorizedArray<Number> *gradients_quad[][dim],
+ VectorizedArray<Number> *scratch_data,
+ const bool integrate_val,
+ const bool integrate_grad)
+ {
+ const EvaluatorVariant variant =
+ EvaluatorSelector<type,(fe_degree+n_q_points_1d>4)>::variant;
+ typedef EvaluatorTensorProduct<variant, dim, fe_degree, n_q_points_1d,
+ VectorizedArray<Number> > Eval;
+ Eval eval (variant == evaluate_evenodd ? shape_info.shape_val_evenodd :
+ shape_info.shape_values,
+ variant == evaluate_evenodd ? shape_info.shape_gra_evenodd :
+ shape_info.shape_gradients,
+ variant == evaluate_evenodd ? shape_info.shape_hes_evenodd :
+ shape_info.shape_hessians,
+ shape_info.fe_degree,
+ shape_info.n_q_points_1d);
+
+ const unsigned int temp_size = Eval::dofs_per_cell == numbers::invalid_unsigned_int ? 0
+ : (Eval::dofs_per_cell > Eval::n_q_points ?
+ Eval::dofs_per_cell : Eval::n_q_points);
+ VectorizedArray<Number> temp_data[(temp_size > 0 && temp_size < 100) ? 2*temp_size : 1];
+ VectorizedArray<Number> *temp1;
+ VectorizedArray<Number> *temp2;
+ if (temp_size == 0)
+ {
+ temp1 = scratch_data;
+ temp2 = temp1 + std::max(Utilities::fixed_power<dim>(shape_info.fe_degree+1),
+ Utilities::fixed_power<dim>(shape_info.n_q_points_1d));
+ }
+ else if (temp_size > 100)
+ {
+ temp1 = scratch_data;
+ temp2 = temp1 + temp_size;
+ }
+ else
+ {
+ temp1 = &temp_data[0];
+ temp2 = temp1 + temp_size;
+ }
+
+ // expand dof_values to tensor product for truncated tensor products
+ VectorizedArray<Number> **values_dofs = values_dofs_actual;
+ VectorizedArray<Number> *expanded_dof_values[n_components];
+ if (type == MatrixFreeFunctions::truncated_tensor)
+ {
+ values_dofs = expanded_dof_values;
+ for (unsigned int c=0; c<n_components; ++c)
+ expanded_dof_values[c] = scratch_data+2*(std::max(shape_info.dofs_per_cell,
+ shape_info.n_q_points)) +
+ c*Utilities::fixed_power<dim>(shape_info.fe_degree+1);
+ }
+
+ // These avoid compiler errors; they are only used in sensible context but
+ // compilers typically cannot detect when we access something like
+ // gradients_quad[2] only for dim==3.
+ const unsigned int d1 = dim>1?1:0;
+ const unsigned int d2 = dim>2?2:0;
+
+ switch (dim)
+ {
+ case 1:
+ for (unsigned int c=0; c<n_components; c++)
+ {
+ if (integrate_val == true)
+ eval.template values<0,false,false> (values_quad[c], values_dofs[c]);
+ if (integrate_grad == true)
+ {
+ if (integrate_val == true)
+ eval.template gradients<0,false,true> (gradients_quad[c][0], values_dofs[c]);
+ else
+ eval.template gradients<0,false,false> (gradients_quad[c][0], values_dofs[c]);
+ }
+ }
+ break;
+
+ case 2:
+ for (unsigned int c=0; c<n_components; c++)
+ {
+ if (integrate_val == true)
+ {
+ // val
+ eval.template values<0,false,false> (values_quad[c], temp1);
+ //grad x
+ if (integrate_grad == true)
+ eval.template gradients<0,false,true> (gradients_quad[c][0], temp1);
+ eval.template values<1,false,false>(temp1, values_dofs[c]);
+ }
+ if (integrate_grad == true)
+ {
+ // grad y
+ eval.template values<0,false,false> (gradients_quad[c][d1], temp1);
+ if (integrate_val == false)
+ {
+ eval.template gradients<1,false,false>(temp1, values_dofs[c]);
+ //grad x
+ eval.template gradients<0,false,false> (gradients_quad[c][0], temp1);
+ eval.template values<1,false,true> (temp1, values_dofs[c]);
+ }
+ else
+ eval.template gradients<1,false,true>(temp1, values_dofs[c]);
+ }
+ }
+ break;
+
+ case 3:
+ for (unsigned int c=0; c<n_components; c++)
+ {
+ if (integrate_val == true)
+ {
+ // val
+ eval.template values<0,false,false> (values_quad[c], temp1);
+ //grad x: can sum to temporary value in temp1
+ if (integrate_grad == true)
+ eval.template gradients<0,false,true> (gradients_quad[c][0], temp1);
+ eval.template values<1,false,false>(temp1, temp2);
+ if (integrate_grad == true)
+ {
+ eval.template values<0,false,false> (gradients_quad[c][d1], temp1);
+ eval.template gradients<1,false,true>(temp1, temp2);
+ }
+ eval.template values<2,false,false> (temp2, values_dofs[c]);
+ }
+ else if (integrate_grad == true)
+ {
+ eval.template gradients<0,false,false>(gradients_quad[c][0], temp1);
+ eval.template values<1,false,false> (temp1, temp2);
+ eval.template values<0,false,false> (gradients_quad[c][d1], temp1);
+ eval.template gradients<1,false,true>(temp1, temp2);
+ eval.template values<2,false,false> (temp2, values_dofs[c]);
+ }
+ if (integrate_grad == true)
+ {
+ // grad z: can sum to temporary x and y value in output
+ eval.template values<0,false,false> (gradients_quad[c][d2], temp1);
+ eval.template values<1,false,false> (temp1, temp2);
+ eval.template gradients<2,false,true> (temp2, values_dofs[c]);
+ }
+ }
+ break;
+
+ default:
+ AssertThrow(false, ExcNotImplemented());
+ }
+
+ // case FE_Q_DG0: add values, gradients and second derivatives are zero
+ if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0)
+ {
+ if (integrate_val)
+ for (unsigned int c=0; c<n_components; ++c)
+ {
+ values_dofs[c][shape_info.dofs_per_cell-1] = values_quad[c][0];
+ for (unsigned int q=1; q<shape_info.n_q_points; ++q)
+ values_dofs[c][shape_info.dofs_per_cell-1] += values_quad[c][q];
+ }
+ else
+ for (unsigned int c=0; c<n_components; ++c)
+ values_dofs[c][shape_info.dofs_per_cell-1] = VectorizedArray<Number>();
+ }
+
+ if (type == MatrixFreeFunctions::truncated_tensor)
+ {
+ unsigned int count_p = 0, count_q = 0;
+ const int degree = fe_degree != -1 ? fe_degree : shape_info.fe_degree;
+ for (int i=0; i<(dim>2?degree+1:1); ++i)
+ {
+ for (int j=0; j<(dim>1?degree+1-i:1); ++j)
+ {
+ for (int k=0; k<degree+1-j-i; ++k, ++count_p, ++count_q)
+ {
+ for (unsigned int c=0; c<n_components; ++c)
+ values_dofs_actual[c][count_p] = expanded_dof_values[c][count_q];
+ }
+ count_q += j+i;
+ }
+ count_q += i*(degree+1);
+ }
+ AssertDimension(count_q, Utilities::fixed_power<dim>(shape_info.fe_degree+1));
+ }
+ }
+
+ // This a specialization for Gauss-Lobatto elements where the 'values'
+ // operation is identity, which allows us to write shorter code.
+ template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
+ struct FEEvaluationImpl<MatrixFreeFunctions::tensor_gausslobatto, dim,
+ fe_degree, n_q_points_1d, n_components, Number>
+ {
+ static
+ void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ VectorizedArray<Number> *values_dofs[],
+ VectorizedArray<Number> *values_quad[],
+ VectorizedArray<Number> *gradients_quad[][dim],
+ VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
+ VectorizedArray<Number> *scratch_data,
+ const bool evaluate_val,
+ const bool evaluate_grad,
+ const bool evaluate_lapl);
+
+ static
+ void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ VectorizedArray<Number> *values_dofs[],
+ VectorizedArray<Number> *values_quad[],
+ VectorizedArray<Number> *gradients_quad[][dim],
+ VectorizedArray<Number> *scratch_data,
+ const bool integrate_val,
+ const bool integrate_grad);
+ };
+
+ template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
+ inline
+ void
+ FEEvaluationImpl<MatrixFreeFunctions::tensor_gausslobatto, dim,
+ fe_degree, n_q_points_1d, n_components, Number>
+ ::evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ VectorizedArray<Number> *values_dofs[],
+ VectorizedArray<Number> *values_quad[],
+ VectorizedArray<Number> *gradients_quad[][dim],
+ VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
+ VectorizedArray<Number> *scratch_data,
+ const bool evaluate_val,
+ const bool evaluate_grad,
+ const bool evaluate_lapl)
+ {
+ typedef EvaluatorTensorProduct<evaluate_evenodd, dim, fe_degree, fe_degree+1,
+ VectorizedArray<Number> > Eval;
+ Eval eval (shape_info.shape_val_evenodd,
+ shape_info.shape_gra_evenodd,
+ shape_info.shape_hes_evenodd,
+ shape_info.fe_degree,
+ shape_info.n_q_points_1d);
+
+ // These avoid compiler errors; they are only used in sensible context but
+ // compilers typically cannot detect when we access something like
+ // gradients_quad[2] only for dim==3.
+ const unsigned int d1 = dim>1?1:0;
+ const unsigned int d2 = dim>2?2:0;
+ const unsigned int d3 = dim>2?3:0;
+ const unsigned int d4 = dim>2?4:0;
+ const unsigned int d5 = dim>2?5:0;
+
+ switch (dim)
+ {
+ case 1:
+ if (evaluate_val == true)
+ std::memcpy (values_quad[0], values_dofs[0],
+ eval.dofs_per_cell * n_components *
+ sizeof (values_dofs[0][0]));
+ for (unsigned int c=0; c<n_components; c++)
+ {
+ if (evaluate_grad == true)
+ eval.template gradients<0,true,false>(values_dofs[c], gradients_quad[c][0]);
+ if (evaluate_lapl == true)
+ eval.template hessians<0,true,false> (values_dofs[c], hessians_quad[c][0]);
+ }
+ break;
+
+ case 2:
+ if (evaluate_val == true)
+ {
+ std::memcpy (values_quad[0], values_dofs[0],
+ Eval::dofs_per_cell * n_components *
+ sizeof (values_dofs[0][0]));
+ }
+ if (evaluate_grad == true)
+ for (unsigned int comp=0; comp<n_components; comp++)
+ {
+ // grad x
+ eval.template gradients<0,true,false> (values_dofs[comp],
+ gradients_quad[comp][0]);
+ // grad y
+ eval.template gradients<1,true,false> (values_dofs[comp],
+ gradients_quad[comp][d1]);
+ }
+ if (evaluate_lapl == true)
+ for (unsigned int comp=0; comp<n_components; comp++)
+ {
+ // hess x
+ eval.template hessians<0,true,false> (values_dofs[comp],
+ hessians_quad[comp][0]);
+ // hess y
+ eval.template hessians<1,true,false> (values_dofs[comp],
+ hessians_quad[comp][d1]);
+
+ // grad x grad y
+ eval.template gradients<0,true,false> (values_dofs[comp], scratch_data);
+ eval.template gradients<1,true,false> (scratch_data, hessians_quad[comp][d1+d1]);
+ }
+ break;
+
+ case 3:
+ if (evaluate_val == true)
+ {
+ std::memcpy (values_quad[0], values_dofs[0],
+ Eval::dofs_per_cell * n_components *
+ sizeof (values_dofs[0][0]));
+ }
+ if (evaluate_grad == true)
+ for (unsigned int comp=0; comp<n_components; comp++)
+ {
+ // grad x
+ eval.template gradients<0,true,false> (values_dofs[comp],
+ gradients_quad[comp][0]);
+ // grad y
+ eval.template gradients<1,true,false> (values_dofs[comp],
+ gradients_quad[comp][d1]);
+ // grad y
+ eval.template gradients<2,true,false> (values_dofs[comp],
+ gradients_quad[comp][d2]);
+ }
+ if (evaluate_lapl == true)
+ for (unsigned int comp=0; comp<n_components; comp++)
+ {
+ // grad x
+ eval.template hessians<0,true,false> (values_dofs[comp],
+ hessians_quad[comp][0]);
+ // grad y
+ eval.template hessians<1,true,false> (values_dofs[comp],
+ hessians_quad[comp][d1]);
+ // grad y
+ eval.template hessians<2,true,false> (values_dofs[comp],
+ hessians_quad[comp][d2]);
+
+ VectorizedArray<Number> *temp1 = scratch_data;
+ // grad xy
+ eval.template gradients<0,true,false> (values_dofs[comp], temp1);
+ eval.template gradients<1,true,false> (temp1, hessians_quad[comp][d3]);
+ // grad xz
+ eval.template gradients<2,true,false> (temp1, hessians_quad[comp][d4]);
+ // grad yz
+ eval.template gradients<1,true,false> (values_dofs[comp], temp1);
+ eval.template gradients<2,true,false> (temp1, hessians_quad[comp][d5]);
+ }
+ break;
+ default:
+ AssertThrow(false, ExcNotImplemented());
+ }
+ }
+
+ template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
+ inline
+ void
+ FEEvaluationImpl<MatrixFreeFunctions::tensor_gausslobatto, dim,
+ fe_degree, n_q_points_1d, n_components, Number>
+ ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ VectorizedArray<Number> *values_dofs[],
+ VectorizedArray<Number> *values_quad[],
+ VectorizedArray<Number> *gradients_quad[][dim],
+ VectorizedArray<Number> *,
+ const bool integrate_val,
+ const bool integrate_grad)
+ {
+ typedef EvaluatorTensorProduct<evaluate_evenodd, dim, fe_degree, fe_degree+1,
+ VectorizedArray<Number> > Eval;
+ Eval eval (shape_info.shape_val_evenodd,
+ shape_info.shape_gra_evenodd,
+ shape_info.shape_hes_evenodd,
+ shape_info.fe_degree,
+ shape_info.n_q_points_1d);
+
+ // These avoid compiler errors; they are only used in sensible context but
+ // compilers typically cannot detect when we access something like
+ // gradients_quad[2] only for dim==3.
+ const unsigned int d1 = dim>1?1:0;
+ const unsigned int d2 = dim>2?2:0;
+
+ if (integrate_val == true)
+ std::memcpy (values_dofs[0], values_quad[0],
+ Eval::dofs_per_cell * n_components *
+ sizeof (values_dofs[0][0]));
+ switch (dim)
+ {
+ case 1:
+ for (unsigned int c=0; c<n_components; c++)
+ {
+ if (integrate_grad == true)
+ {
+ if (integrate_val == true)
+ eval.template gradients<0,false,true> (gradients_quad[c][0],
+ values_dofs[c]);
+ else
+ eval.template gradients<0,false,false> (gradients_quad[c][0],
+ values_dofs[c]);
+ }
+ }
+
+ break;
+ case 2:
+ if (integrate_grad == true)
+ for (unsigned int comp=0; comp<n_components; comp++)
+ {
+ // grad x: If integrate_val == true we have to add to the
+ // previous output
+ if (integrate_val == true)
+ eval.template gradients<0, false, true> (gradients_quad[comp][0],
+ values_dofs[comp]);
+ else
+ eval.template gradients<0, false, false> (gradients_quad[comp][0],
+ values_dofs[comp]);
+
+ // grad y
+ eval.template gradients<1, false, true> (gradients_quad[comp][d1],
+ values_dofs[comp]);
+ }
+ break;
+
+ case 3:
+ if (integrate_grad == true)
+ for (unsigned int comp=0; comp<n_components; comp++)
+ {
+ // grad x: If integrate_val == true we have to add to the
+ // previous output
+ if (integrate_val == true)
+ eval.template gradients<0, false, true> (gradients_quad[comp][0],
+ values_dofs[comp]);
+ else
+ eval.template gradients<0, false, false> (gradients_quad[comp][0],
+ values_dofs[comp]);
+
+ // grad y
+ eval.template gradients<1, false, true> (gradients_quad[comp][d1],
+ values_dofs[comp]);
+
+ // grad z
+ eval.template gradients<2, false, true> (gradients_quad[comp][d2],
+ values_dofs[comp]);
+ }
+ break;
+
+ default:
+ AssertThrow(false, ExcNotImplemented());
+ }
+ }
+
+} // end of namespace internal
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
#include <deal.II/matrix_free/mapping_data_on_the_fly.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/matrix_free/shape_info.h>
+#include <deal.II/matrix_free/evaluation_kernels.h>
+#include <deal.II/matrix_free/tensor_product_kernels.h>
DEAL_II_NAMESPACE_OPEN
-namespace internal
-{
- /**
- * In this namespace, the evaluator routines that evaluate the tensor
- * products are implemented.
- */
- enum EvaluatorVariant
- {
- /**
- * Do not use anything more than the tensor product structure of the
- * finite element.
- */
- evaluate_general,
- /**
- * Perform evaluation by exploiting symmetry in the finite element: i.e.,
- * skip some computations by utilizing the symmetry in the shape functions
- * and quadrature points.
- */
- evaluate_symmetric,
- /**
- * Use symmetry to apply the operator to even and odd parts of the input
- * vector separately: see the documentation of the EvaluatorTensorProduct
- * specialization for more information.
- */
- evaluate_evenodd
- };
-
- /**
- * Generic evaluator framework
- */
- template <EvaluatorVariant variant, int dim, int fe_degree, int n_q_points_1d,
- typename Number>
- struct EvaluatorTensorProduct
- {};
-
- /**
- * Internal evaluator for 1d-3d shape function using the tensor product form
- * of the basis functions
- */
- template <int dim, int fe_degree, int n_q_points_1d, typename Number>
- struct EvaluatorTensorProduct<evaluate_general,dim,fe_degree,n_q_points_1d,Number>
- {
- static const unsigned int dofs_per_cell = Utilities::fixed_int_power<fe_degree+1,dim>::value;
- static const unsigned int n_q_points = Utilities::fixed_int_power<n_q_points_1d,dim>::value;
-
- /**
- * Empty constructor. Does nothing. Be careful when using 'values' and
- * related methods because they need to be filled with the other pointer
- */
- EvaluatorTensorProduct ()
- :
- shape_values (0),
- shape_gradients (0),
- shape_hessians (0)
- {}
-
- /**
- * Constructor, taking the data from ShapeInfo
- */
- EvaluatorTensorProduct (const AlignedVector<Number> &shape_values,
- const AlignedVector<Number> &shape_gradients,
- const AlignedVector<Number> &shape_hessians,
- const unsigned int dummy1 = 0,
- const unsigned int dummy2 = 0)
- :
- shape_values (shape_values.begin()),
- shape_gradients (shape_gradients.begin()),
- shape_hessians (shape_hessians.begin())
- {
- (void)dummy1;
- (void)dummy2;
- }
-
- template <int direction, bool dof_to_quad, bool add>
- void
- values (const Number in [],
- Number out[]) const
- {
- apply<direction,dof_to_quad,add>(shape_values, in, out);
- }
-
- template <int direction, bool dof_to_quad, bool add>
- void
- gradients (const Number in [],
- Number out[]) const
- {
- apply<direction,dof_to_quad,add>(shape_gradients, in, out);
- }
-
- template <int direction, bool dof_to_quad, bool add>
- void
- hessians (const Number in [],
- Number out[]) const
- {
- apply<direction,dof_to_quad,add>(shape_hessians, in, out);
- }
-
- template <int direction, bool dof_to_quad, bool add>
- static void apply (const Number *shape_data,
- const Number in [],
- Number out []);
-
- const Number *shape_values;
- const Number *shape_gradients;
- const Number *shape_hessians;
- };
-
- // evaluates the given shape data in 1d-3d using the tensor 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>
- template <int direction, bool dof_to_quad, bool add>
- inline
- void
- EvaluatorTensorProduct<evaluate_general,dim,fe_degree,n_q_points_1d,Number>
- ::apply (const Number *shape_data,
- const Number in [],
- 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_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<nn; ++col)
- {
- Number val0;
- if (dof_to_quad == true)
- val0 = shape_data[col];
- else
- val0 = shape_data[col*n_q_points_1d];
- Number res0 = val0 * in[0];
- for (int ind=1; ind<mm; ++ind)
- {
- if (dof_to_quad == true)
- val0 = shape_data[ind*n_q_points_1d+col];
- else
- val0 = shape_data[col*n_q_points_1d+ind];
- res0 += val0 * in[stride*ind];
- }
- if (add == false)
- out[stride*col] = res0;
- else
- out[stride*col] += 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);
- }
- }
- }
-
-
-
- // 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:
- ++in;
- ++out;
- // faces 2 and 3 in 3D use local coordinate system zx, which
- // is the other way around compared to the tensor
- // product. Need to take that into account.
- if (dim == 3)
- {
- if (dof_to_quad)
- out += fe_degree;
- else
- in += fe_degree;
- }
- break;
- case 2:
- ++in;
- ++out;
- break;
- default:
- Assert (false, ExcNotImplemented());
- }
- }
- if (face_direction == 1 && dim == 3)
- {
- in += mm*(mm-1);
- out += nn*(nn-1);
- // adjust for local coordinate system zx
- if (dof_to_quad)
- out -= (fe_degree+1)*(fe_degree+1)-1;
- else
- in -= (fe_degree+1)*(fe_degree+1)-1;
- }
- }
- }
-
-
-
- /**
- * Internal evaluator for 1d-3d shape function using the tensor product form
- * of the basis functions. The same as above but without making use of
- * template arguments and rather variable loop bounds.
- */
- template <int dim, typename Number>
- struct EvaluatorTensorProduct<evaluate_general,dim,-1,0,Number>
- {
- static const unsigned int dofs_per_cell = numbers::invalid_unsigned_int;
- static const unsigned int n_q_points = numbers::invalid_unsigned_int;
-
- /**
- * Empty constructor. Does nothing. Be careful when using 'values' and
- * related methods because they need to be filled with the other constructor
- */
- EvaluatorTensorProduct ()
- :
- shape_values (0),
- shape_gradients (0),
- shape_hessians (0),
- fe_degree (numbers::invalid_unsigned_int),
- n_q_points_1d (numbers::invalid_unsigned_int)
- {}
-
- /**
- * Constructor, taking the data from ShapeInfo
- */
- EvaluatorTensorProduct (const AlignedVector<Number> &shape_values,
- const AlignedVector<Number> &shape_gradients,
- const AlignedVector<Number> &shape_hessians,
- const unsigned int fe_degree,
- const unsigned int n_q_points_1d)
- :
- shape_values (shape_values.begin()),
- shape_gradients (shape_gradients.begin()),
- shape_hessians (shape_hessians.begin()),
- fe_degree (fe_degree),
- n_q_points_1d (n_q_points_1d)
- {}
-
- template <int direction, bool dof_to_quad, bool add>
- void
- values (const Number *in,
- Number *out) const
- {
- apply<direction,dof_to_quad,add>(shape_values, in, out);
- }
-
- template <int direction, bool dof_to_quad, bool add>
- void
- gradients (const Number *in,
- Number *out) const
- {
- apply<direction,dof_to_quad,add>(shape_gradients, in, out);
- }
-
- template <int direction, bool dof_to_quad, bool add>
- void
- hessians (const Number *in,
- Number *out) const
- {
- apply<direction,dof_to_quad,add>(shape_hessians, in, out);
- }
-
- template <int direction, bool dof_to_quad, bool add>
- void apply (const Number *shape_data,
- const Number *in,
- Number *out) const;
-
- const Number *shape_values;
- const Number *shape_gradients;
- const Number *shape_hessians;
- const unsigned int fe_degree;
- const unsigned int n_q_points_1d;
- };
-
- // evaluates the given shape data in 1d-3d using the tensor 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, typename Number>
- template <int direction, bool dof_to_quad, bool add>
- inline
- void
- EvaluatorTensorProduct<evaluate_general,dim,-1,0,Number>
- ::apply (const Number *shape_data,
- const Number *in,
- Number *out) const
- {
- 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_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
- const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
- const int stride = direction==0 ? 1 : Utilities::fixed_power<direction>(nn);
-
- for (int i2=0; i2<n_blocks2; ++i2)
- {
- for (int i1=0; i1<n_blocks1; ++i1)
- {
- for (int col=0; col<nn; ++col)
- {
- Number val0;
- if (dof_to_quad == true)
- val0 = shape_data[col];
- else
- val0 = shape_data[col*n_q_points_1d];
- Number res0 = val0 * in[0];
- for (int ind=1; ind<mm; ++ind)
- {
- if (dof_to_quad == true)
- val0 = shape_data[ind*n_q_points_1d+col];
- else
- val0 = shape_data[col*n_q_points_1d+ind];
- res0 += val0 * in[stride*ind];
- }
- if (add == false)
- out[stride*col] = res0;
- else
- out[stride*col] += 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 evaluator for 1d-3d shape function using the tensor product form
- * of the basis functions. This class 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.
- */
- template <int dim, int fe_degree, int n_q_points_1d, typename Number>
- struct EvaluatorTensorProduct<evaluate_symmetric,dim,fe_degree,n_q_points_1d,Number>
- {
- static const unsigned int dofs_per_cell = Utilities::fixed_int_power<fe_degree+1,dim>::value;
- static const unsigned int n_q_points = Utilities::fixed_int_power<n_q_points_1d,dim>::value;
-
- /**
- * Constructor, taking the data from ShapeInfo
- */
- EvaluatorTensorProduct (const AlignedVector<Number> &shape_values,
- const AlignedVector<Number> &shape_gradients,
- const AlignedVector<Number> &shape_hessians,
- const unsigned int dummy1 = 0,
- const unsigned int dummy2 = 0)
- :
- shape_values (shape_values.begin()),
- shape_gradients (shape_gradients.begin()),
- shape_hessians (shape_hessians.begin())
- {
- (void)dummy1;
- (void)dummy2;
- }
-
- template <int direction, bool dof_to_quad, bool add>
- void
- values (const Number in [],
- Number out[]) const;
-
- template <int direction, bool dof_to_quad, bool add>
- void
- gradients (const Number in [],
- Number out[]) const;
-
- template <int direction, bool dof_to_quad, bool add>
- void
- hessians (const Number in [],
- Number out[]) const;
-
- const Number *shape_values;
- const Number *shape_gradients;
- const Number *shape_hessians;
- };
-
-
-
- // In this 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>
- template <int direction, bool dof_to_quad, bool add>
- inline
- void
- EvaluatorTensorProduct<evaluate_symmetric,dim,fe_degree,n_q_points_1d,Number>
- ::values (const Number in [],
- Number out []) const
- {
- 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)
- {
- 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 = 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)
- {
- Number res0;
- 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];
- Number val1 = in[stride*ind] + in[stride*(mm-1-ind)];
- val1 *= val0;
- res0 += val1;
- }
- }
- else
- res0 = Number();
- if (add == false)
- out[stride*n_cols] = res0;
- else
- out[stride*n_cols] += res0;
- }
- else if (dof_to_quad == false && nn%2 == 1)
- {
- Number res0;
- if (mid > 0)
- {
- 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];
- 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);
- }
- }
- }
-
-
-
- // 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>
- template <int direction, bool dof_to_quad, bool add>
- inline
- void
- EvaluatorTensorProduct<evaluate_symmetric,dim,fe_degree,n_q_points_1d,Number>
- ::gradients (const Number in [],
- Number out []) const
- {
- 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)
- {
- 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 = 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 )
- {
- 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];
- 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>
- template <int direction, bool dof_to_quad, bool add>
- inline
- void
- EvaluatorTensorProduct<evaluate_symmetric,dim,fe_degree,n_q_points_1d,Number>
- ::hessians (const Number in [],
- Number out []) const
- {
- 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)
- {
- 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 = 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 )
- {
- 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];
- Number val1 = in[stride*ind] + in[stride*(mm-1-ind)];
- val1 *= val0;
- res0 += val1;
- }
- }
- else
- res0 = 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);
- }
- }
- }
-
-
-
- /**
- * Internal evaluator for 1d-3d shape function using the tensor product form
- * of the basis functions.
- *
- * This class 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>
- struct EvaluatorTensorProduct<evaluate_evenodd,dim,fe_degree,n_q_points_1d,Number>
- {
- static const unsigned int dofs_per_cell = Utilities::fixed_int_power<fe_degree+1,dim>::value;
- static const unsigned int n_q_points = Utilities::fixed_int_power<n_q_points_1d,dim>::value;
-
- /**
- * Empty constructor. Does nothing. Be careful when using 'values' and
- * related methods because they need to be filled with the other pointer
- */
- EvaluatorTensorProduct ()
- :
- shape_values (0),
- shape_gradients (0),
- shape_hessians (0)
- {}
-
- /**
- * Constructor, taking the data from ShapeInfo (using the even-odd
- * variants stored there)
- */
- EvaluatorTensorProduct (const AlignedVector<Number> &shape_values,
- const AlignedVector<Number> &shape_gradients,
- const AlignedVector<Number> &shape_hessians,
- const unsigned int dummy1 = 0,
- const unsigned int dummy2 = 0)
- :
- shape_values (shape_values.begin()),
- shape_gradients (shape_gradients.begin()),
- shape_hessians (shape_hessians.begin())
- {
- (void)dummy1;
- (void)dummy2;
- }
-
- template <int direction, bool dof_to_quad, bool add>
- void
- values (const Number in [],
- Number out[]) const
- {
- apply<direction,dof_to_quad,add,0>(shape_values, in, out);
- }
-
- template <int direction, bool dof_to_quad, bool add>
- void
- gradients (const Number in [],
- Number out[]) const
- {
- apply<direction,dof_to_quad,add,1>(shape_gradients, in, out);
- }
-
- template <int direction, bool dof_to_quad, bool add>
- void
- hessians (const Number in [],
- Number out[]) const
- {
- apply<direction,dof_to_quad,add,2>(shape_hessians, in, out);
- }
-
- template <int direction, bool dof_to_quad, bool add, int type>
- static void apply (const Number *shape_data,
- const Number in [],
- Number out []);
-
- const Number *shape_values;
- const Number *shape_gradients;
- const Number *shape_hessians;
- };
-
-
-
- template <int dim, int fe_degree, int n_q_points_1d, typename Number>
- template <int direction, bool dof_to_quad, bool add, int type>
- inline
- void
- EvaluatorTensorProduct<evaluate_evenodd,dim,fe_degree,n_q_points_1d,Number>
- ::apply (const Number *shapes,
- 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;
-
- const int offset = (n_q_points_1d+1)/2;
-
- // 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>0?mid:1], xm[mid>0?mid:1];
- 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[col] * xp[0];
- r1 = shapes[fe_degree*offset + col] * xm[0];
- }
- else
- {
- r0 = shapes[col*offset] * xp[0];
- r1 = shapes[(fe_degree-col)*offset] * xm[0];
- }
- for (int ind=1; ind<mid; ++ind)
- {
- if (dof_to_quad == true)
- {
- r0 += shapes[ind*offset+col] * xp[ind];
- r1 += shapes[(fe_degree-ind)*offset+col] * xm[ind];
- }
- else
- {
- r0 += shapes[col*offset+ind] * xp[ind];
- r1 += shapes[(fe_degree-col)*offset+ind] * xm[ind];
- }
- }
- }
- else
- r0 = r1 = Number();
- if (mm % 2 == 1 && dof_to_quad == true)
- {
- if (type == 1)
- r1 += shapes[mid*offset+col] * in[stride*mid];
- else
- r0 += shapes[mid*offset+col] * in[stride*mid];
- }
- else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0))
- r0 += shapes[col*offset+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[n_cols] * xp[0];
- for (int ind=1; ind<mid; ++ind)
- r0 += shapes[ind*offset+n_cols] * xp[ind];
- }
- else
- r0 = Number();
- if (type != 1 && mm % 2 == 1)
- r0 += shapes[mid*offset+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*offset] * xm[0];
- for (int ind=1; ind<mid; ++ind)
- r0 += shapes[n_cols*offset+ind] * xm[ind];
- }
- else
- {
- r0 = shapes[n_cols*offset] * xp[0];
- for (int ind=1; ind<mid; ++ind)
- r0 += shapes[n_cols*offset+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*offset+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);
- }
- }
- }
-
-
-
- // Select evaluator type from element shape function type
- template <MatrixFreeFunctions::ElementType element, bool is_long>
- struct EvaluatorSelector {};
-
- template <bool is_long>
- struct EvaluatorSelector<MatrixFreeFunctions::tensor_general,is_long>
- {
- static const EvaluatorVariant variant = evaluate_general;
- };
-
- template <>
- struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric,false>
- {
- static const EvaluatorVariant variant = evaluate_symmetric;
- };
-
- template <> struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric,true>
- {
- static const EvaluatorVariant variant = evaluate_evenodd;
- };
-
- template <bool is_long>
- struct EvaluatorSelector<MatrixFreeFunctions::truncated_tensor,is_long>
- {
- static const EvaluatorVariant variant = evaluate_general;
- };
-
- template <>
- struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0,false>
- {
- static const EvaluatorVariant variant = evaluate_general;
- };
-
- template <>
- struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0,true>
- {
- static const EvaluatorVariant variant = evaluate_evenodd;
- };
-
- template <bool is_long>
- struct EvaluatorSelector<MatrixFreeFunctions::tensor_gausslobatto,is_long>
- {
- static const EvaluatorVariant variant = evaluate_evenodd;
- };
-
-
-
- // This struct performs the evaluation of function values, gradients and
- // Hessians for tensor-product finite elements. The operation is used for
- // both the symmetric and non-symmetric case, which use different apply
- // functions 'values', 'gradients' in the individual coordinate
- // directions. The apply functions for values are provided through one of
- // the template classes EvaluatorTensorProduct which in turn are selected
- // from the MatrixFreeFunctions::ElementType template argument.
- //
- // There is a specialization made for Gauss-Lobatto elements further down
- // where the 'values' operation is identity, which allows us to write
- // shorter code.
- template <MatrixFreeFunctions::ElementType type, int dim, int fe_degree,
- int n_q_points_1d, int n_components, typename Number>
- struct FEEvaluationImpl
- {
- static
- void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
- VectorizedArray<Number> *scratch_data,
- const bool evaluate_val,
- const bool evaluate_grad,
- const bool evaluate_lapl);
-
- static
- void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *scratch_data,
- const bool evaluate_val,
- const bool evaluate_grad);
- };
-
-
- template <MatrixFreeFunctions::ElementType type, int dim, int fe_degree,
- int n_q_points_1d, int n_components, typename Number>
- inline
- void
- FEEvaluationImpl<type,dim,fe_degree,n_q_points_1d,n_components,Number>
- ::evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
- VectorizedArray<Number> *scratch_data,
- const bool evaluate_val,
- const bool evaluate_grad,
- const bool evaluate_lapl)
- {
- if (evaluate_val == false && evaluate_grad == false && evaluate_lapl == false)
- return;
-
- const EvaluatorVariant variant =
- EvaluatorSelector<type,(fe_degree+n_q_points_1d>4)>::variant;
- typedef EvaluatorTensorProduct<variant, dim, fe_degree, n_q_points_1d,
- VectorizedArray<Number> > Eval;
- Eval eval (variant == evaluate_evenodd ? shape_info.shape_val_evenodd :
- shape_info.shape_values,
- variant == evaluate_evenodd ? shape_info.shape_gra_evenodd :
- shape_info.shape_gradients,
- variant == evaluate_evenodd ? shape_info.shape_hes_evenodd :
- shape_info.shape_hessians,
- shape_info.fe_degree,
- shape_info.n_q_points_1d);
-
- const unsigned int temp_size = Eval::dofs_per_cell == numbers::invalid_unsigned_int ? 0
- : (Eval::dofs_per_cell > Eval::n_q_points ?
- Eval::dofs_per_cell : Eval::n_q_points);
- VectorizedArray<Number> temp_data[(temp_size > 0 && temp_size < 100) ? 2*temp_size : 1];
- VectorizedArray<Number> *temp1;
- VectorizedArray<Number> *temp2;
- if (temp_size == 0)
- {
- temp1 = scratch_data;
- temp2 = temp1 + std::max(Utilities::fixed_power<dim>(shape_info.fe_degree+1),
- Utilities::fixed_power<dim>(shape_info.n_q_points_1d));
- }
- else if (temp_size > 100)
- {
- temp1 = scratch_data;
- temp2 = temp1 + temp_size;
- }
- else
- {
- temp1 = &temp_data[0];
- temp2 = temp1 + temp_size;
- }
-
- VectorizedArray<Number> **values_dofs = values_dofs_actual;
- VectorizedArray<Number> *expanded_dof_values[n_components];
- if (type == MatrixFreeFunctions::truncated_tensor)
- {
- values_dofs = expanded_dof_values;
- for (unsigned int c=0; c<n_components; ++c)
- expanded_dof_values[c] = scratch_data+2*(std::max(shape_info.dofs_per_cell,
- shape_info.n_q_points)) +
- c*Utilities::fixed_power<dim>(shape_info.fe_degree+1);
- const int degree = fe_degree != -1 ? fe_degree : shape_info.fe_degree;
- unsigned int count_p = 0, count_q = 0;
- for (int i=0; i<(dim>2?degree+1:1); ++i)
- {
- for (int j=0; j<(dim>1?degree+1-i:1); ++j)
- {
- for (int k=0; k<degree+1-j-i; ++k, ++count_p, ++count_q)
- for (unsigned int c=0; c<n_components; ++c)
- expanded_dof_values[c][count_q] = values_dofs_actual[c][count_p];
- for (int k=degree+1-j-i; k<degree+1; ++k, ++count_q)
- for (unsigned int c=0; c<n_components; ++c)
- expanded_dof_values[c][count_q] = VectorizedArray<Number>();
- }
- for (int j=degree+1-i; j<degree+1; ++j)
- for (int k=0; k<degree+1; ++k, ++count_q)
- for (unsigned int c=0; c<n_components; ++c)
- expanded_dof_values[c][count_q] = VectorizedArray<Number>();
- }
- AssertDimension(count_q, Utilities::fixed_power<dim>(shape_info.fe_degree+1));
- }
-
- // These avoid compiler errors; they are only used in sensible context but
- // compilers typically cannot detect when we access something like
- // gradients_quad[2] only for dim==3.
- const unsigned int d1 = dim>1?1:0;
- const unsigned int d2 = dim>2?2:0;
- const unsigned int d3 = dim>2?3:0;
- const unsigned int d4 = dim>2?4:0;
- const unsigned int d5 = dim>2?5:0;
-
- switch (dim)
- {
- case 1:
- for (unsigned int c=0; c<n_components; c++)
- {
- if (evaluate_val == true)
- eval.template values<0,true,false> (values_dofs[c], values_quad[c]);
- if (evaluate_grad == true)
- eval.template gradients<0,true,false>(values_dofs[c], gradients_quad[c][0]);
- if (evaluate_lapl == true)
- eval.template hessians<0,true,false> (values_dofs[c], hessians_quad[c][0]);
- }
- break;
-
- case 2:
- for (unsigned int c=0; c<n_components; c++)
- {
- // grad x
- if (evaluate_grad == true)
- {
- eval.template gradients<0,true,false> (values_dofs[c], temp1);
- eval.template values<1,true,false> (temp1, gradients_quad[c][0]);
- }
- if (evaluate_lapl == true)
- {
- // grad xy
- if (evaluate_grad == false)
- eval.template gradients<0,true,false>(values_dofs[c], temp1);
- eval.template gradients<1,true,false> (temp1, hessians_quad[c][d1+d1]);
-
- // grad xx
- eval.template hessians<0,true,false>(values_dofs[c], temp1);
- eval.template values<1,true,false> (temp1, hessians_quad[c][0]);
- }
-
- // grad y
- eval.template values<0,true,false> (values_dofs[c], temp1);
- if (evaluate_grad == true)
- eval.template gradients<1,true,false> (temp1, gradients_quad[c][d1]);
-
- // grad yy
- if (evaluate_lapl == true)
- eval.template hessians<1,true,false> (temp1, hessians_quad[c][d1]);
-
- // val: can use values applied in x
- if (evaluate_val == true)
- eval.template values<1,true,false> (temp1, values_quad[c]);
- }
- break;
-
- case 3:
- for (unsigned int c=0; c<n_components; c++)
- {
- if (evaluate_grad == true)
- {
- // grad x
- eval.template gradients<0,true,false> (values_dofs[c], temp1);
- eval.template values<1,true,false> (temp1, temp2);
- eval.template values<2,true,false> (temp2, gradients_quad[c][0]);
- }
-
- if (evaluate_lapl == true)
- {
- // grad xz
- if (evaluate_grad == false)
- {
- eval.template gradients<0,true,false> (values_dofs[c], temp1);
- eval.template values<1,true,false> (temp1, temp2);
- }
- eval.template gradients<2,true,false> (temp2, hessians_quad[c][d4]);
-
- // grad xy
- eval.template gradients<1,true,false> (temp1, temp2);
- eval.template values<2,true,false> (temp2, hessians_quad[c][d3]);
-
- // grad xx
- eval.template hessians<0,true,false>(values_dofs[c], temp1);
- eval.template values<1,true,false> (temp1, temp2);
- eval.template values<2,true,false> (temp2, hessians_quad[c][0]);
- }
-
- // grad y
- eval.template values<0,true,false> (values_dofs[c], temp1);
- if (evaluate_grad == true)
- {
- eval.template gradients<1,true,false>(temp1, temp2);
- eval.template values<2,true,false> (temp2, gradients_quad[c][d1]);
- }
-
- if (evaluate_lapl == true)
- {
- // grad yz
- if (evaluate_grad == false)
- eval.template gradients<1,true,false>(temp1, temp2);
- eval.template gradients<2,true,false> (temp2, hessians_quad[c][d5]);
-
- // grad yy
- eval.template hessians<1,true,false> (temp1, temp2);
- eval.template values<2,true,false> (temp2, hessians_quad[c][d1]);
- }
-
- // grad z: can use the values applied in x direction stored in temp1
- eval.template values<1,true,false> (temp1, temp2);
- if (evaluate_grad == true)
- eval.template gradients<2,true,false> (temp2, gradients_quad[c][d2]);
-
- // grad zz: can use the values applied in x and y direction stored
- // in temp2
- if (evaluate_lapl == true)
- eval.template hessians<2,true,false>(temp2, hessians_quad[c][d2]);
-
- // val: can use the values applied in x & y direction stored in temp2
- if (evaluate_val == true)
- eval.template values<2,true,false> (temp2, values_quad[c]);
- }
- break;
-
- default:
- AssertThrow(false, ExcNotImplemented());
- }
-
- // case additional dof for FE_Q_DG0: add values; gradients and second
- // derivatives evaluate to zero
- if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0 && evaluate_val)
- for (unsigned int c=0; c<n_components; ++c)
- for (unsigned int q=0; q<shape_info.n_q_points; ++q)
- values_quad[c][q] += values_dofs[c][shape_info.dofs_per_cell-1];
- }
-
-
-
- template <MatrixFreeFunctions::ElementType type, int dim, int fe_degree,
- int n_q_points_1d, int n_components, typename Number>
- inline
- void
- FEEvaluationImpl<type,dim,fe_degree,n_q_points_1d,n_components,Number>
- ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *scratch_data,
- const bool integrate_val,
- const bool integrate_grad)
- {
- const EvaluatorVariant variant =
- EvaluatorSelector<type,(fe_degree+n_q_points_1d>4)>::variant;
- typedef EvaluatorTensorProduct<variant, dim, fe_degree, n_q_points_1d,
- VectorizedArray<Number> > Eval;
- Eval eval (variant == evaluate_evenodd ? shape_info.shape_val_evenodd :
- shape_info.shape_values,
- variant == evaluate_evenodd ? shape_info.shape_gra_evenodd :
- shape_info.shape_gradients,
- variant == evaluate_evenodd ? shape_info.shape_hes_evenodd :
- shape_info.shape_hessians,
- shape_info.fe_degree,
- shape_info.n_q_points_1d);
-
- const unsigned int temp_size = Eval::dofs_per_cell == numbers::invalid_unsigned_int ? 0
- : (Eval::dofs_per_cell > Eval::n_q_points ?
- Eval::dofs_per_cell : Eval::n_q_points);
- VectorizedArray<Number> temp_data[(temp_size > 0 && temp_size < 100) ? 2*temp_size : 1];
- VectorizedArray<Number> *temp1;
- VectorizedArray<Number> *temp2;
- if (temp_size == 0)
- {
- temp1 = scratch_data;
- temp2 = temp1 + std::max(Utilities::fixed_power<dim>(shape_info.fe_degree+1),
- Utilities::fixed_power<dim>(shape_info.n_q_points_1d));
- }
- else if (temp_size > 100)
- {
- temp1 = scratch_data;
- temp2 = temp1 + temp_size;
- }
- else
- {
- temp1 = &temp_data[0];
- temp2 = temp1 + temp_size;
- }
-
- // expand dof_values to tensor product for truncated tensor products
- VectorizedArray<Number> **values_dofs = values_dofs_actual;
- VectorizedArray<Number> *expanded_dof_values[n_components];
- if (type == MatrixFreeFunctions::truncated_tensor)
- {
- values_dofs = expanded_dof_values;
- for (unsigned int c=0; c<n_components; ++c)
- expanded_dof_values[c] = scratch_data+2*(std::max(shape_info.dofs_per_cell,
- shape_info.n_q_points)) +
- c*Utilities::fixed_power<dim>(shape_info.fe_degree+1);
- }
-
- // These avoid compiler errors; they are only used in sensible context but
- // compilers typically cannot detect when we access something like
- // gradients_quad[2] only for dim==3.
- const unsigned int d1 = dim>1?1:0;
- const unsigned int d2 = dim>2?2:0;
-
- switch (dim)
- {
- case 1:
- for (unsigned int c=0; c<n_components; c++)
- {
- if (integrate_val == true)
- eval.template values<0,false,false> (values_quad[c], values_dofs[c]);
- if (integrate_grad == true)
- {
- if (integrate_val == true)
- eval.template gradients<0,false,true> (gradients_quad[c][0], values_dofs[c]);
- else
- eval.template gradients<0,false,false> (gradients_quad[c][0], values_dofs[c]);
- }
- }
- break;
-
- case 2:
- for (unsigned int c=0; c<n_components; c++)
- {
- if (integrate_val == true)
- {
- // val
- eval.template values<0,false,false> (values_quad[c], temp1);
- //grad x
- if (integrate_grad == true)
- eval.template gradients<0,false,true> (gradients_quad[c][0], temp1);
- eval.template values<1,false,false>(temp1, values_dofs[c]);
- }
- if (integrate_grad == true)
- {
- // grad y
- eval.template values<0,false,false> (gradients_quad[c][d1], temp1);
- if (integrate_val == false)
- {
- eval.template gradients<1,false,false>(temp1, values_dofs[c]);
- //grad x
- eval.template gradients<0,false,false> (gradients_quad[c][0], temp1);
- eval.template values<1,false,true> (temp1, values_dofs[c]);
- }
- else
- eval.template gradients<1,false,true>(temp1, values_dofs[c]);
- }
- }
- break;
-
- case 3:
- for (unsigned int c=0; c<n_components; c++)
- {
- if (integrate_val == true)
- {
- // val
- eval.template values<0,false,false> (values_quad[c], temp1);
- //grad x: can sum to temporary value in temp1
- if (integrate_grad == true)
- eval.template gradients<0,false,true> (gradients_quad[c][0], temp1);
- eval.template values<1,false,false>(temp1, temp2);
- if (integrate_grad == true)
- {
- eval.template values<0,false,false> (gradients_quad[c][d1], temp1);
- eval.template gradients<1,false,true>(temp1, temp2);
- }
- eval.template values<2,false,false> (temp2, values_dofs[c]);
- }
- else if (integrate_grad == true)
- {
- eval.template gradients<0,false,false>(gradients_quad[c][0], temp1);
- eval.template values<1,false,false> (temp1, temp2);
- eval.template values<0,false,false> (gradients_quad[c][d1], temp1);
- eval.template gradients<1,false,true>(temp1, temp2);
- eval.template values<2,false,false> (temp2, values_dofs[c]);
- }
- if (integrate_grad == true)
- {
- // grad z: can sum to temporary x and y value in output
- eval.template values<0,false,false> (gradients_quad[c][d2], temp1);
- eval.template values<1,false,false> (temp1, temp2);
- eval.template gradients<2,false,true> (temp2, values_dofs[c]);
- }
- }
- break;
-
- default:
- AssertThrow(false, ExcNotImplemented());
- }
-
- // case FE_Q_DG0: add values, gradients and second derivatives are zero
- if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0)
- {
- if (integrate_val)
- for (unsigned int c=0; c<n_components; ++c)
- {
- values_dofs[c][shape_info.dofs_per_cell-1] = values_quad[c][0];
- for (unsigned int q=1; q<shape_info.n_q_points; ++q)
- values_dofs[c][shape_info.dofs_per_cell-1] += values_quad[c][q];
- }
- else
- for (unsigned int c=0; c<n_components; ++c)
- values_dofs[c][shape_info.dofs_per_cell-1] = VectorizedArray<Number>();
- }
-
- if (type == MatrixFreeFunctions::truncated_tensor)
- {
- unsigned int count_p = 0, count_q = 0;
- const int degree = fe_degree != -1 ? fe_degree : shape_info.fe_degree;
- for (int i=0; i<(dim>2?degree+1:1); ++i)
- {
- for (int j=0; j<(dim>1?degree+1-i:1); ++j)
- {
- for (int k=0; k<degree+1-j-i; ++k, ++count_p, ++count_q)
- {
- for (unsigned int c=0; c<n_components; ++c)
- values_dofs_actual[c][count_p] = expanded_dof_values[c][count_q];
- }
- count_q += j+i;
- }
- count_q += i*(degree+1);
- }
- AssertDimension(count_q, Utilities::fixed_power<dim>(shape_info.fe_degree+1));
- }
- }
-
- // This a specialization for Gauss-Lobatto elements where the 'values'
- // operation is identity, which allows us to write shorter code.
- template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
- struct FEEvaluationImpl<MatrixFreeFunctions::tensor_gausslobatto, dim,
- fe_degree, n_q_points_1d, n_components, Number>
- {
- static
- void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
- VectorizedArray<Number> *values_dofs[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
- VectorizedArray<Number> *scratch_data,
- const bool evaluate_val,
- const bool evaluate_grad,
- const bool evaluate_lapl);
-
- static
- void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
- VectorizedArray<Number> *values_dofs[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *scratch_data,
- const bool integrate_val,
- const bool integrate_grad);
- };
-
- template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
- inline
- void
- FEEvaluationImpl<MatrixFreeFunctions::tensor_gausslobatto, dim,
- fe_degree, n_q_points_1d, n_components, Number>
- ::evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
- VectorizedArray<Number> *values_dofs[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
- VectorizedArray<Number> *scratch_data,
- const bool evaluate_val,
- const bool evaluate_grad,
- const bool evaluate_lapl)
- {
- typedef EvaluatorTensorProduct<evaluate_evenodd, dim, fe_degree, fe_degree+1,
- VectorizedArray<Number> > Eval;
- Eval eval (shape_info.shape_val_evenodd,
- shape_info.shape_gra_evenodd,
- shape_info.shape_hes_evenodd,
- shape_info.fe_degree,
- shape_info.n_q_points_1d);
-
- // These avoid compiler errors; they are only used in sensible context but
- // compilers typically cannot detect when we access something like
- // gradients_quad[2] only for dim==3.
- const unsigned int d1 = dim>1?1:0;
- const unsigned int d2 = dim>2?2:0;
- const unsigned int d3 = dim>2?3:0;
- const unsigned int d4 = dim>2?4:0;
- const unsigned int d5 = dim>2?5:0;
-
- switch (dim)
- {
- case 1:
- if (evaluate_val == true)
- std::memcpy (values_quad[0], values_dofs[0],
- eval.dofs_per_cell * n_components *
- sizeof (values_dofs[0][0]));
- for (unsigned int c=0; c<n_components; c++)
- {
- if (evaluate_grad == true)
- eval.template gradients<0,true,false>(values_dofs[c], gradients_quad[c][0]);
- if (evaluate_lapl == true)
- eval.template hessians<0,true,false> (values_dofs[c], hessians_quad[c][0]);
- }
- break;
-
- case 2:
- if (evaluate_val == true)
- {
- std::memcpy (values_quad[0], values_dofs[0],
- Eval::dofs_per_cell * n_components *
- sizeof (values_dofs[0][0]));
- }
- if (evaluate_grad == true)
- for (unsigned int comp=0; comp<n_components; comp++)
- {
- // grad x
- eval.template gradients<0,true,false> (values_dofs[comp],
- gradients_quad[comp][0]);
- // grad y
- eval.template gradients<1,true,false> (values_dofs[comp],
- gradients_quad[comp][d1]);
- }
- if (evaluate_lapl == true)
- for (unsigned int comp=0; comp<n_components; comp++)
- {
- // hess x
- eval.template hessians<0,true,false> (values_dofs[comp],
- hessians_quad[comp][0]);
- // hess y
- eval.template hessians<1,true,false> (values_dofs[comp],
- hessians_quad[comp][d1]);
-
- // grad x grad y
- eval.template gradients<0,true,false> (values_dofs[comp], scratch_data);
- eval.template gradients<1,true,false> (scratch_data, hessians_quad[comp][d1+d1]);
- }
- break;
-
- case 3:
- if (evaluate_val == true)
- {
- std::memcpy (values_quad[0], values_dofs[0],
- Eval::dofs_per_cell * n_components *
- sizeof (values_dofs[0][0]));
- }
- if (evaluate_grad == true)
- for (unsigned int comp=0; comp<n_components; comp++)
- {
- // grad x
- eval.template gradients<0,true,false> (values_dofs[comp],
- gradients_quad[comp][0]);
- // grad y
- eval.template gradients<1,true,false> (values_dofs[comp],
- gradients_quad[comp][d1]);
- // grad y
- eval.template gradients<2,true,false> (values_dofs[comp],
- gradients_quad[comp][d2]);
- }
- if (evaluate_lapl == true)
- for (unsigned int comp=0; comp<n_components; comp++)
- {
- // grad x
- eval.template hessians<0,true,false> (values_dofs[comp],
- hessians_quad[comp][0]);
- // grad y
- eval.template hessians<1,true,false> (values_dofs[comp],
- hessians_quad[comp][d1]);
- // grad y
- eval.template hessians<2,true,false> (values_dofs[comp],
- hessians_quad[comp][d2]);
-
- VectorizedArray<Number> *temp1 = scratch_data;
- // grad xy
- eval.template gradients<0,true,false> (values_dofs[comp], temp1);
- eval.template gradients<1,true,false> (temp1, hessians_quad[comp][d3]);
- // grad xz
- eval.template gradients<2,true,false> (temp1, hessians_quad[comp][d4]);
- // grad yz
- eval.template gradients<1,true,false> (values_dofs[comp], temp1);
- eval.template gradients<2,true,false> (temp1, hessians_quad[comp][d5]);
- }
- break;
- default:
- AssertThrow(false, ExcNotImplemented());
- }
- }
-
- template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
- inline
- void
- FEEvaluationImpl<MatrixFreeFunctions::tensor_gausslobatto, dim,
- fe_degree, n_q_points_1d, n_components, Number>
- ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
- VectorizedArray<Number> *values_dofs[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *,
- const bool integrate_val,
- const bool integrate_grad)
- {
- typedef EvaluatorTensorProduct<evaluate_evenodd, dim, fe_degree, fe_degree+1,
- VectorizedArray<Number> > Eval;
- Eval eval (shape_info.shape_val_evenodd,
- shape_info.shape_gra_evenodd,
- shape_info.shape_hes_evenodd,
- shape_info.fe_degree,
- shape_info.n_q_points_1d);
-
- // These avoid compiler errors; they are only used in sensible context but
- // compilers typically cannot detect when we access something like
- // gradients_quad[2] only for dim==3.
- const unsigned int d1 = dim>1?1:0;
- const unsigned int d2 = dim>2?2:0;
-
- if (integrate_val == true)
- std::memcpy (values_dofs[0], values_quad[0],
- Eval::dofs_per_cell * n_components *
- sizeof (values_dofs[0][0]));
- switch (dim)
- {
- case 1:
- for (unsigned int c=0; c<n_components; c++)
- {
- if (integrate_grad == true)
- {
- if (integrate_val == true)
- eval.template gradients<0,false,true> (gradients_quad[c][0],
- values_dofs[c]);
- else
- eval.template gradients<0,false,false> (gradients_quad[c][0],
- values_dofs[c]);
- }
- }
-
- break;
- case 2:
- if (integrate_grad == true)
- for (unsigned int comp=0; comp<n_components; comp++)
- {
- // grad x: If integrate_val == true we have to add to the
- // previous output
- if (integrate_val == true)
- eval.template gradients<0, false, true> (gradients_quad[comp][0],
- values_dofs[comp]);
- else
- eval.template gradients<0, false, false> (gradients_quad[comp][0],
- values_dofs[comp]);
-
- // grad y
- eval.template gradients<1, false, true> (gradients_quad[comp][d1],
- values_dofs[comp]);
- }
- break;
-
- case 3:
- if (integrate_grad == true)
- for (unsigned int comp=0; comp<n_components; comp++)
- {
- // grad x: If integrate_val == true we have to add to the
- // previous output
- if (integrate_val == true)
- eval.template gradients<0, false, true> (gradients_quad[comp][0],
- values_dofs[comp]);
- else
- eval.template gradients<0, false, false> (gradients_quad[comp][0],
- values_dofs[comp]);
-
- // grad y
- eval.template gradients<1, false, true> (gradients_quad[comp][d1],
- values_dofs[comp]);
-
- // grad z
- eval.template gradients<2, false, true> (gradients_quad[comp][d2],
- values_dofs[comp]);
- }
- break;
-
- default:
- AssertThrow(false, ExcNotImplemented());
- }
- }
-
-} // end of namespace internal
-
-
/*-------------------------- FEEvaluation -----------------------------------*/
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii__matrix_free_tensor_product_kernels_h
+#define dealii__matrix_free_tensor_product_kernels_h
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/aligned_vector.h>
+#include <deal.II/base/utilities.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+
+namespace internal
+{
+ /**
+ * In this namespace, the evaluator routines that evaluate the tensor
+ * products are implemented.
+ */
+ enum EvaluatorVariant
+ {
+ /**
+ * Do not use anything more than the tensor product structure of the
+ * finite element.
+ */
+ evaluate_general,
+ /**
+ * Perform evaluation by exploiting symmetry in the finite element: i.e.,
+ * skip some computations by utilizing the symmetry in the shape functions
+ * and quadrature points.
+ */
+ evaluate_symmetric,
+ /**
+ * Use symmetry to apply the operator to even and odd parts of the input
+ * vector separately: see the documentation of the EvaluatorTensorProduct
+ * specialization for more information.
+ */
+ evaluate_evenodd
+ };
+
+ /**
+ * Generic evaluator framework
+ */
+ template <EvaluatorVariant variant, int dim, int fe_degree, int n_q_points_1d,
+ typename Number>
+ struct EvaluatorTensorProduct
+ {};
+
+ /**
+ * Internal evaluator for 1d-3d shape function using the tensor product form
+ * of the basis functions
+ */
+ template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+ struct EvaluatorTensorProduct<evaluate_general,dim,fe_degree,n_q_points_1d,Number>
+ {
+ static const unsigned int dofs_per_cell = Utilities::fixed_int_power<fe_degree+1,dim>::value;
+ static const unsigned int n_q_points = Utilities::fixed_int_power<n_q_points_1d,dim>::value;
+
+ /**
+ * Empty constructor. Does nothing. Be careful when using 'values' and
+ * related methods because they need to be filled with the other pointer
+ */
+ EvaluatorTensorProduct ()
+ :
+ shape_values (0),
+ shape_gradients (0),
+ shape_hessians (0)
+ {}
+
+ /**
+ * Constructor, taking the data from ShapeInfo
+ */
+ EvaluatorTensorProduct (const AlignedVector<Number> &shape_values,
+ const AlignedVector<Number> &shape_gradients,
+ const AlignedVector<Number> &shape_hessians,
+ const unsigned int dummy1 = 0,
+ const unsigned int dummy2 = 0)
+ :
+ shape_values (shape_values.begin()),
+ shape_gradients (shape_gradients.begin()),
+ shape_hessians (shape_hessians.begin())
+ {
+ (void)dummy1;
+ (void)dummy2;
+ }
+
+ template <int direction, bool dof_to_quad, bool add>
+ void
+ values (const Number in [],
+ Number out[]) const
+ {
+ apply<direction,dof_to_quad,add>(shape_values, in, out);
+ }
+
+ template <int direction, bool dof_to_quad, bool add>
+ void
+ gradients (const Number in [],
+ Number out[]) const
+ {
+ apply<direction,dof_to_quad,add>(shape_gradients, in, out);
+ }
+
+ template <int direction, bool dof_to_quad, bool add>
+ void
+ hessians (const Number in [],
+ Number out[]) const
+ {
+ apply<direction,dof_to_quad,add>(shape_hessians, in, out);
+ }
+
+ template <int direction, bool dof_to_quad, bool add>
+ static void apply (const Number *shape_data,
+ const Number in [],
+ Number out []);
+
+ const Number *shape_values;
+ const Number *shape_gradients;
+ const Number *shape_hessians;
+ };
+
+ // evaluates the given shape data in 1d-3d using the tensor 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>
+ template <int direction, bool dof_to_quad, bool add>
+ inline
+ void
+ EvaluatorTensorProduct<evaluate_general,dim,fe_degree,n_q_points_1d,Number>
+ ::apply (const Number *shape_data,
+ const Number in [],
+ 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_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<nn; ++col)
+ {
+ Number val0;
+ if (dof_to_quad == true)
+ val0 = shape_data[col];
+ else
+ val0 = shape_data[col*n_q_points_1d];
+ Number res0 = val0 * in[0];
+ for (int ind=1; ind<mm; ++ind)
+ {
+ if (dof_to_quad == true)
+ val0 = shape_data[ind*n_q_points_1d+col];
+ else
+ val0 = shape_data[col*n_q_points_1d+ind];
+ res0 += val0 * in[stride*ind];
+ }
+ if (add == false)
+ out[stride*col] = res0;
+ else
+ out[stride*col] += 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);
+ }
+ }
+ }
+
+
+
+ // 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:
+ ++in;
+ ++out;
+ // faces 2 and 3 in 3D use local coordinate system zx, which
+ // is the other way around compared to the tensor
+ // product. Need to take that into account.
+ if (dim == 3)
+ {
+ if (dof_to_quad)
+ out += fe_degree;
+ else
+ in += fe_degree;
+ }
+ break;
+ case 2:
+ ++in;
+ ++out;
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+ }
+ if (face_direction == 1 && dim == 3)
+ {
+ in += mm*(mm-1);
+ out += nn*(nn-1);
+ // adjust for local coordinate system zx
+ if (dof_to_quad)
+ out -= (fe_degree+1)*(fe_degree+1)-1;
+ else
+ in -= (fe_degree+1)*(fe_degree+1)-1;
+ }
+ }
+ }
+
+
+
+ /**
+ * Internal evaluator for 1d-3d shape function using the tensor product form
+ * of the basis functions. The same as above but without making use of
+ * template arguments and rather variable loop bounds.
+ */
+ template <int dim, typename Number>
+ struct EvaluatorTensorProduct<evaluate_general,dim,-1,0,Number>
+ {
+ static const unsigned int dofs_per_cell = numbers::invalid_unsigned_int;
+ static const unsigned int n_q_points = numbers::invalid_unsigned_int;
+
+ /**
+ * Empty constructor. Does nothing. Be careful when using 'values' and
+ * related methods because they need to be filled with the other constructor
+ */
+ EvaluatorTensorProduct ()
+ :
+ shape_values (0),
+ shape_gradients (0),
+ shape_hessians (0),
+ fe_degree (numbers::invalid_unsigned_int),
+ n_q_points_1d (numbers::invalid_unsigned_int)
+ {}
+
+ /**
+ * Constructor, taking the data from ShapeInfo
+ */
+ EvaluatorTensorProduct (const AlignedVector<Number> &shape_values,
+ const AlignedVector<Number> &shape_gradients,
+ const AlignedVector<Number> &shape_hessians,
+ const unsigned int fe_degree,
+ const unsigned int n_q_points_1d)
+ :
+ shape_values (shape_values.begin()),
+ shape_gradients (shape_gradients.begin()),
+ shape_hessians (shape_hessians.begin()),
+ fe_degree (fe_degree),
+ n_q_points_1d (n_q_points_1d)
+ {}
+
+ template <int direction, bool dof_to_quad, bool add>
+ void
+ values (const Number *in,
+ Number *out) const
+ {
+ apply<direction,dof_to_quad,add>(shape_values, in, out);
+ }
+
+ template <int direction, bool dof_to_quad, bool add>
+ void
+ gradients (const Number *in,
+ Number *out) const
+ {
+ apply<direction,dof_to_quad,add>(shape_gradients, in, out);
+ }
+
+ template <int direction, bool dof_to_quad, bool add>
+ void
+ hessians (const Number *in,
+ Number *out) const
+ {
+ apply<direction,dof_to_quad,add>(shape_hessians, in, out);
+ }
+
+ template <int direction, bool dof_to_quad, bool add>
+ void apply (const Number *shape_data,
+ const Number *in,
+ Number *out) const;
+
+ const Number *shape_values;
+ const Number *shape_gradients;
+ const Number *shape_hessians;
+ const unsigned int fe_degree;
+ const unsigned int n_q_points_1d;
+ };
+
+ // evaluates the given shape data in 1d-3d using the tensor 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, typename Number>
+ template <int direction, bool dof_to_quad, bool add>
+ inline
+ void
+ EvaluatorTensorProduct<evaluate_general,dim,-1,0,Number>
+ ::apply (const Number *shape_data,
+ const Number *in,
+ Number *out) const
+ {
+ 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_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
+ const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
+ const int stride = direction==0 ? 1 : Utilities::fixed_power<direction>(nn);
+
+ for (int i2=0; i2<n_blocks2; ++i2)
+ {
+ for (int i1=0; i1<n_blocks1; ++i1)
+ {
+ for (int col=0; col<nn; ++col)
+ {
+ Number val0;
+ if (dof_to_quad == true)
+ val0 = shape_data[col];
+ else
+ val0 = shape_data[col*n_q_points_1d];
+ Number res0 = val0 * in[0];
+ for (int ind=1; ind<mm; ++ind)
+ {
+ if (dof_to_quad == true)
+ val0 = shape_data[ind*n_q_points_1d+col];
+ else
+ val0 = shape_data[col*n_q_points_1d+ind];
+ res0 += val0 * in[stride*ind];
+ }
+ if (add == false)
+ out[stride*col] = res0;
+ else
+ out[stride*col] += 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 evaluator for 1d-3d shape function using the tensor product form
+ * of the basis functions. This class 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.
+ */
+ template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+ struct EvaluatorTensorProduct<evaluate_symmetric,dim,fe_degree,n_q_points_1d,Number>
+ {
+ static const unsigned int dofs_per_cell = Utilities::fixed_int_power<fe_degree+1,dim>::value;
+ static const unsigned int n_q_points = Utilities::fixed_int_power<n_q_points_1d,dim>::value;
+
+ /**
+ * Constructor, taking the data from ShapeInfo
+ */
+ EvaluatorTensorProduct (const AlignedVector<Number> &shape_values,
+ const AlignedVector<Number> &shape_gradients,
+ const AlignedVector<Number> &shape_hessians,
+ const unsigned int dummy1 = 0,
+ const unsigned int dummy2 = 0)
+ :
+ shape_values (shape_values.begin()),
+ shape_gradients (shape_gradients.begin()),
+ shape_hessians (shape_hessians.begin())
+ {
+ (void)dummy1;
+ (void)dummy2;
+ }
+
+ template <int direction, bool dof_to_quad, bool add>
+ void
+ values (const Number in [],
+ Number out[]) const;
+
+ template <int direction, bool dof_to_quad, bool add>
+ void
+ gradients (const Number in [],
+ Number out[]) const;
+
+ template <int direction, bool dof_to_quad, bool add>
+ void
+ hessians (const Number in [],
+ Number out[]) const;
+
+ const Number *shape_values;
+ const Number *shape_gradients;
+ const Number *shape_hessians;
+ };
+
+
+
+ // In this 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>
+ template <int direction, bool dof_to_quad, bool add>
+ inline
+ void
+ EvaluatorTensorProduct<evaluate_symmetric,dim,fe_degree,n_q_points_1d,Number>
+ ::values (const Number in [],
+ Number out []) const
+ {
+ 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)
+ {
+ 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 = 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)
+ {
+ Number res0;
+ 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];
+ Number val1 = in[stride*ind] + in[stride*(mm-1-ind)];
+ val1 *= val0;
+ res0 += val1;
+ }
+ }
+ else
+ res0 = Number();
+ if (add == false)
+ out[stride*n_cols] = res0;
+ else
+ out[stride*n_cols] += res0;
+ }
+ else if (dof_to_quad == false && nn%2 == 1)
+ {
+ Number res0;
+ if (mid > 0)
+ {
+ 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];
+ 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);
+ }
+ }
+ }
+
+
+
+ // 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>
+ template <int direction, bool dof_to_quad, bool add>
+ inline
+ void
+ EvaluatorTensorProduct<evaluate_symmetric,dim,fe_degree,n_q_points_1d,Number>
+ ::gradients (const Number in [],
+ Number out []) const
+ {
+ 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)
+ {
+ 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 = 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 )
+ {
+ 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];
+ 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>
+ template <int direction, bool dof_to_quad, bool add>
+ inline
+ void
+ EvaluatorTensorProduct<evaluate_symmetric,dim,fe_degree,n_q_points_1d,Number>
+ ::hessians (const Number in [],
+ Number out []) const
+ {
+ 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)
+ {
+ 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 = 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 )
+ {
+ 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];
+ Number val1 = in[stride*ind] + in[stride*(mm-1-ind)];
+ val1 *= val0;
+ res0 += val1;
+ }
+ }
+ else
+ res0 = 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);
+ }
+ }
+ }
+
+
+
+ /**
+ * Internal evaluator for 1d-3d shape function using the tensor product form
+ * of the basis functions.
+ *
+ * This class 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>
+ struct EvaluatorTensorProduct<evaluate_evenodd,dim,fe_degree,n_q_points_1d,Number>
+ {
+ static const unsigned int dofs_per_cell = Utilities::fixed_int_power<fe_degree+1,dim>::value;
+ static const unsigned int n_q_points = Utilities::fixed_int_power<n_q_points_1d,dim>::value;
+
+ /**
+ * Empty constructor. Does nothing. Be careful when using 'values' and
+ * related methods because they need to be filled with the other pointer
+ */
+ EvaluatorTensorProduct ()
+ :
+ shape_values (0),
+ shape_gradients (0),
+ shape_hessians (0)
+ {}
+
+ /**
+ * Constructor, taking the data from ShapeInfo (using the even-odd
+ * variants stored there)
+ */
+ EvaluatorTensorProduct (const AlignedVector<Number> &shape_values,
+ const AlignedVector<Number> &shape_gradients,
+ const AlignedVector<Number> &shape_hessians,
+ const unsigned int dummy1 = 0,
+ const unsigned int dummy2 = 0)
+ :
+ shape_values (shape_values.begin()),
+ shape_gradients (shape_gradients.begin()),
+ shape_hessians (shape_hessians.begin())
+ {
+ (void)dummy1;
+ (void)dummy2;
+ }
+
+ template <int direction, bool dof_to_quad, bool add>
+ void
+ values (const Number in [],
+ Number out[]) const
+ {
+ apply<direction,dof_to_quad,add,0>(shape_values, in, out);
+ }
+
+ template <int direction, bool dof_to_quad, bool add>
+ void
+ gradients (const Number in [],
+ Number out[]) const
+ {
+ apply<direction,dof_to_quad,add,1>(shape_gradients, in, out);
+ }
+
+ template <int direction, bool dof_to_quad, bool add>
+ void
+ hessians (const Number in [],
+ Number out[]) const
+ {
+ apply<direction,dof_to_quad,add,2>(shape_hessians, in, out);
+ }
+
+ template <int direction, bool dof_to_quad, bool add, int type>
+ static void apply (const Number *shape_data,
+ const Number in [],
+ Number out []);
+
+ const Number *shape_values;
+ const Number *shape_gradients;
+ const Number *shape_hessians;
+ };
+
+
+
+ template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+ template <int direction, bool dof_to_quad, bool add, int type>
+ inline
+ void
+ EvaluatorTensorProduct<evaluate_evenodd,dim,fe_degree,n_q_points_1d,Number>
+ ::apply (const Number *shapes,
+ 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;
+
+ const int offset = (n_q_points_1d+1)/2;
+
+ // 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>0?mid:1], xm[mid>0?mid:1];
+ 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[col] * xp[0];
+ r1 = shapes[fe_degree*offset + col] * xm[0];
+ }
+ else
+ {
+ r0 = shapes[col*offset] * xp[0];
+ r1 = shapes[(fe_degree-col)*offset] * xm[0];
+ }
+ for (int ind=1; ind<mid; ++ind)
+ {
+ if (dof_to_quad == true)
+ {
+ r0 += shapes[ind*offset+col] * xp[ind];
+ r1 += shapes[(fe_degree-ind)*offset+col] * xm[ind];
+ }
+ else
+ {
+ r0 += shapes[col*offset+ind] * xp[ind];
+ r1 += shapes[(fe_degree-col)*offset+ind] * xm[ind];
+ }
+ }
+ }
+ else
+ r0 = r1 = Number();
+ if (mm % 2 == 1 && dof_to_quad == true)
+ {
+ if (type == 1)
+ r1 += shapes[mid*offset+col] * in[stride*mid];
+ else
+ r0 += shapes[mid*offset+col] * in[stride*mid];
+ }
+ else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0))
+ r0 += shapes[col*offset+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[n_cols] * xp[0];
+ for (int ind=1; ind<mid; ++ind)
+ r0 += shapes[ind*offset+n_cols] * xp[ind];
+ }
+ else
+ r0 = Number();
+ if (type != 1 && mm % 2 == 1)
+ r0 += shapes[mid*offset+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*offset] * xm[0];
+ for (int ind=1; ind<mid; ++ind)
+ r0 += shapes[n_cols*offset+ind] * xm[ind];
+ }
+ else
+ {
+ r0 = shapes[n_cols*offset] * xp[0];
+ for (int ind=1; ind<mid; ++ind)
+ r0 += shapes[n_cols*offset+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*offset+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);
+ }
+ }
+ }
+
+} // end of namespace internal
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
#include <deal.II/multigrid/mg_transfer_matrix_free.h>
#include <deal.II/multigrid/mg_transfer_internal.h>
-#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/tensor_product_kernels.h>
#include <algorithm>
projects / dealii.git / commitdiff