From: Peter Munch Date: Tue, 24 Aug 2021 08:34:09 +0000 (+0200) Subject: Fix HangingNodes::setup_constraints() for multiple components X-Git-Tag: v9.4.0-rc1~1033^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F12702%2Fhead;p=dealii.git Fix HangingNodes::setup_constraints() for multiple components --- diff --git a/include/deal.II/matrix_free/hanging_nodes_internal.h b/include/deal.II/matrix_free/hanging_nodes_internal.h index 88a1c7c1e9..225b654d04 100644 --- a/include/deal.II/matrix_free/hanging_nodes_internal.h +++ b/include/deal.II/matrix_free/hanging_nodes_internal.h @@ -346,18 +346,22 @@ namespace internal auto &mask = masks[comp]; mask = ConstraintKinds::unconstrained; - const auto &fe = cell->get_fe().base_element(base_element_index); + const auto &fe_base = + cell->get_fe().base_element(base_element_index); - if (dim == 1 || dynamic_cast *>(&fe) == nullptr) + if (dim == 1 || + dynamic_cast *>(&fe_base) == nullptr) continue; - const unsigned int fe_degree = fe.tensor_degree(); + const unsigned int fe_degree = fe_base.tensor_degree(); const unsigned int n_dofs_1d = fe_degree + 1; const unsigned int dofs_per_face = Utilities::fixed_power(n_dofs_1d); std::vector neighbor_dofs_all( idx_offset.back()); + std::vector neighbor_dofs_all_temp( + idx_offset.back()); std::vector neighbor_dofs(dofs_per_face); @@ -659,20 +663,24 @@ namespace internal ExcInternalError(); // Copy the unconstrained values - neighbor_dofs.resize(n_dofs_1d * n_dofs_1d * - n_dofs_1d); DoFCellAccessor( &neighbor_cell->get_triangulation(), neighbor_cell->level(), neighbor_cell->index(), &cell->get_dof_handler()) - .get_dof_indices(neighbor_dofs); + .get_dof_indices(neighbor_dofs_all); // If the vector is distributed, we need to // transform the global indices to local ones. if (partitioner) - for (auto &index : neighbor_dofs) + for (auto &index : neighbor_dofs_all) index = partitioner->global_to_local(index); + for (unsigned int i = 0; + i < neighbor_dofs_all_temp.size(); + ++i) + neighbor_dofs_all_temp[i] = + neighbor_dofs_all[lexicographic_mapping[i]]; + for (unsigned int i = 0; i < n_dofs_1d; ++i) { // Get local dof index along line @@ -680,14 +688,12 @@ namespace internal line_dof_idx(local_line, i, n_dofs_1d); dof_indices[idx + idx_offset[comp]] = - neighbor_dofs - [lexicographic_mapping - [fe.component_to_system_index( - comp, - line_dof_idx( - local_line_neighbor, - flipped ? fe_degree - i : i, - n_dofs_1d))]]; + neighbor_dofs_all_temp + [line_dof_idx(local_line_neighbor, + flipped ? fe_degree - i : + i, + n_dofs_1d) + + idx_offset[comp]]; } // Stop looping over edge neighbors diff --git a/tests/matrix_free/matrix_vector_stokes_noflux2.cc b/tests/matrix_free/matrix_vector_stokes_noflux2.cc new file mode 100644 index 0000000000..c3d2ea5d69 --- /dev/null +++ b/tests/matrix_free/matrix_vector_stokes_noflux2.cc @@ -0,0 +1,340 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2013 - 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. +// +// --------------------------------------------------------------------- + + + +// this function tests the correctness of the implementation of matrix free +// matrix-vector products by comparing with the result of deal.II sparse +// matrix. The implementation is similar to matrix_vector_stokes_noflux.cc +// but tests edge hanging-node constraints. + +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include + +#include +#include + +#include "../tests.h" + + + +template +class MatrixFreeTest +{ +public: + using CellIterator = typename DoFHandler::active_cell_iterator; + using Number = double; + + MatrixFreeTest(const MatrixFree &data_in) + : data(data_in){}; + + void + local_apply(const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &cell_range) const + { + using vector_t = VectorizedArray; + FEEvaluation velocity(data, + 0); + FEEvaluation pressure(data, 1); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + velocity.reinit(cell); + velocity.read_dof_values(src.block(0)); + velocity.evaluate(EvaluationFlags::gradients); + pressure.reinit(cell); + pressure.read_dof_values(src.block(1)); + pressure.evaluate(EvaluationFlags::values); + + for (unsigned int q = 0; q < velocity.n_q_points; ++q) + { + SymmetricTensor<2, dim, vector_t> sym_grad_u = + velocity.get_symmetric_gradient(q); + vector_t pres = pressure.get_value(q); + vector_t div = -trace(sym_grad_u); + pressure.submit_value(div, q); + + // subtract p * I + for (unsigned int d = 0; d < dim; ++d) + sym_grad_u[d][d] -= pres; + + velocity.submit_symmetric_gradient(sym_grad_u, q); + } + + velocity.integrate(EvaluationFlags::gradients); + velocity.distribute_local_to_global(dst.block(0)); + pressure.integrate(EvaluationFlags::values); + pressure.distribute_local_to_global(dst.block(1)); + } + } + + + void + vmult(VectorType &dst, const VectorType &src) const + { + dst = 0; + data.cell_loop(&MatrixFreeTest::local_apply, + this, + dst, + src); + }; + +private: + const MatrixFree &data; +}; + + + +template +void +test() +{ + Triangulation triangulation; + GridGenerator::hyper_cube(triangulation); + triangulation.refine_global(); + for (const auto &cell : triangulation.active_cell_iterators()) + if (cell->index() > 0) + cell->set_refine_flag(); + triangulation.execute_coarsening_and_refinement(); + + MappingQ mapping(3); + FE_Q fe_u_scal(fe_degree + 1); + FESystem fe_u(fe_u_scal, dim); + FE_Q fe_p(fe_degree); + FESystem fe(fe_u_scal, dim, fe_p, 1); + DoFHandler dof_handler_u(triangulation); + DoFHandler dof_handler_p(triangulation); + DoFHandler dof_handler(triangulation); + + MatrixFree mf_data; + + AffineConstraints constraints, constraints_u, constraints_p; + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + + BlockVector solution; + BlockVector system_rhs; + BlockVector mf_solution; + + dof_handler.distribute_dofs(fe); + dof_handler_u.distribute_dofs(fe_u); + dof_handler_p.distribute_dofs(fe_p); + std::vector stokes_sub_blocks(dim + 1, 0); + stokes_sub_blocks[dim] = 1; + DoFRenumbering::component_wise(dof_handler, stokes_sub_blocks); + + std::set no_normal_flux_boundaries; + no_normal_flux_boundaries.insert(0); + no_normal_flux_boundaries.insert(1); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + VectorTools::compute_no_normal_flux_constraints( + dof_handler, 0, no_normal_flux_boundaries, constraints, mapping); + constraints.close(); + DoFTools::make_hanging_node_constraints(dof_handler_u, constraints_u); + VectorTools::compute_no_normal_flux_constraints( + dof_handler_u, 0, no_normal_flux_boundaries, constraints_u, mapping); + constraints_u.close(); + DoFTools::make_hanging_node_constraints(dof_handler_p, constraints_p); + constraints_p.close(); + + const std::vector dofs_per_block = + DoFTools::count_dofs_per_fe_block(dof_handler, stokes_sub_blocks); + + // std::cout << "Number of active cells: " + // << triangulation.n_active_cells() + // << std::endl + // << "Number of degrees of freedom: " + // << dof_handler.n_dofs() + // << " (" << n_u << '+' << n_p << ')' + // << std::endl; + + { + BlockDynamicSparsityPattern csp(2, 2); + + for (unsigned int d = 0; d < 2; ++d) + for (unsigned int e = 0; e < 2; ++e) + csp.block(d, e).reinit(dofs_per_block[d], dofs_per_block[e]); + + csp.collect_sizes(); + + DoFTools::make_sparsity_pattern(dof_handler, csp, constraints, false); + sparsity_pattern.copy_from(csp); + } + + system_matrix.reinit(sparsity_pattern); + + // this is from step-22 + { + QGauss quadrature_formula(fe_degree + 2); + + FEValues fe_values(mapping, + fe, + quadrature_formula, + update_values | update_JxW_values | + update_gradients); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix local_matrix(dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + const FEValuesExtractors::Vector velocities(0); + const FEValuesExtractors::Scalar pressure(dim); + + std::vector> phi_grads_u(dofs_per_cell); + std::vector div_phi_u(dofs_per_cell); + std::vector phi_p(dofs_per_cell); + + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell != endc; ++cell) + { + fe_values.reinit(cell); + local_matrix = 0; + + for (unsigned int q = 0; q < n_q_points; ++q) + { + for (unsigned int k = 0; k < dofs_per_cell; ++k) + { + phi_grads_u[k] = fe_values[velocities].symmetric_gradient(k, q); + div_phi_u[k] = fe_values[velocities].divergence(k, q); + phi_p[k] = fe_values[pressure].value(k, q); + } + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j <= i; ++j) + { + local_matrix(i, j) += + (phi_grads_u[i] * phi_grads_u[j] - + div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j]) * + fe_values.JxW(q); + } + } + } + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = i + 1; j < dofs_per_cell; ++j) + local_matrix(i, j) = local_matrix(j, i); + + cell->get_dof_indices(local_dof_indices); + constraints.distribute_local_to_global(local_matrix, + local_dof_indices, + system_matrix); + } + } + + + solution.reinit(2); + for (unsigned int d = 0; d < 2; ++d) + solution.block(d).reinit(dofs_per_block[d]); + solution.collect_sizes(); + + system_rhs.reinit(solution); + mf_solution.reinit(solution); + + // fill system_rhs with random numbers + for (unsigned int j = 0; j < system_rhs.block(0).size(); ++j) + if (constraints_u.is_constrained(j) == false) + { + const double val = -1 + 2. * random_value(); + system_rhs.block(0)(j) = val; + } + for (unsigned int j = 0; j < system_rhs.block(1).size(); ++j) + if (constraints_p.is_constrained(j) == false) + { + const double val = -1 + 2. * random_value(); + system_rhs.block(1)(j) = val; + } + + // setup matrix-free structure + { + std::vector *> dofs; + dofs.push_back(&dof_handler_u); + dofs.push_back(&dof_handler_p); + std::vector *> constraints; + constraints.push_back(&constraints_u); + constraints.push_back(&constraints_p); + QGauss<1> quad(fe_degree + 2); + // no parallelism + mf_data.reinit(mapping, + dofs, + constraints, + quad, + typename MatrixFree::AdditionalData( + MatrixFree::AdditionalData::none)); + } + + system_matrix.vmult(solution, system_rhs); + + MatrixFreeTest> mf(mf_data); + mf.vmult(mf_solution, system_rhs); + + // Verification + mf_solution -= solution; + const double error = mf_solution.linfty_norm(); + const double relative = solution.linfty_norm(); + deallog << "Verification fe degree " << fe_degree << ": " << error / relative + << std::endl + << std::endl; +} + + + +int +main() +{ + initlog(); + deallog << std::setprecision(3); + + { + deallog << std::endl << "Test with doubles" << std::endl << std::endl; + deallog.push("2d"); + test<2, 1>(); + test<2, 2>(); + test<2, 3>(); + deallog.pop(); + deallog.push("3d"); + test<3, 1>(); + deallog.pop(); + } +} diff --git a/tests/matrix_free/matrix_vector_stokes_noflux2.output b/tests/matrix_free/matrix_vector_stokes_noflux2.output new file mode 100644 index 0000000000..0434059632 --- /dev/null +++ b/tests/matrix_free/matrix_vector_stokes_noflux2.output @@ -0,0 +1,12 @@ + +DEAL:: +DEAL::Test with doubles +DEAL:: +DEAL:2d::Verification fe degree 1: 0 +DEAL:2d:: +DEAL:2d::Verification fe degree 2: 0 +DEAL:2d:: +DEAL:2d::Verification fe degree 3: 0 +DEAL:2d:: +DEAL:3d::Verification fe degree 1: 0 +DEAL:3d::