From: Martin Kronbichler Date: Thu, 17 Sep 2020 14:20:16 +0000 (+0200) Subject: Use new factory class within library X-Git-Tag: v9.3.0-rc1~1087^2~7 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=42e7433ca9a2f8c9b24e718f68f3afa9a5bc56e0;p=dealii.git Use new factory class within library --- diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index 99b0349a47..c7561253db 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -28,8 +28,7 @@ #include #include -#include -#include +#include #include DEAL_II_NAMESPACE_OPEN @@ -1302,7 +1301,7 @@ namespace internal start_indices, cell_points.data()); - SelectEvaluator::evaluate( + FEEvaluationFactory::evaluate( dim, EvaluationFlags::values | EvaluationFlags::gradients | (update_flags_cells & update_jacobian_grads ? @@ -2126,21 +2125,21 @@ namespace internal // now let the matrix-free evaluators provide us with the // data on faces - FEFaceEvaluationImplEvaluateSelector:: - template run<-1, 0>(dim, - shape_info, - cell_points.data(), - face_quads.data(), - face_grads.data(), - scratch_data.data(), - true, - true, - face_no, - GeometryInfo::max_children_per_cell, - faces[face].face_orientation > 8 ? - faces[face].face_orientation - 8 : - 0, - my_data.descriptor[0].face_orientations); + FEFaceEvaluationFactory::evaluate( + dim, + shape_info, + cell_points.data(), + face_quads.data(), + face_grads.data(), + scratch_data.data(), + true, + true, + face_no, + GeometryInfo::max_children_per_cell, + faces[face].face_orientation > 8 ? + faces[face].face_orientation - 8 : + 0, + my_data.descriptor[0].face_orientations); if (update_flags_faces & update_quadrature_points) @@ -2243,22 +2242,21 @@ namespace internal start_indices, cell_points.data()); - FEFaceEvaluationImplEvaluateSelector:: - template run<-1, 0>( - dim, - shape_info, - cell_points.data(), - face_quads.data(), - face_grads.data(), - scratch_data.data(), - false, - true, - faces[face].exterior_face_no, - faces[face].subface_index, - faces[face].face_orientation < 8 ? - faces[face].face_orientation : - 0, - my_data.descriptor[0].face_orientations); + FEFaceEvaluationFactory:: + evaluate(dim, + shape_info, + cell_points.data(), + face_quads.data(), + face_grads.data(), + scratch_data.data(), + false, + true, + faces[face].exterior_face_no, + faces[face].subface_index, + faces[face].face_orientation < 8 ? + faces[face].face_orientation : + 0, + my_data.descriptor[0].face_orientations); for (unsigned int q = 0; q < n_points_compute; ++q) { diff --git a/include/deal.II/numerics/vector_tools_project.templates.h b/include/deal.II/numerics/vector_tools_project.templates.h index 6db682e4e1..45d366dd92 100644 --- a/include/deal.II/numerics/vector_tools_project.templates.h +++ b/include/deal.II/numerics/vector_tools_project.templates.h @@ -149,13 +149,9 @@ namespace VectorTools /* * MatrixFree implementation of project() for an arbitrary number of - * components and arbitrary degree of the FiniteElement. + * components of the FiniteElement. */ - template + template void project_matrix_free( const Mapping & mapping, @@ -180,9 +176,6 @@ namespace VectorTools Assert(dof.get_fe(0).n_components() == function.n_components, ExcDimensionMismatch(dof.get_fe(0).n_components(), function.n_components)); - Assert(fe_degree == -1 || - dof.get_fe().degree == static_cast(fe_degree), - ExcDimensionMismatch(fe_degree, dof.get_fe().degree)); Assert(dof.get_fe(0).n_components() == components, ExcDimensionMismatch(components, dof.get_fe(0).n_components())); @@ -201,8 +194,8 @@ namespace VectorTools additional_data); using MatrixType = MatrixFreeOperators::MassOperator< dim, - fe_degree, - fe_degree + 2, + -1, + 0, components, LinearAlgebra::distributed::Vector>; MatrixType mass_matrix; @@ -222,8 +215,7 @@ namespace VectorTools // account for inhomogeneous constraints inhomogeneities.update_ghost_values(); - FEEvaluation phi( - *matrix_free); + FEEvaluation phi(*matrix_free); for (unsigned int cell = 0; cell < matrix_free->n_macro_cells(); ++cell) { phi.reinit(cell); @@ -255,81 +247,6 @@ namespace VectorTools - /** - * Helper interface. After figuring out the number of components in - * project_matrix_free_component, we determine the degree of the - * FiniteElement and call project_matrix_free with the appropriate - * template arguments. - */ - template - void - project_matrix_free_degree( - const Mapping & mapping, - const DoFHandler &dof, - const AffineConstraints &constraints, - const Quadrature & quadrature, - const Function< - spacedim, - typename LinearAlgebra::distributed::Vector::value_type> - & function, - LinearAlgebra::distributed::Vector &work_result, - const bool enforce_zero_boundary, - const Quadrature & q_boundary, - const bool project_to_boundary_first) - { - switch (dof.get_fe().degree) - { - case 1: - project_matrix_free(mapping, - dof, - constraints, - quadrature, - function, - work_result, - enforce_zero_boundary, - q_boundary, - project_to_boundary_first); - break; - - case 2: - project_matrix_free(mapping, - dof, - constraints, - quadrature, - function, - work_result, - enforce_zero_boundary, - q_boundary, - project_to_boundary_first); - break; - - case 3: - project_matrix_free(mapping, - dof, - constraints, - quadrature, - function, - work_result, - enforce_zero_boundary, - q_boundary, - project_to_boundary_first); - break; - - default: - project_matrix_free(mapping, - dof, - constraints, - quadrature, - function, - work_result, - enforce_zero_boundary, - q_boundary, - project_to_boundary_first); - } - } - - - // Helper interface for the matrix-free implementation of project(). // Used to determine the number of components. template @@ -351,51 +268,51 @@ namespace VectorTools switch (dof.get_fe(0).n_components()) { case 1: - project_matrix_free_degree<1>(mapping, - dof, - constraints, - quadrature, - function, - work_result, - enforce_zero_boundary, - q_boundary, - project_to_boundary_first); + project_matrix_free<1>(mapping, + dof, + constraints, + quadrature, + function, + work_result, + enforce_zero_boundary, + q_boundary, + project_to_boundary_first); break; case 2: - project_matrix_free_degree<2>(mapping, - dof, - constraints, - quadrature, - function, - work_result, - enforce_zero_boundary, - q_boundary, - project_to_boundary_first); + project_matrix_free<2>(mapping, + dof, + constraints, + quadrature, + function, + work_result, + enforce_zero_boundary, + q_boundary, + project_to_boundary_first); break; case 3: - project_matrix_free_degree<3>(mapping, - dof, - constraints, - quadrature, - function, - work_result, - enforce_zero_boundary, - q_boundary, - project_to_boundary_first); + project_matrix_free<3>(mapping, + dof, + constraints, + quadrature, + function, + work_result, + enforce_zero_boundary, + q_boundary, + project_to_boundary_first); break; case 4: - project_matrix_free_degree<4>(mapping, - dof, - constraints, - quadrature, - function, - work_result, - enforce_zero_boundary, - q_boundary, - project_to_boundary_first); + project_matrix_free<4>(mapping, + dof, + constraints, + quadrature, + function, + work_result, + enforce_zero_boundary, + q_boundary, + project_to_boundary_first); break; default: @@ -581,7 +498,7 @@ namespace VectorTools vec_result); } - template + template void project_parallel( const Mapping & mapping, @@ -598,9 +515,6 @@ namespace VectorTools ExcDimensionMismatch(dof.get_fe(0).n_components(), 1)); Assert(vec_result.size() == dof.n_dofs(), ExcDimensionMismatch(vec_result.size(), dof.n_dofs())); - Assert(fe_degree == -1 || - dof.get_fe().degree == static_cast(fe_degree), - ExcDimensionMismatch(fe_degree, dof.get_fe().degree)); // set up mass matrix and right hand side typename MatrixFree::AdditionalData additional_data; @@ -615,12 +529,8 @@ namespace VectorTools constraints, QGauss<1>(dof.get_fe().degree + 2), additional_data); - using MatrixType = MatrixFreeOperators::MassOperator< - dim, - fe_degree, - fe_degree + 2, - 1, - LinearAlgebra::distributed::Vector>; + using MatrixType = MatrixFreeOperators:: + MassOperator>; MatrixType mass_matrix; mass_matrix.initialize(matrix_free); mass_matrix.compute_diagonal(); @@ -696,11 +606,7 @@ namespace VectorTools - template + template void project_parallel( std::shared_ptr> @@ -720,16 +626,9 @@ namespace VectorTools ExcDimensionMismatch(dof.get_fe(0).n_components(), 1)); Assert(vec_result.size() == dof.n_dofs(), ExcDimensionMismatch(vec_result.size(), dof.n_dofs())); - Assert(fe_degree == -1 || - dof.get_fe().degree == static_cast(fe_degree), - ExcDimensionMismatch(fe_degree, dof.get_fe().degree)); - using MatrixType = MatrixFreeOperators::MassOperator< - dim, - fe_degree, - n_q_points_1d, - 1, - LinearAlgebra::distributed::Vector>; + using MatrixType = MatrixFreeOperators:: + MassOperator>; MatrixType mass_matrix; mass_matrix.initialize(matrix_free, {fe_component}); mass_matrix.compute_diagonal(); @@ -744,9 +643,8 @@ namespace VectorTools // assemble right hand side: { - FEEvaluation fe_eval( - *matrix_free, fe_component); - const unsigned int n_cells = matrix_free->n_macro_cells(); + FEEvaluation fe_eval(*matrix_free, fe_component); + const unsigned int n_cells = matrix_free->n_cell_batches(); const unsigned int n_q_points = fe_eval.n_q_points; for (unsigned int cell = 0; cell < n_cells; ++cell) @@ -857,24 +755,8 @@ namespace VectorTools const unsigned int)> & func, VectorType & vec_result) { - switch (dof.get_fe().degree) - { - case 1: - internal::project_parallel( - mapping, dof, constraints, quadrature, func, vec_result); - break; - case 2: - internal::project_parallel( - mapping, dof, constraints, quadrature, func, vec_result); - break; - case 3: - internal::project_parallel( - mapping, dof, constraints, quadrature, func, vec_result); - break; - default: - internal::project_parallel( - mapping, dof, constraints, quadrature, func, vec_result); - } + internal::project_parallel( + mapping, dof, constraints, quadrature, func, vec_result); } @@ -886,38 +768,15 @@ namespace VectorTools typename VectorType::value_type, VectorizedArray>> matrix_free, const AffineConstraints &constraints, - const unsigned int n_q_points_1d, + const unsigned int, const std::function( const unsigned int, const unsigned int)> &func, VectorType & vec_result, const unsigned int fe_component) { - const unsigned int fe_degree = - matrix_free->get_dof_handler(fe_component).get_fe().degree; - - if (fe_degree + 1 == n_q_points_1d) - switch (fe_degree) - { - case 1: - internal::project_parallel( - matrix_free, constraints, func, vec_result, fe_component); - break; - case 2: - internal::project_parallel( - matrix_free, constraints, func, vec_result, fe_component); - break; - case 3: - internal::project_parallel( - matrix_free, constraints, func, vec_result, fe_component); - break; - default: - internal::project_parallel( - matrix_free, constraints, func, vec_result, fe_component); - } - else - internal::project_parallel( - matrix_free, constraints, func, vec_result, fe_component); + internal::project_parallel( + matrix_free, constraints, func, vec_result, fe_component); } diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index fa08b3466a..3d878bbc04 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -36,8 +36,7 @@ #include #include -#include -#include +#include #include #include @@ -1541,15 +1540,17 @@ namespace internal } // do the actual tensorized evaluation - SelectEvaluator>::evaluate( - n_comp, - evaluation_flag, - data.shape_info, - data.values_dofs.begin(), - data.values_quad.begin(), - data.gradients_quad.begin(), - data.hessians_quad.begin(), - data.scratch.begin()); + internal::FEEvaluationFactory< + dim, + double, + VectorizedArray>::evaluate(n_comp, + evaluation_flag, + data.shape_info, + data.values_dofs.begin(), + data.values_quad.begin(), + data.gradients_quad.begin(), + data.hessians_quad.begin(), + data.scratch.begin()); } // do the postprocessing