From: Peter Munch Date: Sun, 13 Sep 2020 09:43:45 +0000 (+0200) Subject: Fix categorize_by_boundary_ids for PBC X-Git-Tag: v9.3.0-rc1~1121^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=77dc76a5d6997783799468aa094e090131f802f3;p=dealii.git Fix categorize_by_boundary_ids for PBC --- diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h index ffcdd5b790..b898a5bf32 100644 --- a/include/deal.II/matrix_free/tools.h +++ b/include/deal.II/matrix_free/tools.h @@ -39,11 +39,24 @@ namespace MatrixFreeTools categorize_by_boundary_ids( const Triangulation &tria, AdditionalData & additional_data, - const unsigned int level = numbers::invalid_unsigned_int) + const unsigned int level = numbers::invalid_unsigned_int); + + + + // implementations + + + + template + void + categorize_by_boundary_ids(const Triangulation &tria, + AdditionalData & additional_data, + const unsigned int level) { - bool is_mg = (level != numbers::invalid_unsigned_int); + // ... determine if we are on an active or a multigrid level + const bool is_mg = (level != numbers::invalid_unsigned_int); - // ... create list for the category of each cell + // ... create empty list for the category of each cell if (is_mg) additional_data.cell_vectorization_category.resize( std::distance(tria.begin(level), tria.end(level))); @@ -51,28 +64,30 @@ namespace MatrixFreeTools additional_data.cell_vectorization_category.resize(tria.n_active_cells()); // ... set up scaling factor - std::vector factors(dim * 2); - - const auto bids = tria.get_boundary_ids(); + std::vector factors(GeometryInfo::faces_per_cell); - std::map bid_map; - for (const auto id : tria.get_boundary_ids()) - bid_map[id] = bid_map.size() + 1; + auto bids = tria.get_boundary_ids(); + std::sort(bids.begin(), bids.end()); { unsigned int n_bids = bids.size() + 1; int offset = 1; - for (unsigned int i = 0; i < dim * 2; i++, offset = offset * n_bids) + for (unsigned int i = 0; i < GeometryInfo::faces_per_cell; + i++, offset = offset * n_bids) factors[i] = offset; } - auto to_category = [&](auto &cell) { + const auto to_category = [&](const auto &cell) { unsigned int c_num = 0; - for (unsigned int i = 0; i < dim * 2; i++) + for (unsigned int i = 0; i < GeometryInfo::faces_per_cell; i++) { auto &face = *cell->face(i); - if (face.at_boundary()) - c_num += factors[i] * bid_map[face.boundary_id()]; + if (face.at_boundary() && !cell->has_periodic_neighbor(i)) + c_num += + factors[i] * (1 + std::distance(bids.begin(), + std::find(bids.begin(), + bids.end(), + face.boundary_id()))); } return c_num; };