Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Provide more general factory class for FEEvaluation #10927

Merged
merged 9 commits into from Sep 22, 2020
Merged

Provide more general factory class for FEEvaluation #10927

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
dc4efca
Implement factory class for FEEvaluation and FEFaceEvaluation
kronbichler Sep 17, 2020
42e7433
Use new factory class within library
kronbichler Sep 17, 2020
5e67295
Add 3/2 over-integration to list of instantiated quadratures
kronbichler Sep 18, 2020
d4697e2
Add inverse mass matrix with degree -1
kronbichler Sep 18, 2020
67101d5
Convert a few tests to template parameter -1
kronbichler Sep 21, 2020
d96e09d
New test cases for CellwiseInverseMassMatrix<dim,-1>
kronbichler Sep 21, 2020
4f7d0ca
Remove obsolete code
kronbichler Sep 21, 2020
e6d7131
Changelog
kronbichler Sep 21, 2020
1b5806b
Add documentation
kronbichler Sep 22, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
104 changes: 93 additions & 11 deletions include/deal.II/matrix_free/evaluation_kernels.h
View file
Open in desktop
Expand Up @@ -1500,8 +1500,7 @@ namespace internal
hessians_quad,
scratch_data);
}
else if (shape_info.element_type ==
internal::MatrixFreeFunctions::tensor_general)
else
{
internal::FEEvaluationImpl<
internal::MatrixFreeFunctions::tensor_general,
Expand All @@ -1517,8 +1516,6 @@ namespace internal
hessians_quad,
scratch_data);
}
else
AssertThrow(false, ExcNotImplemented());

return false;
}
Expand Down Expand Up @@ -1649,8 +1646,7 @@ namespace internal
scratch_data,
sum_into_values_array);
}
else if (shape_info.element_type ==
internal::MatrixFreeFunctions::tensor_general)
else
{
internal::FEEvaluationImpl<
internal::MatrixFreeFunctions::tensor_general,
Expand All @@ -1666,8 +1662,6 @@ namespace internal
scratch_data,
sum_into_values_array);
}
else
AssertThrow(false, ExcNotImplemented());

return false;
}
Expand Down Expand Up @@ -3535,7 +3529,8 @@ namespace internal
false,
Number> &fe_eval,
const Number * in_array,
Number * out_array)
Number * out_array,
typename std::enable_if<fe_degree != -1>::type * = nullptr)
{
constexpr unsigned int dofs_per_component =
Utilities::pow(fe_degree + 1, dim);
Expand Down Expand Up @@ -3584,6 +3579,61 @@ namespace internal
}
return false;
}

template <int fe_degree, int = 0>
static bool
run(const unsigned int n_components,
const FEEvaluationBaseData<dim,
typename Number::value_type,
false,
Number> &fe_eval,
const Number * in_array,
Number * out_array,
typename std::enable_if<fe_degree == -1>::type * = nullptr)
{
static_assert(fe_degree == -1, "Only usable for degree -1");
const unsigned int dofs_per_component =
fe_eval.get_shape_info().dofs_per_component_on_cell;

Assert(dim >= 1 || dim <= 3, ExcNotImplemented());

internal::
EvaluatorTensorProduct<internal::evaluate_general, dim, 0, 0, Number>
evaluator(fe_eval.get_shape_info().data.front().inverse_shape_values,
AlignedVector<Number>(),
AlignedVector<Number>(),
fe_eval.get_shape_info().data.front().fe_degree + 1,
fe_eval.get_shape_info().data.front().fe_degree + 1);

for (unsigned int d = 0; d < n_components; ++d)
{
const Number *in = in_array + d * dofs_per_component;
Number * out = out_array + d * dofs_per_component;
// Need to select 'apply' method with hessian slot because values
// assume symmetries that do not exist in the inverse shapes
evaluator.template values<0, true, false>(in, out);
if (dim > 1)
evaluator.template values<1, true, false>(out, out);
if (dim > 2)
evaluator.template values<2, true, false>(out, out);
}
for (unsigned int q = 0; q < dofs_per_component; ++q)
{
const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
for (unsigned int d = 0; d < n_components; ++d)
out_array[q + d * dofs_per_component] *= inverse_JxW_q;
}
for (unsigned int d = 0; d < n_components; ++d)
{
Number *out = out_array + d * dofs_per_component;
if (dim > 2)
evaluator.template values<2, false, false>(out, out);
if (dim > 1)
evaluator.template values<1, false, false>(out, out);
evaluator.template values<0, false, false>(out, out);
}
return false;
}
};


Expand All @@ -3601,7 +3651,8 @@ namespace internal
const AlignedVector<Number> &inverse_shape,
const AlignedVector<Number> &inverse_coefficients,
const Number * in_array,
Number * out_array)
Number * out_array,
typename std::enable_if<fe_degree != -1>::type * = nullptr)
{
constexpr unsigned int dofs_per_component =
Utilities::pow(fe_degree + 1, dim);
Expand Down Expand Up @@ -3653,6 +3704,23 @@ namespace internal
}
return false;
}

/**
* Version for degree = -1
*/
template <int fe_degree, int = 0>
static bool
run(const unsigned int,
const AlignedVector<Number> &,
const AlignedVector<Number> &,
const Number *,
Number *,
typename std::enable_if<fe_degree == -1>::type * = nullptr)
{
static_assert(fe_degree == -1, "Only usable for degree -1");
Assert(false, ExcNotImplemented());
return false;
}
};


Expand All @@ -3669,7 +3737,8 @@ namespace internal
run(const unsigned int n_desired_components,
const AlignedVector<Number> &inverse_shape,
const Number * in_array,
Number * out_array)
Number * out_array,
typename std::enable_if<fe_degree != -1>::type * = nullptr)
{
constexpr unsigned int dofs_per_cell = Utilities::pow(fe_degree + 1, dim);
internal::EvaluatorTensorProduct<internal::evaluate_evenodd,
Expand Down Expand Up @@ -3702,6 +3771,19 @@ namespace internal
}
return false;
}

