const unsigned int n_components,
const EvaluationFlags::EvaluationFlags evaluation_flag,
typename Processor::Number2_ * global_vector_ptr,
- const std::vector<ArrayView<const typename Processor::Number_>> *sm_ptr,
- const EvaluationData & fe_eval,
- typename Processor::VectorizedArrayType_ * temp1)
+ const std::vector<ArrayView<const typename Processor::Number2_>> *sm_ptr,
+ const EvaluationData & fe_eval,
+ typename Processor::VectorizedArrayType_ * temp1)
{
constexpr int dim = Processor::dim_;
constexpr int fe_degree = Processor::fe_degree_;
dof_info
.dof_indices_contiguous_sm[dof_access_index]
[cell * n_lanes + v];
- vector_ptrs[v] = const_cast<Number *>(
+ vector_ptrs[v] = const_cast<Number2_ *>(
sm_ptr->operator[](temp.first).data() +
temp.second + index_offset);
}
dof_info
.dof_indices_contiguous_sm[dof_access_index]
[cells[v]];
- vector_ptrs[v] = const_cast<Number *>(
+ vector_ptrs[v] = const_cast<Number2_ *>(
sm_ptr->operator[](temp.first).data() +
temp.second + index_offset);
}
- template <int dim,
- typename Number,
- typename VectorizedArrayType,
- typename Number2 = Number>
+ template <int dim, typename Number2, typename VectorizedArrayType>
struct FEFaceEvaluationImplGatherEvaluateSelector
{
+ using Number = typename VectorizedArrayType::value_type;
+
template <int fe_degree, int n_q_points_1d>
static bool
run(const unsigned int n_components,
const EvaluationFlags::EvaluationFlags evaluation_flag,
const Number2 * src_ptr,
- const std::vector<ArrayView<const Number>> * sm_ptr,
+ const std::vector<ArrayView<const Number2>> * sm_ptr,
FEEvaluationData<dim, VectorizedArrayType, true> &fe_eval)
{
Assert(fe_degree > -1, ExcInternalError());
supports(
const EvaluationFlags::EvaluationFlags evaluation_flag,
const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &shape_info,
- const Number * vector_ptr,
+ const Number2 * vector_ptr,
MatrixFreeFunctions::DoFInfo::IndexStorageVariants storage)
{
const unsigned int fe_degree = shape_info.data.front().fe_degree;
hermite_grad(T0 & temp_1,
T0 & temp_2,
const T1 &src_ptr_1,
- const T2 &src_ptr_2,
+ const T1 &src_ptr_2,
const T2 &grad_weight)
{
// case 3a)
- template <int dim,
- typename Number,
- typename VectorizedArrayType,
- typename Number2 = Number>
+ template <int dim, typename Number2, typename VectorizedArrayType>
struct FEFaceEvaluationImplIntegrateScatterSelector
{
+ using Number = typename VectorizedArrayType::value_type;
+
template <int fe_degree, int n_q_points_1d>
static bool
run(const unsigned int n_components,
this->data->data.front().n_q_points_1d) &&
internal::FEFaceEvaluationImplGatherEvaluateSelector<
dim,
- Number,
- VectorizedArrayType>::supports(evaluation_flag,
- *this->data,
- internal::get_beginning<Number>(
- input_vector),
- this->dof_info->index_storage_variants
- [this->dof_access_index][this->cell]))
+ typename VectorType::value_type,
+ VectorizedArrayType>::
+ supports(evaluation_flag,
+ *this->data,
+ internal::get_beginning<typename VectorType::value_type>(
+ input_vector),
+ this->dof_info->index_storage_variants[this->dof_access_index]
+ [this->cell]))
{
if (fe_degree > -1)
{
internal::FEFaceEvaluationImplGatherEvaluateSelector<
dim,
- Number,
+ typename VectorType::value_type,
VectorizedArrayType>::template run<fe_degree,
n_q_points_1d>(
n_components,
evaluation_flag,
- internal::get_beginning<Number>(input_vector),
+ internal::get_beginning<typename VectorType::value_type>(
+ input_vector),
shared_vector_data,
*this);
}
{
internal::FEFaceEvaluationGatherFactory<
dim,
- Number,
+ typename VectorType::value_type,
VectorizedArrayType>::evaluate(n_components,
evaluation_flag,
- internal::get_beginning<Number>(
+ internal::get_beginning<
+ typename VectorType::value_type>(
input_vector),
shared_vector_data,
*this);
this->data->data.front().n_q_points_1d) &&
internal::FEFaceEvaluationImplGatherEvaluateSelector<
dim,
- Number,
- VectorizedArrayType>::supports(integration_flag,
- *this->data,
- internal::get_beginning<Number>(
- destination),
- this->dof_info->index_storage_variants
- [this->dof_access_index][this->cell]))
+ typename VectorType::value_type,
+ VectorizedArrayType>::
+ supports(integration_flag,
+ *this->data,
+ internal::get_beginning<typename VectorType::value_type>(
+ destination),
+ this->dof_info->index_storage_variants[this->dof_access_index]
+ [this->cell]))
{
if (fe_degree > -1)
{
internal::FEFaceEvaluationImplIntegrateScatterSelector<
dim,
- Number,
+ typename VectorType::value_type,
VectorizedArrayType>::template run<fe_degree,
n_q_points_1d>(
n_components,
integration_flag,
- internal::get_beginning<Number>(destination),
+ internal::get_beginning<typename VectorType::value_type>(
+ destination),
shared_vector_data,
*this);
}
{
internal::FEFaceEvaluationGatherFactory<
dim,
- Number,
+ typename VectorType::value_type,
VectorizedArrayType>::integrate(n_components,
integration_flag,
- internal::get_beginning<Number>(
+ internal::get_beginning<
+ typename VectorType::value_type>(
destination),
shared_vector_data,
*this);
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 - 2021 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// like faces_value_optimization but mixed types
+
+#include <deal.II/base/function.h>
+
+#include <deal.II/fe/fe_dgq.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/affine_constraints.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include "../tests.h"
+
+
+template <int dim, int fe_degree, typename number>
+class MatrixFreeTest
+{
+public:
+ MatrixFreeTest(const MatrixFree<dim, number> &data)
+ : data(data)
+ , error(2)
+ {}
+
+ ~MatrixFreeTest()
+ {
+ deallog << "Error between variants: " << error[0] / error[1] << std::endl;
+ }
+
+ void
+ check_error(Vector<number> &src) const
+ {
+ data.loop(&MatrixFreeTest::local_apply,
+ &MatrixFreeTest::local_apply_face,
+ &MatrixFreeTest::local_apply_face,
+ this,
+ src,
+ src);
+ }
+
+private:
+ void
+ local_apply(const MatrixFree<dim, number> &,
+ Vector<number> &,
+ const Vector<number> &,
+ const std::pair<unsigned int, unsigned int> &) const
+ {}
+
+ void
+ local_apply_face(
+ const MatrixFree<dim, number> &data,
+ Vector<number> &,
+ const Vector<number> & src_,
+ const std::pair<unsigned int, unsigned int> &face_range) const
+ {
+ FEFaceEvaluation<dim, fe_degree, fe_degree + 1, 1, number> ref(data, true);
+ FEFaceEvaluation<dim, fe_degree, fe_degree + 1, 1, number> check(data,
+ true);
+
+ Vector<float> src;
+ src = src_;
+
+ for (unsigned int face = face_range.first; face < face_range.second; ++face)
+ {
+ ref.reinit(face);
+ check.reinit(face);
+
+ ref.read_dof_values(src);
+ ref.evaluate(EvaluationFlags::values);
+ check.gather_evaluate(src, EvaluationFlags::values);
+
+ for (unsigned int q = 0; q < ref.n_q_points; ++q)
+ {
+ VectorizedArray<number> diff =
+ ref.get_value(q) - check.get_value(q);
+ for (unsigned int v = 0; v < VectorizedArray<number>::size(); ++v)
+ {
+ if (std::abs(diff[v]) > 1e-12)
+ {
+ deallog << "Error detected on interior face " << face
+ << ", v=" << v << ", q=" << q << "!" << std::endl;
+ deallog << "ref: ";
+ for (unsigned int i = 0; i < ref.dofs_per_cell; ++i)
+ deallog << ref.get_dof_value(i)[v] << " ";
+ deallog << std::endl;
+ deallog << "done: " << check.get_value(q)[v]
+ << " instead of " << ref.get_value(q)[v]
+ << std::endl;
+
+ deallog
+ << data.get_face_info(face).cells_interior[v] << " "
+ << (int)data.get_face_info(face).interior_face_no << " "
+ << (int)data.get_face_info(face).face_orientation << " "
+ << (int)data.get_face_info(face).subface_index
+ << std::endl;
+ deallog << std::endl;
+ }
+ error[0] += std::abs(diff[v]);
+ error[1] += std::abs(ref.get_value(q)[v]);
+ }
+ }
+ }
+
+ FEFaceEvaluation<dim, fe_degree, fe_degree + 1, 1, number> refr(data,
+ false);
+ FEFaceEvaluation<dim, fe_degree, fe_degree + 1, 1, number> checkr(data,
+ false);
+
+ for (unsigned int face = face_range.first;
+ face < std::min(data.n_inner_face_batches(), face_range.second);
+ face++)
+ {
+ refr.reinit(face);
+ checkr.reinit(face);
+
+ refr.read_dof_values(src);
+ refr.evaluate(EvaluationFlags::values);
+ checkr.gather_evaluate(src, EvaluationFlags::values);
+
+ for (unsigned int q = 0; q < ref.n_q_points; ++q)
+ {
+ VectorizedArray<number> diff =
+ (refr.get_value(q) - checkr.get_value(q));
+ for (unsigned int v = 0; v < VectorizedArray<number>::size(); ++v)
+ {
+ if (std::abs(diff[v]) > 1e-12)
+ {
+ deallog << "Error detected on exterior face " << face
+ << ", v=" << v << ", q=" << q << "!" << std::endl;
+ deallog << "ref: ";
+ for (unsigned int i = 0; i < ref.dofs_per_cell; ++i)
+ deallog << refr.get_dof_value(i)[v] << " ";
+ deallog << std::endl;
+ deallog << "done: " << checkr.get_value(q)[v]
+ << " instead of " << refr.get_value(q)[v]
+ << std::endl;
+
+ deallog
+ << data.get_face_info(face).cells_exterior[v] << " "
+ << (int)data.get_face_info(face).exterior_face_no << " "
+ << (int)data.get_face_info(face).face_orientation << " "
+ << (int)data.get_face_info(face).subface_index
+ << std::endl;
+ deallog << std::endl;
+ }
+ error[0] += std::abs(diff[v]);
+ error[1] += std::abs(refr.get_value(q)[v]);
+ }
+ }
+ }
+ }
+
+ const MatrixFree<dim, number> &data;
+ mutable std::vector<double> error;
+};
+
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_ball(tria);
+ tria.begin_active()->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+ typename Triangulation<dim>::active_cell_iterator cell, endc;
+ cell = tria.begin_active();
+ endc = tria.end();
+ for (; cell != endc; ++cell)
+ if (cell->center().norm() < 0.5)
+ cell->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+ tria.begin(tria.n_levels() - 1)->set_refine_flag();
+ tria.last()->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+ tria.refine_global(1);
+ cell = tria.begin_active();
+ for (unsigned int i = 0; i < 10 - 3 * dim; ++i)
+ {
+ cell = tria.begin_active();
+ endc = tria.end();
+ unsigned int counter = 0;
+ for (; cell != endc; ++cell, ++counter)
+ if (counter % (7 - i) == 0)
+ cell->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+ }
+
+ FE_DGQ<dim> fe(fe_degree);
+ DoFHandler<dim> dof(tria);
+ dof.distribute_dofs(fe);
+ AffineConstraints<double> constraints;
+ constraints.close();
+
+ using number = double;
+ using number2 = float;
+
+ MatrixFree<dim, number> mf_data;
+ {
+ const QGauss<1> quad(fe_degree + 1);
+ typename MatrixFree<dim, number>::AdditionalData data;
+ data.tasks_parallel_scheme = MatrixFree<dim, number>::AdditionalData::none;
+ data.tasks_block_size = 3;
+ data.mapping_update_flags_inner_faces =
+ (update_gradients | update_JxW_values);
+ data.mapping_update_flags_boundary_faces =
+ (update_gradients | update_JxW_values);
+
+ mf_data.reinit(dof, constraints, quad, data);
+ }
+
+ MatrixFreeTest<dim, fe_degree, number> mf(mf_data);
+
+ Vector<number> in(dof.n_dofs());
+
+ for (unsigned int i = 0; i < dof.n_dofs(); ++i)
+ {
+ if (constraints.is_constrained(i))
+ continue;
+ const double entry = Testing::rand() / (double)RAND_MAX;
+ in(i) = entry;
+ }
+
+ mf.check_error(in);
+}
+
+
+
+int
+main()
+{
+ initlog();
+
+ deallog << std::setprecision(14);
+
+ {
+ deallog.push("2d");
+ test<2, 1>();
+ test<2, 2>();
+ deallog.pop();
+ deallog.push("3d");
+ test<3, 1>();
+ test<3, 2>();
+ deallog.pop();
+ }
+}
projects / dealii.git / commitdiff