From: Martin Kronbichler Date: Thu, 2 Mar 2017 09:02:27 +0000 (+0100) Subject: Move tensor product kernels from fe_evaluation.h into separate file X-Git-Tag: v8.5.0-rc1~79^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F4018%2Fhead;p=dealii.git Move tensor product kernels from fe_evaluation.h into separate file --- diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h new file mode 100644 index 0000000000..15c0ac4b08 --- /dev/null +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -0,0 +1,787 @@ +// --------------------------------------------------------------------- +// +// 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 +#include +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + + + +namespace internal +{ + // Select evaluator type from element shape function type + template + struct EvaluatorSelector {}; + + template + struct EvaluatorSelector + { + static const EvaluatorVariant variant = evaluate_general; + }; + + template <> + struct EvaluatorSelector + { + static const EvaluatorVariant variant = evaluate_symmetric; + }; + + template <> struct EvaluatorSelector + { + static const EvaluatorVariant variant = evaluate_evenodd; + }; + + template + struct EvaluatorSelector + { + static const EvaluatorVariant variant = evaluate_general; + }; + + template <> + struct EvaluatorSelector + { + static const EvaluatorVariant variant = evaluate_general; + }; + + template <> + struct EvaluatorSelector + { + static const EvaluatorVariant variant = evaluate_evenodd; + }; + + template + struct EvaluatorSelector + { + 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 + struct FEEvaluationImpl + { + static + void evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, + VectorizedArray *values_dofs_actual[], + VectorizedArray *values_quad[], + VectorizedArray *gradients_quad[][dim], + VectorizedArray *hessians_quad[][(dim*(dim+1))/2], + VectorizedArray *scratch_data, + const bool evaluate_val, + const bool evaluate_grad, + const bool evaluate_lapl); + + static + void integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, + VectorizedArray *values_dofs_actual[], + VectorizedArray *values_quad[], + VectorizedArray *gradients_quad[][dim], + VectorizedArray *scratch_data, + const bool evaluate_val, + const bool evaluate_grad); + }; + + + template + inline + void + FEEvaluationImpl + ::evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, + VectorizedArray *values_dofs_actual[], + VectorizedArray *values_quad[], + VectorizedArray *gradients_quad[][dim], + VectorizedArray *hessians_quad[][(dim*(dim+1))/2], + VectorizedArray *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 = + EvaluatorSelector4)>::variant; + typedef EvaluatorTensorProduct > 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 temp_data[(temp_size > 0 && temp_size < 100) ? 2*temp_size : 1]; + VectorizedArray *temp1; + VectorizedArray *temp2; + if (temp_size == 0) + { + temp1 = scratch_data; + temp2 = temp1 + std::max(Utilities::fixed_power(shape_info.fe_degree+1), + Utilities::fixed_power(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 **values_dofs = values_dofs_actual; + VectorizedArray *expanded_dof_values[n_components]; + if (type == MatrixFreeFunctions::truncated_tensor) + { + values_dofs = expanded_dof_values; + for (unsigned int c=0; c(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(); + } + for (int j=degree+1-i; j(); + } + AssertDimension(count_q, Utilities::fixed_power(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 (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 (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 (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 + inline + void + FEEvaluationImpl + ::integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, + VectorizedArray *values_dofs_actual[], + VectorizedArray *values_quad[], + VectorizedArray *gradients_quad[][dim], + VectorizedArray *scratch_data, + const bool integrate_val, + const bool integrate_grad) + { + const EvaluatorVariant variant = + EvaluatorSelector4)>::variant; + typedef EvaluatorTensorProduct > 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 temp_data[(temp_size > 0 && temp_size < 100) ? 2*temp_size : 1]; + VectorizedArray *temp1; + VectorizedArray *temp2; + if (temp_size == 0) + { + temp1 = scratch_data; + temp2 = temp1 + std::max(Utilities::fixed_power(shape_info.fe_degree+1), + Utilities::fixed_power(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 **values_dofs = values_dofs_actual; + VectorizedArray *expanded_dof_values[n_components]; + if (type == MatrixFreeFunctions::truncated_tensor) + { + values_dofs = expanded_dof_values; + for (unsigned int c=0; c(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 (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 (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 (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(); + } + + 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(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 + struct FEEvaluationImpl + { + static + void evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, + VectorizedArray *values_dofs[], + VectorizedArray *values_quad[], + VectorizedArray *gradients_quad[][dim], + VectorizedArray *hessians_quad[][(dim*(dim+1))/2], + VectorizedArray *scratch_data, + const bool evaluate_val, + const bool evaluate_grad, + const bool evaluate_lapl); + + static + void integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, + VectorizedArray *values_dofs[], + VectorizedArray *values_quad[], + VectorizedArray *gradients_quad[][dim], + VectorizedArray *scratch_data, + const bool integrate_val, + const bool integrate_grad); + }; + + template + inline + void + FEEvaluationImpl + ::evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, + VectorizedArray *values_dofs[], + VectorizedArray *values_quad[], + VectorizedArray *gradients_quad[][dim], + VectorizedArray *hessians_quad[][(dim*(dim+1))/2], + VectorizedArray *scratch_data, + const bool evaluate_val, + const bool evaluate_grad, + const bool evaluate_lapl) + { + typedef EvaluatorTensorProduct > 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(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 (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 (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 (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 (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 *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 + inline + void + FEEvaluationImpl + ::integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, + VectorizedArray *values_dofs[], + VectorizedArray *values_quad[], + VectorizedArray *gradients_quad[][dim], + VectorizedArray *, + const bool integrate_val, + const bool integrate_grad) + { + typedef EvaluatorTensorProduct > 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 (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 (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 (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 diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 0213174a93..7945b28a32 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -28,6 +28,8 @@ #include #include #include +#include +#include DEAL_II_NAMESPACE_OPEN @@ -4995,2024 +4997,6 @@ FEEvaluationAccess<1,1,Number> -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 - struct EvaluatorTensorProduct - {}; - - /** - * Internal evaluator for 1d-3d shape function using the tensor product form - * of the basis functions - */ - template - struct EvaluatorTensorProduct - { - static const unsigned int dofs_per_cell = Utilities::fixed_int_power::value; - static const unsigned int n_q_points = Utilities::fixed_int_power::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 &shape_values, - const AlignedVector &shape_gradients, - const AlignedVector &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 - void - values (const Number in [], - Number out[]) const - { - apply(shape_values, in, out); - } - - template - void - gradients (const Number in [], - Number out[]) const - { - apply(shape_gradients, in, out); - } - - template - void - hessians (const Number in [], - Number out[]) const - { - apply(shape_hessians, in, out); - } - - template - 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 - template - inline - void - EvaluatorTensorProduct - ::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::value; - - for (int i2=0; i2 - 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::value; - - for (int i2=0; i2 - struct EvaluatorTensorProduct - { - 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 &shape_values, - const AlignedVector &shape_gradients, - const AlignedVector &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 - void - values (const Number *in, - Number *out) const - { - apply(shape_values, in, out); - } - - template - void - gradients (const Number *in, - Number *out) const - { - apply(shape_gradients, in, out); - } - - template - void - hessians (const Number *in, - Number *out) const - { - apply(shape_hessians, in, out); - } - - template - 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 - template - inline - void - EvaluatorTensorProduct - ::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(nn); - - for (int i2=0; i2 - struct EvaluatorTensorProduct - { - static const unsigned int dofs_per_cell = Utilities::fixed_int_power::value; - static const unsigned int n_q_points = Utilities::fixed_int_power::value; - - /** - * Constructor, taking the data from ShapeInfo - */ - EvaluatorTensorProduct (const AlignedVector &shape_values, - const AlignedVector &shape_gradients, - const AlignedVector &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 - void - values (const Number in [], - Number out[]) const; - - template - void - gradients (const Number in [], - Number out[]) const; - - template - 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 - template - inline - void - EvaluatorTensorProduct - ::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::value; - - for (int i2=0; i2 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 0) - { - res0 = in[0] + in[stride*(mm-1)]; - res0 *= val0; - for (int ind=1; ind 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 [-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 - template - inline - void - EvaluatorTensorProduct - ::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::value; - - for (int i2=0; i2 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 - template - inline - void - EvaluatorTensorProduct - ::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::value; - - for (int i2=0; i2 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 0) - { - res0 = in[0] + in[stride*(mm-1)]; - res0 *= val0; - for (int ind=1; ind - struct EvaluatorTensorProduct - { - static const unsigned int dofs_per_cell = Utilities::fixed_int_power::value; - static const unsigned int n_q_points = Utilities::fixed_int_power::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 &shape_values, - const AlignedVector &shape_gradients, - const AlignedVector &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 - void - values (const Number in [], - Number out[]) const - { - apply(shape_values, in, out); - } - - template - void - gradients (const Number in [], - Number out[]) const - { - apply(shape_gradients, in, out); - } - - template - void - hessians (const Number in [], - Number out[]) const - { - apply(shape_hessians, in, out); - } - - template - 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 - template - inline - void - EvaluatorTensorProduct - ::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::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; i20?mid:1], xm[mid>0?mid:1]; - for (int i=0; i 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 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 0) - { - if (type == 1) - { - r0 = shapes[n_cols*offset] * xm[0]; - for (int ind=1; ind - struct EvaluatorSelector {}; - - template - struct EvaluatorSelector - { - static const EvaluatorVariant variant = evaluate_general; - }; - - template <> - struct EvaluatorSelector - { - static const EvaluatorVariant variant = evaluate_symmetric; - }; - - template <> struct EvaluatorSelector - { - static const EvaluatorVariant variant = evaluate_evenodd; - }; - - template - struct EvaluatorSelector - { - static const EvaluatorVariant variant = evaluate_general; - }; - - template <> - struct EvaluatorSelector - { - static const EvaluatorVariant variant = evaluate_general; - }; - - template <> - struct EvaluatorSelector - { - static const EvaluatorVariant variant = evaluate_evenodd; - }; - - template - struct EvaluatorSelector - { - 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 - struct FEEvaluationImpl - { - static - void evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *hessians_quad[][(dim*(dim+1))/2], - VectorizedArray *scratch_data, - const bool evaluate_val, - const bool evaluate_grad, - const bool evaluate_lapl); - - static - void integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *scratch_data, - const bool evaluate_val, - const bool evaluate_grad); - }; - - - template - inline - void - FEEvaluationImpl - ::evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *hessians_quad[][(dim*(dim+1))/2], - VectorizedArray *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 = - EvaluatorSelector4)>::variant; - typedef EvaluatorTensorProduct > 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 temp_data[(temp_size > 0 && temp_size < 100) ? 2*temp_size : 1]; - VectorizedArray *temp1; - VectorizedArray *temp2; - if (temp_size == 0) - { - temp1 = scratch_data; - temp2 = temp1 + std::max(Utilities::fixed_power(shape_info.fe_degree+1), - Utilities::fixed_power(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 **values_dofs = values_dofs_actual; - VectorizedArray *expanded_dof_values[n_components]; - if (type == MatrixFreeFunctions::truncated_tensor) - { - values_dofs = expanded_dof_values; - for (unsigned int c=0; c(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(); - } - for (int j=degree+1-i; j(); - } - AssertDimension(count_q, Utilities::fixed_power(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 (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 (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 (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 - inline - void - FEEvaluationImpl - ::integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *scratch_data, - const bool integrate_val, - const bool integrate_grad) - { - const EvaluatorVariant variant = - EvaluatorSelector4)>::variant; - typedef EvaluatorTensorProduct > 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 temp_data[(temp_size > 0 && temp_size < 100) ? 2*temp_size : 1]; - VectorizedArray *temp1; - VectorizedArray *temp2; - if (temp_size == 0) - { - temp1 = scratch_data; - temp2 = temp1 + std::max(Utilities::fixed_power(shape_info.fe_degree+1), - Utilities::fixed_power(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 **values_dofs = values_dofs_actual; - VectorizedArray *expanded_dof_values[n_components]; - if (type == MatrixFreeFunctions::truncated_tensor) - { - values_dofs = expanded_dof_values; - for (unsigned int c=0; c(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 (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 (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 (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(); - } - - 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(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 - struct FEEvaluationImpl - { - static - void evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, - VectorizedArray *values_dofs[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *hessians_quad[][(dim*(dim+1))/2], - VectorizedArray *scratch_data, - const bool evaluate_val, - const bool evaluate_grad, - const bool evaluate_lapl); - - static - void integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, - VectorizedArray *values_dofs[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *scratch_data, - const bool integrate_val, - const bool integrate_grad); - }; - - template - inline - void - FEEvaluationImpl - ::evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, - VectorizedArray *values_dofs[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *hessians_quad[][(dim*(dim+1))/2], - VectorizedArray *scratch_data, - const bool evaluate_val, - const bool evaluate_grad, - const bool evaluate_lapl) - { - typedef EvaluatorTensorProduct > 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(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 (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 (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 (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 (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 *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 - inline - void - FEEvaluationImpl - ::integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, - VectorizedArray *values_dofs[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *, - const bool integrate_val, - const bool integrate_grad) - { - typedef EvaluatorTensorProduct > 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 (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 (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 (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 -----------------------------------*/ diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h new file mode 100644 index 0000000000..3e3683ecb2 --- /dev/null +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -0,0 +1,1298 @@ +// --------------------------------------------------------------------- +// +// 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 +#include +#include + + +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 + struct EvaluatorTensorProduct + {}; + + /** + * Internal evaluator for 1d-3d shape function using the tensor product form + * of the basis functions + */ + template + struct EvaluatorTensorProduct + { + static const unsigned int dofs_per_cell = Utilities::fixed_int_power::value; + static const unsigned int n_q_points = Utilities::fixed_int_power::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 &shape_values, + const AlignedVector &shape_gradients, + const AlignedVector &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 + void + values (const Number in [], + Number out[]) const + { + apply(shape_values, in, out); + } + + template + void + gradients (const Number in [], + Number out[]) const + { + apply(shape_gradients, in, out); + } + + template + void + hessians (const Number in [], + Number out[]) const + { + apply(shape_hessians, in, out); + } + + template + 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 + template + inline + void + EvaluatorTensorProduct + ::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::value; + + for (int i2=0; i2 + 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::value; + + for (int i2=0; i2 + struct EvaluatorTensorProduct + { + 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 &shape_values, + const AlignedVector &shape_gradients, + const AlignedVector &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 + void + values (const Number *in, + Number *out) const + { + apply(shape_values, in, out); + } + + template + void + gradients (const Number *in, + Number *out) const + { + apply(shape_gradients, in, out); + } + + template + void + hessians (const Number *in, + Number *out) const + { + apply(shape_hessians, in, out); + } + + template + 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 + template + inline + void + EvaluatorTensorProduct + ::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(nn); + + for (int i2=0; i2 + struct EvaluatorTensorProduct + { + static const unsigned int dofs_per_cell = Utilities::fixed_int_power::value; + static const unsigned int n_q_points = Utilities::fixed_int_power::value; + + /** + * Constructor, taking the data from ShapeInfo + */ + EvaluatorTensorProduct (const AlignedVector &shape_values, + const AlignedVector &shape_gradients, + const AlignedVector &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 + void + values (const Number in [], + Number out[]) const; + + template + void + gradients (const Number in [], + Number out[]) const; + + template + 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 + template + inline + void + EvaluatorTensorProduct + ::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::value; + + for (int i2=0; i2 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 0) + { + res0 = in[0] + in[stride*(mm-1)]; + res0 *= val0; + for (int ind=1; ind 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 [-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 + template + inline + void + EvaluatorTensorProduct + ::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::value; + + for (int i2=0; i2 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 + template + inline + void + EvaluatorTensorProduct + ::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::value; + + for (int i2=0; i2 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 0) + { + res0 = in[0] + in[stride*(mm-1)]; + res0 *= val0; + for (int ind=1; ind + struct EvaluatorTensorProduct + { + static const unsigned int dofs_per_cell = Utilities::fixed_int_power::value; + static const unsigned int n_q_points = Utilities::fixed_int_power::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 &shape_values, + const AlignedVector &shape_gradients, + const AlignedVector &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 + void + values (const Number in [], + Number out[]) const + { + apply(shape_values, in, out); + } + + template + void + gradients (const Number in [], + Number out[]) const + { + apply(shape_gradients, in, out); + } + + template + void + hessians (const Number in [], + Number out[]) const + { + apply(shape_hessians, in, out); + } + + template + 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 + template + inline + void + EvaluatorTensorProduct + ::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::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; i20?mid:1], xm[mid>0?mid:1]; + for (int i=0; i 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 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 0) + { + if (type == 1) + { + r0 = shapes[n_cols*offset] * xm[0]; + for (int ind=1; ind #include -#include +#include #include