template <int fe_degree, int = 0>
static bool
run(const unsigned int,
const AlignedVector<Number> &,
const Number *,
Number *,
typename std::enable_if<fe_degree == -1>::type * = nullptr)
{
static_assert(fe_degree == -1, "Only usable for degree -1");
Assert(false, ExcNotImplemented());
return false;
}
};

} // end of namespace internal
Expand Down
180 changes: 180 additions & 0 deletions include/deal.II/matrix_free/evaluation_template_factory.h
View file
Open in desktop
@@ -0,0 +1,180 @@
// ---------------------------------------------------------------------
//
// Copyright (C) 2020 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------


#ifndef dealii_matrix_free_evaluation_template_factory_h
#define dealii_matrix_free_evaluation_template_factory_h


#include <deal.II/base/config.h>

#include <deal.II/matrix_free/dof_info.h>
#include <deal.II/matrix_free/evaluation_flags.h>
#include <deal.II/matrix_free/shape_info.h>


DEAL_II_NAMESPACE_OPEN


template <int, typename, bool, typename>
class FEEvaluationBaseData;


namespace internal
{
template <int dim,
typename Number,
typename VectorizedArrayType = VectorizedArray<Number>>
struct FEEvaluationFactory
{
static void
evaluate(
const unsigned int n_components,
const EvaluationFlags::EvaluationFlags evaluation_flag,
const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &shape_info,
VectorizedArrayType *values_dofs_actual,
VectorizedArrayType *values_quad,
VectorizedArrayType *gradients_quad,
VectorizedArrayType *hessians_quad,
VectorizedArrayType *scratch_data);

static void
integrate(
const unsigned int n_components,
const EvaluationFlags::EvaluationFlags integration_flag,
const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &shape_info,
VectorizedArrayType *values_dofs_actual,
VectorizedArrayType *values_quad,
VectorizedArrayType *gradients_quad,
VectorizedArrayType *scratch_data,
const bool sum_into_values_array);
};



template <int dim,
typename Number,
typename VectorizedArrayType = VectorizedArray<Number>>
struct FEFaceEvaluationFactory
{
static void
evaluate(const unsigned int n_components,
const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
const VectorizedArrayType * values_array,
VectorizedArrayType * values_quad,
VectorizedArrayType * gradients_quad,
VectorizedArrayType * scratch_data,
const bool evaluate_values,
const bool evaluate_gradients,
const unsigned int face_no,
const unsigned int subface_index,
const unsigned int face_orientation,
const Table<2, unsigned int> &orientation_map);

static void
integrate(const unsigned int n_components,
const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
VectorizedArrayType * values_array,
VectorizedArrayType * values_quad,
VectorizedArrayType * gradients_quad,
VectorizedArrayType * scratch_data,
const bool integrate_values,
const bool integrate_gradients,
const unsigned int face_no,
const unsigned int subface_index,
const unsigned int face_orientation,
const Table<2, unsigned int> &orientation_map);

template <std::size_t n_face_orientations>
static bool
gather_evaluate(
const unsigned int n_components,
const Number * src_ptr,
const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
const MatrixFreeFunctions::DoFInfo & dof_info,
VectorizedArrayType * values_quad,
VectorizedArrayType * gradients_quad,
VectorizedArrayType * scratch_data,
const bool evaluate_values,
const bool evaluate_gradients,
const unsigned int active_fe_index,
const unsigned int first_selected_component,
const std::array<unsigned int, n_face_orientations> cells,
const std::array<unsigned int, n_face_orientations> face_nos,
const unsigned int subface_index,
const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
const std::array<unsigned int, n_face_orientations> face_orientations,
const Table<2, unsigned int> & orientation_map);

template <std::size_t n_face_orientations>
static bool
integrate_scatter(
const unsigned int n_components,
Number * dst_ptr,
const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
const MatrixFreeFunctions::DoFInfo & dof_info,
VectorizedArrayType * values_array,
VectorizedArrayType * values_quad,
VectorizedArrayType * gradients_quad,
VectorizedArrayType * scratch_data,
const bool integrate_values,
const bool integrate_gradients,
const unsigned int active_fe_index,
const unsigned int first_selected_component,
const std::array<unsigned int, n_face_orientations> cells,
const std::array<unsigned int, n_face_orientations> face_nos,
const unsigned int subface_index,
const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
const std::array<unsigned int, n_face_orientations> face_orientations,
const Table<2, unsigned int> & orientation_map);
};



template <int dim,
typename Number,
typename VectorizedArrayType = VectorizedArray<Number>>
struct CellwiseInverseMassFactory
{
static void
apply(const unsigned int n_components,
const unsigned int fe_degree,
const FEEvaluationBaseData<dim, Number, false, VectorizedArrayType>
& fe_eval,
const VectorizedArrayType *in_array,
VectorizedArrayType * out_array);

static void
apply(const unsigned int n_components,
const unsigned int fe_degree,
const AlignedVector<VectorizedArrayType> &inverse_shape,
const AlignedVector<VectorizedArrayType> &inverse_coefficients,
const VectorizedArrayType * in_array,
VectorizedArrayType * out_array);

static void
transform_from_q_points_to_basis(
const unsigned int n_components,
const unsigned int fe_degree,
const AlignedVector<VectorizedArrayType> &inverse_shape,
const VectorizedArrayType * in_array,
VectorizedArrayType * out_array);
};

} // end of namespace internal

DEAL_II_NAMESPACE_CLOSE

#endif