From: Peter Munch Date: Wed, 8 Dec 2021 21:22:03 +0000 (+0100) Subject: Use FEEvaluationHangingNodes in MGTwoLevelTransfer X-Git-Tag: v9.4.0-rc1~631^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F13048%2Fhead;p=dealii.git Use FEEvaluationHangingNodes in MGTwoLevelTransfer --- diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.h index 35767cde1a..4b323f6a88 100644 --- a/include/deal.II/multigrid/mg_transfer_global_coarsening.h +++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.h @@ -24,6 +24,9 @@ #include #include +#include +#include + #include DEAL_II_NAMESPACE_OPEN @@ -358,6 +361,13 @@ private: * 1D restriction matrix for tensor-product elements. */ AlignedVector> restriction_matrix_1d; + + /** + * ShapeInfo description of the coarse cell. Needed during the + * fast application of hanging-node constraints. + */ + internal::MatrixFreeFunctions::ShapeInfo> + shape_info_coarse; }; /** @@ -432,6 +442,12 @@ private: */ std::vector level_dof_indices_coarse; + /** + * DoF indices of the coarse cells, expressed in indices local to the MPI + * rank. + */ + std::vector level_dof_indices_coarse_plain; + /** * DoF indices of the fine cells, expressed in indices local to the MPI * rank. @@ -443,7 +459,28 @@ private: */ unsigned int n_components; + /** + * Refinement configuration of the coarse cells, needed for fast + * application of hanging-node constraints. If no hanging-node + * constraints have to be applied or the fast algorithm is not + * applicable, the vector is empty. + */ + std::vector + coarse_cell_refinement_configurations; + friend class internal::MGTwoLevelTransferImplementation; + + /** + * Apply hanging-node constrains with fast algorithm. + */ + void + apply_hanging_node_constraints( + const MGTransferScheme &scheme, + const internal::MatrixFreeFunctions::ConstraintKinds + * coarse_cell_refinement_configurations_ptr, + const unsigned int n_lanes_filled, + const bool transpose, + AlignedVector> &evaluation_data_coarse) const; }; diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h index b3dfc48ecc..cebc27b6ad 100644 --- a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h +++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h @@ -40,6 +40,7 @@ #include #include +#include #include #include @@ -325,7 +326,7 @@ namespace internal const unsigned int fe_shift_1d, const unsigned int fe_degree) { - std::vector> cell_local_chilren_indices( + std::vector> cell_local_children_indices( GeometryInfo::max_children_per_cell, std::vector(n_dofs_per_cell_coarse)); for (unsigned int c = 0; c < GeometryInfo::max_children_per_cell; @@ -333,23 +334,23 @@ namespace internal get_child_offset(c, fe_shift_1d, fe_degree, - cell_local_chilren_indices[c]); - return cell_local_chilren_indices; + cell_local_children_indices[c]); + return cell_local_children_indices; } template std::vector> get_child_offsets_general(const unsigned int n_dofs_per_cell_coarse) { - std::vector> cell_local_chilren_indices( + std::vector> cell_local_children_indices( GeometryInfo::max_children_per_cell, std::vector(n_dofs_per_cell_coarse)); for (unsigned int c = 0, k = 0; c < GeometryInfo::max_children_per_cell; c++) for (unsigned int d = 0; d < n_dofs_per_cell_coarse; ++d, ++k) - cell_local_chilren_indices[c][d] = k; - return cell_local_chilren_indices; + cell_local_children_indices[c][d] = k; + return cell_local_children_indices; } template @@ -1267,9 +1268,16 @@ namespace internal AssertDimension(min_active_fe_indices[0], max_active_fe_indices[0]); AssertDimension(min_active_fe_indices[1], max_active_fe_indices[1]); - const auto &fe_fine = dof_handler_fine.get_fe(min_active_fe_indices[0]); - const auto &fe_coarse = - dof_handler_coarse.get_fe(min_active_fe_indices[1]); + // set up two mg-schemes + // (0) no refinement -> identity + // (1) h-refinement + transfer.schemes.resize(2); + + const unsigned int fe_index_fine = min_active_fe_indices[0]; + const unsigned int fe_index_coarse = min_active_fe_indices[1]; + + const auto &fe_fine = dof_handler_fine.get_fe(fe_index_fine); + const auto &fe_coarse = dof_handler_coarse.get_fe(fe_index_coarse); // extract number of components AssertDimension(fe_fine.n_components(), fe_coarse.n_components()); @@ -1278,6 +1286,110 @@ namespace internal const auto reference_cell = dof_handler_fine.get_fe(0).reference_cell(); + std::unique_ptr> + hanging_nodes; + std::vector> lexicographic_mappings; + + { + // algorithm is only needed on active levels + bool use_fast_hanging_node_algorithm = + mg_level_coarse == numbers::invalid_unsigned_int; + + const auto &tria = dof_handler_coarse.get_triangulation(); + + // algorithm can be only used on meshes consisting of hypercube and + // simplices + if (use_fast_hanging_node_algorithm) + { + const auto &reference_cells = tria.get_reference_cells(); + use_fast_hanging_node_algorithm = + std::all_of(reference_cells.begin(), + reference_cells.end(), + [](const auto &r) { + return r.is_hyper_cube() || r.is_simplex(); + }); + } + + // local p-refinement is not supported + if (use_fast_hanging_node_algorithm) + { + const auto &fes = dof_handler_coarse.get_fe_collection(); + + use_fast_hanging_node_algorithm &= + std::all_of(fes.begin(), fes.end(), [&fes](const auto &fe) { + return fes[0].compare_for_domination(fe) == + FiniteElementDomination::Domination:: + either_element_can_dominate; + }); + } + + // check that all components are either supported or not + if (use_fast_hanging_node_algorithm) + { + const std::vector> supported_components = + internal::MatrixFreeFunctions::HangingNodes< + dim>::compute_supported_components(dof_handler_coarse + .get_fe_collection()); + + use_fast_hanging_node_algorithm &= std::any_of( + supported_components.begin(), + supported_components.end(), + [](const auto &supported_components_per_fe) { + return std::all_of(supported_components_per_fe.begin(), + supported_components_per_fe.end(), + [](const auto &a) { return a == true; }); + }); + + use_fast_hanging_node_algorithm &= + std::all_of(supported_components.begin(), + supported_components.end(), + [](const auto &supported_components_per_fe) { + return std::all_of( + supported_components_per_fe.begin(), + supported_components_per_fe.end(), + [&supported_components_per_fe](const auto &a) { + return a == supported_components_per_fe[0]; + }); + }); + } + + use_fast_hanging_node_algorithm = false; + + // create helper class + if (use_fast_hanging_node_algorithm) + { + hanging_nodes = std::make_unique< + internal::MatrixFreeFunctions::HangingNodes>(tria); + + const auto &fes = dof_handler_coarse.get_fe_collection(); + lexicographic_mappings.resize(fes.size()); + + for (unsigned int i = 0; i < fes.size(); ++i) + if (fes[i].reference_cell().is_hyper_cube()) + { + const Quadrature<1> dummy_quadrature( + std::vector>(1, Point<1>())); + internal::MatrixFreeFunctions::ShapeInfo shape_info; + shape_info.reinit(dummy_quadrature, fes[i], 0); + lexicographic_mappings[i] = + shape_info.lexicographic_numbering; + } + + if (fes[fe_index_coarse].reference_cell().is_hyper_cube()) + { + const Quadrature<1> dummy_quadrature( + std::vector>(1, Point<1>())); + internal::MatrixFreeFunctions::ShapeInfo< + VectorizedArray> + shape_info; + shape_info.reinit(dummy_quadrature, fes[fe_index_coarse], 0); + + transfer.schemes[0].shape_info_coarse = + transfer.schemes[1].shape_info_coarse = shape_info; + } + } + } + // create partitioners and vectors for internal purposes { // ... for fine mesh @@ -1352,11 +1464,6 @@ namespace internal }); }; - // set up two mg-schemes - // (0) no refinement -> identity - // (1) h-refinement - transfer.schemes.resize(2); - // check if FE is the same AssertDimension(fe_coarse.n_dofs_per_cell(), fe_fine.n_dofs_per_cell()); @@ -1402,7 +1509,7 @@ namespace internal } - const auto cell_local_chilren_indices = + const auto cell_local_children_indices = (reference_cell == ReferenceCells::get_hypercube()) ? get_child_offsets(transfer.schemes[0].n_dofs_per_cell_coarse, is_feq ? fe_fine.degree : (fe_fine.degree + 1), @@ -1432,7 +1539,10 @@ namespace internal // indices { transfer.level_dof_indices_fine.resize(n_dof_indices_fine.back()); - transfer.level_dof_indices_coarse.resize(n_dof_indices_coarse.back()); + if (hanging_nodes) + transfer.level_dof_indices_coarse.resize(n_dof_indices_coarse.back()); + transfer.level_dof_indices_coarse_plain.resize( + n_dof_indices_coarse.back()); std::vector local_dof_indices( transfer.schemes[0].n_dofs_per_cell_coarse); @@ -1464,6 +1574,8 @@ namespace internal // ------------------------------ indices ------------------------------ unsigned int *level_dof_indices_coarse_0 = transfer.level_dof_indices_coarse.data(); + unsigned int *level_dof_indices_coarse_plain_0 = + transfer.level_dof_indices_coarse_plain.data(); unsigned int *level_dof_indices_fine_0 = transfer.level_dof_indices_fine.data(); @@ -1471,10 +1583,37 @@ namespace internal level_dof_indices_coarse_0 + transfer.schemes[0].n_dofs_per_cell_coarse * transfer.schemes[0].n_coarse_cells; + unsigned int *level_dof_indices_coarse_plain_1 = + level_dof_indices_coarse_plain_0 + + transfer.schemes[0].n_dofs_per_cell_coarse * + transfer.schemes[0].n_coarse_cells; unsigned int *level_dof_indices_fine_1 = level_dof_indices_fine_0 + transfer.schemes[0].n_dofs_per_cell_fine * transfer.schemes[0].n_coarse_cells; + if (hanging_nodes) + transfer.coarse_cell_refinement_configurations.resize( + transfer.schemes[0].n_coarse_cells + + transfer.schemes[1].n_coarse_cells); + + internal::MatrixFreeFunctions::ConstraintKinds + *coarse_cell_refinement_configurations_0 = + transfer.coarse_cell_refinement_configurations.data(); + internal::MatrixFreeFunctions::ConstraintKinds + *coarse_cell_refinement_configurations_1 = + coarse_cell_refinement_configurations_0 + + transfer.schemes[0].n_coarse_cells; + + std::vector local_dof_indices_lex( + local_dof_indices.size()); + + const auto get_lexicographic_mapping = [&](const unsigned int c) { + if (lexicographic_mappings.size() > 0) + return lexicographic_mappings[fe_index_coarse][c]; + else + return lexicographic_numbering_coarse[c]; + }; + process_cells( [&](const auto &cell_coarse, const auto &cell_fine) { // parent @@ -1484,12 +1623,39 @@ namespace internal else cell_coarse->get_mg_dof_indices(local_dof_indices); + for (unsigned int i = 0; i < local_dof_indices.size(); ++i) + local_dof_indices_lex[i] = + local_dof_indices[get_lexicographic_mapping(i)]; + + for (auto &i : local_dof_indices_lex) + i = transfer.partitioner_coarse->global_to_local(i); + + if (hanging_nodes) + { + std::vector + mask(transfer.n_components); + + local_dof_indices = local_dof_indices_lex; + + hanging_nodes->setup_constraints(cell_coarse, + transfer.partitioner_coarse, + lexicographic_mappings, + local_dof_indices, + mask); + + for (unsigned int i = 0; + i < transfer.schemes[0].n_dofs_per_cell_coarse; + i++) + level_dof_indices_coarse_0[i] = local_dof_indices[i]; + + coarse_cell_refinement_configurations_0[0] = mask[0]; + coarse_cell_refinement_configurations_0 += 1; + } + for (unsigned int i = 0; i < transfer.schemes[0].n_dofs_per_cell_coarse; i++) - level_dof_indices_coarse_0[i] = - transfer.partitioner_coarse->global_to_local( - local_dof_indices[lexicographic_numbering_coarse[i]]); + level_dof_indices_coarse_plain_0[i] = local_dof_indices_lex[i]; } // child @@ -1507,6 +1673,8 @@ namespace internal { level_dof_indices_coarse_0 += transfer.schemes[0].n_dofs_per_cell_coarse; + level_dof_indices_coarse_plain_0 += + transfer.schemes[0].n_dofs_per_cell_coarse; level_dof_indices_fine_0 += transfer.schemes[0].n_dofs_per_cell_fine; } @@ -1520,12 +1688,41 @@ namespace internal else cell_coarse->get_mg_dof_indices(local_dof_indices); + for (unsigned int i = 0; i < local_dof_indices.size(); ++i) + local_dof_indices_lex[i] = + local_dof_indices[get_lexicographic_mapping(i)]; + + for (auto &i : local_dof_indices_lex) + i = transfer.partitioner_coarse->global_to_local(i); + + if (hanging_nodes) + { + std::vector + mask(transfer.n_components); + + local_dof_indices = local_dof_indices_lex; + + hanging_nodes->setup_constraints( + cell_coarse, + transfer.partitioner_coarse, + lexicographic_mappings, + local_dof_indices, + mask); + + for (unsigned int i = 0; + i < transfer.schemes[1].n_dofs_per_cell_coarse; + ++i) + level_dof_indices_coarse_1[i] = local_dof_indices[i]; + + coarse_cell_refinement_configurations_1[0] = mask[0]; + coarse_cell_refinement_configurations_1 += 1; + } + for (unsigned int i = 0; i < transfer.schemes[1].n_dofs_per_cell_coarse; - i++) - level_dof_indices_coarse_1[i] = - transfer.partitioner_coarse->global_to_local( - local_dof_indices[lexicographic_numbering_coarse[i]]); + ++i) + level_dof_indices_coarse_plain_1[i] = + local_dof_indices_lex[i]; } // child @@ -1533,8 +1730,8 @@ namespace internal cell_fine.get_dof_indices(local_dof_indices); for (unsigned int i = 0; i < transfer.schemes[1].n_dofs_per_cell_coarse; - i++) - level_dof_indices_fine_1[cell_local_chilren_indices[c][i]] = + ++i) + level_dof_indices_fine_1[cell_local_children_indices[c][i]] = transfer.partitioner_fine->global_to_local( local_dof_indices[lexicographic_numbering_fine[i]]); } @@ -1544,12 +1741,24 @@ namespace internal { level_dof_indices_coarse_1 += transfer.schemes[1].n_dofs_per_cell_coarse; + level_dof_indices_coarse_plain_1 += + transfer.schemes[1].n_dofs_per_cell_coarse; level_dof_indices_fine_1 += transfer.schemes[1].n_dofs_per_cell_fine; } }); + + if (hanging_nodes && + std::all_of(transfer.coarse_cell_refinement_configurations.begin(), + transfer.coarse_cell_refinement_configurations.end(), + [](const auto i) { + return i == internal::MatrixFreeFunctions:: + ConstraintKinds::unconstrained; + })) + transfer.coarse_cell_refinement_configurations.clear(); } + // ------------- prolongation matrix (0) -> identity matrix -------------- // nothing to do since for identity prolongation matrices a short-cut @@ -1685,7 +1894,7 @@ namespace internal [&](const auto &, const auto &) { for (unsigned int i = 0; i < transfer.schemes[0].n_dofs_per_cell_fine; - i++) + ++i) weights_0[i] = weight_vector.local_element(dof_indices_fine_0[i]); @@ -1695,10 +1904,10 @@ namespace internal [&](const auto &, const auto &, const auto c) { for (unsigned int i = 0; i < transfer.schemes[1].n_dofs_per_cell_coarse; - i++) - weights_1[cell_local_chilren_indices[c][i]] = + ++i) + weights_1[cell_local_children_indices[c][i]] = weight_vector.local_element( - dof_indices_fine_1[cell_local_chilren_indices[c][i]]); + dof_indices_fine_1[cell_local_children_indices[c][i]]); // move pointers (only once at the end) if (c + 1 == GeometryInfo::max_children_per_cell) @@ -1967,7 +2176,8 @@ namespace internal } transfer.level_dof_indices_fine.resize(n_dof_indices_fine.back()); - transfer.level_dof_indices_coarse.resize(n_dof_indices_coarse.back()); + transfer.level_dof_indices_coarse_plain.resize( + n_dof_indices_coarse.back()); // ------------------------------ indices ----------------------------- std::vector level_dof_indices_coarse_( @@ -1980,7 +2190,7 @@ namespace internal level_dof_indices_fine_[i] = transfer.level_dof_indices_fine.data() + n_dof_indices_fine[i]; level_dof_indices_coarse_[i] = - transfer.level_dof_indices_coarse.data() + + transfer.level_dof_indices_coarse_plain.data() + n_dof_indices_coarse[i]; } @@ -2489,9 +2699,13 @@ MGTwoLevelTransfer>:: AlignedVector evaluation_data_fine; AlignedVector evaluation_data_coarse; - const unsigned int *indices_coarse = level_dof_indices_coarse.data(); - const unsigned int *indices_fine = level_dof_indices_fine.data(); - const Number * weights = nullptr; + const unsigned int *indices_coarse = level_dof_indices_coarse.size() > 0 ? + level_dof_indices_coarse.data() : + level_dof_indices_coarse_plain.data(); + const unsigned int *indices_coarse_plain = + level_dof_indices_coarse_plain.data(); + const unsigned int *indices_fine = level_dof_indices_fine.data(); + const Number * weights = nullptr; const std::array, Utilities::pow(3, dim)> *weights_compressed = nullptr; @@ -2501,56 +2715,16 @@ MGTwoLevelTransfer>:: weights_compressed = this->weights_compressed.data(); } + const internal::MatrixFreeFunctions::ConstraintKinds + *coarse_cell_refinement_configurations_ptr = + coarse_cell_refinement_configurations.data(); + for (const auto &scheme : schemes) { - // identity -> take short cut and work directly on global vectors - if (scheme.prolongation_matrix.size() == 0 && - scheme.prolongation_matrix_1d.size() == 0) - { - for (unsigned int cell = 0; cell < scheme.n_coarse_cells; - cell += n_lanes) - { - const unsigned int n_lanes_filled = - (cell + n_lanes > scheme.n_coarse_cells) ? - (scheme.n_coarse_cells - cell) : - n_lanes; - - // read from source vector - for (unsigned int v = 0; v < n_lanes_filled; ++v) - { - if ((scheme.n_dofs_per_cell_fine != 0) && - (scheme.n_dofs_per_cell_coarse != 0)) - { - if (fine_element_is_continuous) - for (unsigned int i = 0; - i < scheme.n_dofs_per_cell_fine; - ++i) - vec_fine_ptr->local_element(indices_fine[i]) += - read_dof_values(indices_coarse[i], vec_coarse) * - weights[i]; - else - for (unsigned int i = 0; - i < scheme.n_dofs_per_cell_fine; - ++i) - vec_fine_ptr->local_element(indices_fine[i]) += - read_dof_values(indices_coarse[i], vec_coarse); - } + const bool needs_interpolation = + (scheme.prolongation_matrix.size() == 0 && + scheme.prolongation_matrix_1d.size() == 0) == false; - indices_fine += scheme.n_dofs_per_cell_fine; - indices_coarse += scheme.n_dofs_per_cell_coarse; - - if (fine_element_is_continuous) - weights += scheme.n_dofs_per_cell_fine; - } - - if (fine_element_is_continuous) - weights_compressed += 1; - } - - continue; - } - - // general case -> local restriction is needed evaluation_data_fine.resize(scheme.n_dofs_per_cell_fine); evaluation_data_coarse.resize(scheme.n_dofs_per_cell_fine); @@ -2569,30 +2743,45 @@ MGTwoLevelTransfer>:: (scheme.n_coarse_cells - cell) : n_lanes; - // read from source vector + // read from source vector, apply general-purpose constraints, and ... for (unsigned int v = 0; v < n_lanes_filled; ++v) { for (unsigned int i = 0; i < scheme.n_dofs_per_cell_coarse; ++i) evaluation_data_coarse[i][v] = read_dof_values(indices_coarse[i], this->vec_coarse); indices_coarse += scheme.n_dofs_per_cell_coarse; + indices_coarse_plain += scheme.n_dofs_per_cell_coarse; } + // ... fast hanging-node-constraints algorithm + apply_hanging_node_constraints( + scheme, + coarse_cell_refinement_configurations_ptr, + n_lanes_filled, + false, + evaluation_data_coarse); + + if (coarse_cell_refinement_configurations.size() > 0) + coarse_cell_refinement_configurations_ptr += n_lanes_filled; + // ---------------------------- coarse ------------------------------- - for (int c = n_components - 1; c >= 0; --c) - { - CellProlongator cell_prolongator( - scheme.prolongation_matrix, - scheme.prolongation_matrix_1d, - evaluation_data_coarse.begin() + c * n_scalar_dofs_coarse, - evaluation_data_fine.begin() + c * n_scalar_dofs_fine); - - if (scheme.prolongation_matrix_1d.size() > 0) - cell_transfer.run(cell_prolongator); - else - cell_prolongator.run_full(n_scalar_dofs_fine, - n_scalar_dofs_coarse); - } + if (needs_interpolation) + for (int c = n_components - 1; c >= 0; --c) + { + CellProlongator cell_prolongator( + scheme.prolongation_matrix, + scheme.prolongation_matrix_1d, + evaluation_data_coarse.begin() + c * n_scalar_dofs_coarse, + evaluation_data_fine.begin() + c * n_scalar_dofs_fine); + + if (scheme.prolongation_matrix_1d.size() > 0) + cell_transfer.run(cell_prolongator); + else + cell_prolongator.run_full(n_scalar_dofs_fine, + n_scalar_dofs_coarse); + } + else + evaluation_data_fine = evaluation_data_coarse; // TODO // ------------------------------ fine ------------------------------- // weight and write into dst vector @@ -2679,9 +2868,13 @@ MGTwoLevelTransfer>:: value * distribute_local_to_global_values[j]; }; - const unsigned int *indices_coarse = level_dof_indices_coarse.data(); - const unsigned int *indices_fine = level_dof_indices_fine.data(); - const Number * weights = nullptr; + const unsigned int *indices_coarse = level_dof_indices_coarse.size() > 0 ? + level_dof_indices_coarse.data() : + level_dof_indices_coarse_plain.data(); + const unsigned int *indices_coarse_plain = + level_dof_indices_coarse_plain.data(); + const unsigned int *indices_fine = level_dof_indices_fine.data(); + const Number * weights = nullptr; const std::array, Utilities::pow(3, dim)> *weights_compressed = nullptr; @@ -2691,59 +2884,16 @@ MGTwoLevelTransfer>:: weights_compressed = this->weights_compressed.data(); } + const internal::MatrixFreeFunctions::ConstraintKinds + *coarse_cell_refinement_configurations_ptr = + coarse_cell_refinement_configurations.data(); + for (const auto &scheme : schemes) { - // identity -> take short cut and work directly on global vectors - if (scheme.prolongation_matrix.size() == 0 && - scheme.prolongation_matrix_1d.size() == 0) - { - for (unsigned int cell = 0; cell < scheme.n_coarse_cells; - cell += n_lanes) - { - const unsigned int n_lanes_filled = - (cell + n_lanes > scheme.n_coarse_cells) ? - (scheme.n_coarse_cells - cell) : - n_lanes; - - for (unsigned int v = 0; v < n_lanes_filled; ++v) - { - if ((scheme.n_dofs_per_cell_fine != 0) && - (scheme.n_dofs_per_cell_coarse != 0)) - { - if (fine_element_is_continuous) - for (unsigned int i = 0; - i < scheme.n_dofs_per_cell_fine; - ++i) - distribute_local_to_global( - indices_coarse[i], - vec_fine_ptr->local_element(indices_fine[i]) * - weights[i], - this->vec_coarse); - else - for (unsigned int i = 0; - i < scheme.n_dofs_per_cell_fine; - ++i) - distribute_local_to_global( - indices_coarse[i], - vec_fine_ptr->local_element(indices_fine[i]), - this->vec_coarse); - } - - indices_fine += scheme.n_dofs_per_cell_fine; - indices_coarse += scheme.n_dofs_per_cell_coarse; - - if (fine_element_is_continuous) - weights += scheme.n_dofs_per_cell_fine; - } + const bool needs_interpolation = + (scheme.prolongation_matrix.size() == 0 && + scheme.prolongation_matrix_1d.size() == 0) == false; - if (fine_element_is_continuous) - weights_compressed += 1; - } - - continue; - } - - // general case -> local restriction is needed evaluation_data_fine.resize(scheme.n_dofs_per_cell_fine); evaluation_data_coarse.resize(scheme.n_dofs_per_cell_fine); @@ -2792,23 +2942,38 @@ MGTwoLevelTransfer>:: } // ------------------------------ fine ------------------------------- - for (int c = n_components - 1; c >= 0; --c) - { - CellRestrictor cell_restrictor( - scheme.prolongation_matrix, - scheme.prolongation_matrix_1d, - evaluation_data_fine.begin() + c * n_scalar_dofs_fine, - evaluation_data_coarse.begin() + c * n_scalar_dofs_coarse); - - if (scheme.prolongation_matrix_1d.size() > 0) - cell_transfer.run(cell_restrictor); - else - cell_restrictor.run_full(n_scalar_dofs_fine, - n_scalar_dofs_coarse); - } + if (needs_interpolation) + for (int c = n_components - 1; c >= 0; --c) + { + CellRestrictor cell_restrictor( + scheme.prolongation_matrix, + scheme.prolongation_matrix_1d, + evaluation_data_fine.begin() + c * n_scalar_dofs_fine, + evaluation_data_coarse.begin() + c * n_scalar_dofs_coarse); + + if (scheme.prolongation_matrix_1d.size() > 0) + cell_transfer.run(cell_restrictor); + else + cell_restrictor.run_full(n_scalar_dofs_fine, + n_scalar_dofs_coarse); + } + else + evaluation_data_coarse = evaluation_data_fine; // TODO // ----------------------------- coarse ------------------------------ - // write into dst vector + // apply fast hanging-node-constraints algorithm, ... + apply_hanging_node_constraints( + scheme, + coarse_cell_refinement_configurations_ptr, + n_lanes_filled, + true, + evaluation_data_coarse); + + + if (coarse_cell_refinement_configurations.size() > 0) + coarse_cell_refinement_configurations_ptr += n_lanes_filled; + + // ... apply general-purpose constraints, and write into dst vector for (unsigned int v = 0; v < n_lanes_filled; ++v) { for (unsigned int i = 0; i < scheme.n_dofs_per_cell_coarse; ++i) @@ -2816,6 +2981,7 @@ MGTwoLevelTransfer>:: evaluation_data_coarse[i][v], this->vec_coarse); indices_coarse += scheme.n_dofs_per_cell_coarse; + indices_coarse_plain += scheme.n_dofs_per_cell_coarse; } } } @@ -2857,8 +3023,9 @@ MGTwoLevelTransfer>:: AlignedVector evaluation_data_fine; AlignedVector evaluation_data_coarse; - const unsigned int *indices_coarse = level_dof_indices_coarse.data(); - const unsigned int *indices_fine = level_dof_indices_fine.data(); + const unsigned int *indices_coarse_plain = + level_dof_indices_coarse_plain.data(); + const unsigned int *indices_fine = level_dof_indices_fine.data(); for (const auto &scheme : schemes) { @@ -2871,11 +3038,11 @@ MGTwoLevelTransfer>:: if ((scheme.n_dofs_per_cell_fine != 0) && (scheme.n_dofs_per_cell_coarse != 0)) for (unsigned int i = 0; i < scheme.n_dofs_per_cell_fine; ++i) - this->vec_coarse.local_element(indices_coarse[i]) = + this->vec_coarse.local_element(indices_coarse_plain[i]) = vec_fine_ptr->local_element(indices_fine[i]); indices_fine += scheme.n_dofs_per_cell_fine; - indices_coarse += scheme.n_dofs_per_cell_coarse; + indices_coarse_plain += scheme.n_dofs_per_cell_coarse; } continue; @@ -2931,9 +3098,9 @@ MGTwoLevelTransfer>:: for (unsigned int v = 0; v < n_lanes_filled; ++v) { for (unsigned int i = 0; i < scheme.n_dofs_per_cell_coarse; ++i) - this->vec_coarse.local_element(indices_coarse[i]) = + this->vec_coarse.local_element(indices_coarse_plain[i]) = evaluation_data_coarse[i][v]; - indices_coarse += scheme.n_dofs_per_cell_coarse; + indices_coarse_plain += scheme.n_dofs_per_cell_coarse; } } } @@ -2949,6 +3116,58 @@ MGTwoLevelTransfer>:: +template +void +MGTwoLevelTransfer>:: + apply_hanging_node_constraints( + const MGTransferScheme &scheme, + const internal::MatrixFreeFunctions::ConstraintKinds + * coarse_cell_refinement_configurations_ptr, + const unsigned int n_lanes_filled, + const bool transpose, + AlignedVector> &evaluation_data_coarse) const +{ + if (coarse_cell_refinement_configurations.size() == 0) + return; + + using VectorizedArrayType = VectorizedArray; + + std::array + constraint_mask; + + bool hn_available = false; + + for (unsigned int v = 0; v < n_lanes_filled; ++v) + { + const auto mask = coarse_cell_refinement_configurations_ptr[v]; + + constraint_mask[v] = mask; + + hn_available |= + (mask != internal::MatrixFreeFunctions::ConstraintKinds::unconstrained); + } + + if (hn_available == true) + { + for (unsigned int v = n_lanes_filled; v < VectorizedArrayType::size(); + ++v) + constraint_mask[v] = + internal::MatrixFreeFunctions::ConstraintKinds::unconstrained; + + internal::FEEvaluationHangingNodesFactory< + dim, + Number, + VectorizedArrayType>::apply(n_components, + scheme.degree_coarse, + scheme.shape_info_coarse, + transpose, + constraint_mask, + evaluation_data_coarse.begin()); + } +} + + template void MGTwoLevelTransfer>:: @@ -3125,6 +3344,7 @@ MGTwoLevelTransfer>:: size += MemoryConsumption::memory_consumption(distribute_local_to_global_ptr); size += MemoryConsumption::memory_consumption(weights); size += MemoryConsumption::memory_consumption(level_dof_indices_coarse); + size += MemoryConsumption::memory_consumption(level_dof_indices_coarse_plain); size += MemoryConsumption::memory_consumption(level_dof_indices_fine); return size;