From: Martin Kronbichler Date: Mon, 16 Apr 2018 10:00:27 +0000 (+0200) Subject: Add internal function for identification of faces. X-Git-Tag: v9.0.0-rc1~148^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F6210%2Fhead;p=dealii.git Add internal function for identification of faces. --- diff --git a/include/deal.II/matrix_free/face_setup_internal.h b/include/deal.II/matrix_free/face_setup_internal.h new file mode 100644 index 0000000000..368364ca4a --- /dev/null +++ b/include/deal.II/matrix_free/face_setup_internal.h @@ -0,0 +1,1068 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +#ifndef dealii_face_setup_internal_h +#define dealii_face_setup_internal_h + +#include +#include +#include +#include +#include +#include + +#include + + +DEAL_II_NAMESPACE_OPEN + + +namespace internal +{ + namespace MatrixFreeFunctions + { + /** + * A struct that is used to represent a collection of faces of a process + * with one of its neighbor within the setup done in struct FaceInfo. + * + * @author Katharina Kormann, Martin Kronbichler, 2018 + */ + struct FaceIdentifier + { + FaceIdentifier() + : + n_hanging_faces_smaller_subdomain (0), + n_hanging_faces_larger_subdomain (0) + {} + + std::vector > shared_faces; + unsigned int n_hanging_faces_smaller_subdomain; + unsigned int n_hanging_faces_larger_subdomain; + }; + + + + /** + * A struct that extracts the faces relevant to a given set of cells, + * including the assignment of which of the two neighboring processors at + * a subdomain boundary with MPI should do the integration (from both + * sides). This data structure is used for the setup of the connectivity + * between faces and cells and for identification of the dof indices to be + * used for face integrals. + * + * @author Katharina Kormann, Martin Kronbichler, 2018 + */ + template + struct FaceSetup + { + FaceSetup(); + + /** + * Perform the initial detection of faces before reading the indices on + * the cells. This does not add the faces yet but only decides on + * whether some of the faces should be considered for processing + * locally. + */ + template + void initialize(const dealii::Triangulation &triangulation, + const MFAddData &additional_data, + std::vector > &cell_levels); + + /** + * Upon completion of the dof indices, this function extracts the + * information relevant for FaceToCellTopology and categorizes the faces + * into interior faces, boundary faces, and ghost faces (not processed + * locally but adjacent to some of the cells present locally). + */ + void + generate_faces(const dealii::Triangulation &triangulation, + const std::vector > &cell_levels, + TaskInfo &task_info); + + /** + * Fills the information about the cell, the face number, and numbers + * within the plain array representation in MatrixFree into + * FaceToCellTopology (without vectorization, which is something applied + * later). + */ + FaceToCellTopology<1> + create_face(const unsigned int face_no, + const typename dealii::Triangulation::cell_iterator &cell, + const unsigned int number_cell_interior, + const typename dealii::Triangulation::cell_iterator &neighbor, + const unsigned int number_cell_exterior); + + bool use_active_cells; + + /** + * A type that categorizes faces in the first initialize() function such + * that we can later get their correct value in generate_faces(). + */ + enum class FaceCategory : char + { + locally_active_at_boundary, + locally_active_done_here, + locally_active_done_elsewhere, + ghosted, + multigrid_refinement_edge + }; + + std::vector face_is_owned; + std::vector at_processor_boundary; + std::vector cells_close_to_boundary; + std::vector > inner_faces; + std::vector > boundary_faces; + std::vector > inner_ghost_faces; + std::vector > refinement_edge_faces; + }; + + + + /** + * Actually form the batches for vectorized execution of face integrals. + */ + template + void + collect_faces_vectorization + (const std::vector > &faces_in, + const std::vector &hard_vectorization_boundary, + std::vector &face_partition_data, + std::vector > &faces_out); + + + + /* -------------------------------------------------------------------- */ + +#ifndef DOXYGEN + + template + FaceSetup::FaceSetup() + : + use_active_cells (true) + {} + + + + template + template + void + FaceSetup::initialize(const dealii::Triangulation &triangulation, + const MFAddData &additional_data, + std::vector > &cell_levels) + { + use_active_cells = additional_data.level_mg_handler == numbers::invalid_unsigned_int; + +#ifdef DEBUG + // safety check + if (use_active_cells) + for (unsigned int i=0; i::cell_iterator + dcell(&triangulation, cell_levels[i].first, + cell_levels[i].second); + Assert(dcell->active(), ExcInternalError()); + } +#endif + + // step 1: add ghost cells for those cells that we identify as + // interesting + + at_processor_boundary.resize(cell_levels.size(), false); + cells_close_to_boundary.clear(); + face_is_owned.resize(dim > 1 ? triangulation.n_raw_faces() : + triangulation.n_vertices(), FaceCategory::locally_active_done_elsewhere); + + // go through the mesh and divide the faces on the processor + // boundaries as evenly as possible between the processors + std::map inner_faces_at_proc_boundary; + if (dynamic_cast*>(&triangulation)) + { + const types::subdomain_id my_domain = triangulation.locally_owned_subdomain(); + for (unsigned int i=0; i0 && cell_levels[i] == cell_levels[i-1]) + continue; + typename dealii::Triangulation::cell_iterator + dcell(&triangulation, cell_levels[i].first, + cell_levels[i].second); + for (unsigned int f=0; f::faces_per_cell; ++f) + { + if (dcell->at_boundary(f) && !dcell->has_periodic_neighbor(f)) + continue; + typename dealii::Triangulation::cell_iterator + neighbor = dcell->neighbor_or_periodic_neighbor(f); + + // faces at hanging nodes are always treated by the processor + // who owns the element on the fine side. but we need to count + // the number of inner faces in order to balance the remaining + // faces properly + const CellId id_mine = dcell->id(); + if (use_active_cells && neighbor->has_children()) + for (unsigned int c=0; cface(f)->n_children(); ++c) + { + typename dealii::Triangulation::cell_iterator + neighbor_c = dcell->at_boundary(f) ? + dcell->periodic_neighbor_child_on_subface(f, c) : + dcell->neighbor_child_on_subface(f, c); + const types::subdomain_id neigh_domain = neighbor_c->subdomain_id(); + if (my_domain < neigh_domain) + inner_faces_at_proc_boundary[neigh_domain].n_hanging_faces_larger_subdomain++; + else if (my_domain > neigh_domain) + inner_faces_at_proc_boundary[neigh_domain].n_hanging_faces_smaller_subdomain++; + } + else + { + const types::subdomain_id neigh_domain = use_active_cells ? + neighbor->subdomain_id() : + neighbor->level_subdomain_id(); + if (neighbor->level() < dcell->level() && + use_active_cells) + { + if (my_domain < neigh_domain) + inner_faces_at_proc_boundary[neigh_domain].n_hanging_faces_smaller_subdomain++; + else if (my_domain > neigh_domain) + inner_faces_at_proc_boundary[neigh_domain].n_hanging_faces_larger_subdomain++; + } + else if (neighbor->level() == dcell->level() && + my_domain != neigh_domain) + { + // always list the cell whose owner has the lower + // subdomain id first. this applies to both processors + // involved, so both processors will generate the same + // list that we will later order + const CellId id_neigh = neighbor->id(); + inner_faces_at_proc_boundary[neigh_domain].shared_faces. + push_back(my_domain < neigh_domain + ? + std::make_pair(id_mine, id_neigh) + : + std::make_pair(id_neigh, id_mine)); + } + } + } + } + + // sort the cell ids related to each neighboring processor. This + // algorithm is symmetric so every processor combination should + // arrive here and no deadlock should be possible + for (std::map::iterator it + = inner_faces_at_proc_boundary.begin(); + it != inner_faces_at_proc_boundary.end(); ++it) + { + Assert(it->first != my_domain, + ExcInternalError("Should not send info to myself")); + std::sort(it->second.shared_faces.begin(), it->second.shared_faces.end()); + it->second.shared_faces.erase(std::unique(it->second.shared_faces.begin(), + it->second.shared_faces.end()), + it->second.shared_faces.end()); + + // safety check: both involved processors should see the same list + // because the pattern of ghosting is symmetric. We test this by + // looking at the length of the lists of faces +#if defined(DEAL_II_WITH_MPI) && defined(DEBUG) + MPI_Comm comm = MPI_COMM_SELF; + if (const parallel::Triangulation *ptria = + dynamic_cast *>(&triangulation)) + comm = ptria->get_communicator(); + + MPI_Status status; + unsigned int mysize = it->second.shared_faces.size(); + unsigned int othersize = numbers::invalid_unsigned_int; + MPI_Sendrecv(&mysize, 1, MPI_UNSIGNED, it->first, 600+my_domain, + &othersize, 1, MPI_UNSIGNED, it->first, 600+it->first, + comm, &status); + AssertDimension(mysize, othersize); + mysize = it->second.n_hanging_faces_smaller_subdomain; + MPI_Sendrecv(&mysize, 1, MPI_UNSIGNED, it->first, 700+my_domain, + &othersize, 1, MPI_UNSIGNED, it->first, 700+it->first, + comm, &status); + AssertDimension(mysize, othersize); + mysize = it->second.n_hanging_faces_larger_subdomain; + MPI_Sendrecv(&mysize, 1, MPI_UNSIGNED, it->first, 800+my_domain, + &othersize, 1, MPI_UNSIGNED, it->first, 800+it->first, + comm, &status); + AssertDimension(mysize, othersize); +#endif + + // Arrange the face "ownership" such that cells that are access + // by more than one face (think of a cell in a corner) get + // ghosted. This arrangement has the advantage that we need to + // send less data because the same data is used twice. The + // strategy applied here is to ensure the same order of face + // pairs on both processors that share some faces, and make the + // same decision on both sides. + + // Create a vector with cell ids sorted over the processor with + // the larger rank. In the code below we need to be able to + // identify the same cell once for the processor with higher + // rank and once for the processor with the lower rank. The + // format for the processor with the higher rank is already + // contained in `shared_faces`, whereas we need a copy that we + // sort differently for the other way around. + std::vector > other_range(it->second.shared_faces.size()); + for (unsigned int i=0; isecond.shared_faces[i].second, + it->second.shared_faces[i].first, + i); + std::sort(other_range.begin(), other_range.end()); + + // the vector 'assignment' sets whether a particular cell + // appears more often and acts as a pre-selection of the rank. A + // value of 1 means that the process with the higher rank gets + // those faces, a value -1 means that the process with the lower + // rank gets it, whereas a value 0 means that the decision can + // be made in an arbitrary way. + unsigned int n_faces_lower_proc = 0, n_faces_higher_proc = 0; + std::vector assignment(other_range.size(), 0); + if (it->second.shared_faces.size() > 0) + { + // identify faces that go to the processor with the higher + // rank + unsigned int count = 0; + for (unsigned int i=1; isecond.shared_faces.size(); ++i) + if (it->second.shared_faces[i].first == + it->second.shared_faces[i-1-count].first) + ++count; + else + { + AssertThrow(count < 2*dim, ExcInternalError()); + if (count > 0) + { + for (unsigned int k=0; k<=count; ++k) + assignment[i-1-k] = 1; + n_faces_higher_proc += count+1; + } + count = 0; + } + + // identify faces that definitely go to the processor with + // the lower rank - this must use the sorting of CellId + // variables from the processor with the higher rank, i.e., + // other_range rather than `shared_faces`. + count = 0; + for (unsigned int i=1; i(other_range[i]) == + std::get<0>(other_range[i-1-count])) + ++count; + else + { + AssertThrow(count < 2*dim, ExcInternalError()); + if (count > 0) + { + for (unsigned int k=0; k<=count; ++k) + { + Assert(it->second.shared_faces[std::get<2>(other_range[i-1])].second == + it->second.shared_faces[std::get<2>(other_range[i-1-k])].second, + ExcInternalError()); + // only assign to -1 if higher rank was not + // yet set + if (assignment[std::get<2>(other_range[i-1-k])] == 0) + { + assignment[std::get<2>(other_range[i-1-k])] = -1; + ++n_faces_lower_proc; + } + } + } + count = 0; + } + } + + + // divide the faces evenly between the two processors. the + // processor with small rank takes the first half, the processor + // with larger rank the second half. Adjust for the hanging + // faces that always get assigned to one side, and the faces we + // have already assigned due to the criterion above + n_faces_lower_proc += it->second.n_hanging_faces_smaller_subdomain; + n_faces_higher_proc += it->second.n_hanging_faces_larger_subdomain; + const unsigned int n_total_faces_at_proc_boundary = + (it->second.shared_faces.size() + + it->second.n_hanging_faces_smaller_subdomain + + it->second.n_hanging_faces_larger_subdomain); + unsigned int split_index = n_total_faces_at_proc_boundary/2; + if (split_index < n_faces_lower_proc) + split_index = 0; + else if (split_index < n_total_faces_at_proc_boundary - + n_faces_higher_proc) + split_index -= n_faces_lower_proc; + else + split_index = n_total_faces_at_proc_boundary - n_faces_higher_proc - n_faces_lower_proc; + + // make sure the splitting is consistent between both sides +#if defined(DEAL_II_WITH_MPI) && defined(DEBUG) + MPI_Sendrecv(&split_index, 1, MPI_UNSIGNED, it->first, 900+my_domain, + &othersize, 1, MPI_UNSIGNED, it->first, 900+it->first, + comm, &status); + AssertDimension(split_index, othersize); + MPI_Sendrecv(&n_faces_lower_proc, 1, MPI_UNSIGNED, it->first, 1000+my_domain, + &othersize, 1, MPI_UNSIGNED, it->first, 1000+it->first, + comm, &status); + AssertDimension(n_faces_lower_proc, othersize); + MPI_Sendrecv(&n_faces_higher_proc, 1, MPI_UNSIGNED, it->first, 1100+my_domain, + &othersize, 1, MPI_UNSIGNED, it->first, 1100+it->first, + comm, &status); + AssertDimension(n_faces_higher_proc, othersize); +#endif + + // collect the faces on both sides + std::vector > owned_faces_lower, owned_faces_higher; + for (unsigned int i=0; isecond.shared_faces[i]); + else if (assignment[i] > 0) + owned_faces_higher.push_back(it->second.shared_faces[i]); + AssertIndexRange(split_index, it->second.shared_faces.size() + 1 - + owned_faces_lower.size() - owned_faces_higher.size()); + + unsigned int i=0, c=0; + for (; isecond.shared_faces[i]); + ++c; + } + for ( ; isecond.shared_faces[i]); + } + +#ifdef DEBUG + // check consistency of faces on both sides + std::vector > check_faces; + check_faces.insert(check_faces.end(), owned_faces_lower.begin(), + owned_faces_lower.end()); + check_faces.insert(check_faces.end(), owned_faces_higher.begin(), + owned_faces_higher.end()); + std::sort(check_faces.begin(), check_faces.end()); + AssertDimension(check_faces.size(), it->second.shared_faces.size()); + for (unsigned int i=0; isecond.shared_faces[i], ExcInternalError()); +#endif + + // now only set half of the faces as the ones to keep + if (my_domain < it->first) + it->second.shared_faces.swap(owned_faces_lower); + else + it->second.shared_faces.swap(owned_faces_higher); + + std::sort(it->second.shared_faces.begin(), it->second.shared_faces.end()); + } + } + + // fill in the additional cells that we need access to via ghosting to + // cell_levels + std::set > ghost_cells; + for (unsigned int i=0; i::cell_iterator + dcell(&triangulation, cell_levels[i].first, + cell_levels[i].second); + if (use_active_cells) + Assert(dcell->active(), ExcNotImplemented()); + for (unsigned int f=0; f::faces_per_cell; ++f) + { + if (dcell->at_boundary(f) && !dcell->has_periodic_neighbor(f)) + face_is_owned[dcell->face(f)->index()] = FaceCategory::locally_active_at_boundary; + + // treat boundaries of cells of different refinement level + // inside the domain in case of multigrid separately + else if ((dcell->at_boundary(f) == false || + dcell->has_periodic_neighbor(f)) && + additional_data.level_mg_handler != numbers::invalid_unsigned_int + && + dcell->neighbor_or_periodic_neighbor(f)->level() < dcell->level()) + { + face_is_owned[dcell->face(f)->index()] = FaceCategory::multigrid_refinement_edge; + } + else + { + typename dealii::Triangulation::cell_iterator neighbor = + dcell->neighbor_or_periodic_neighbor(f); + + // neighbor is refined -> face will be treated by neighbor + if (use_active_cells && neighbor->has_children() && + additional_data.hold_all_faces_to_owned_cells == false) + continue; + + bool add_to_ghost = false; + const types::subdomain_id + id1 = use_active_cells ? dcell->subdomain_id() : dcell->level_subdomain_id(), + id2 = use_active_cells ? (neighbor->has_children() ? + dcell->neighbor_child_on_subface(f,0)->subdomain_id() : + neighbor->subdomain_id()) : + neighbor->level_subdomain_id(); + + // Check whether the current face should be processed + // locally (instead of being processed from the other + // side). We process a face locally when we are more refined + // (in the active cell case) or when the face is listed in + // the `shared_faces` data structure that we built above. + if ((id1 == id2 && (use_active_cells == false || + neighbor->has_children() == false)) || + dcell->level() > neighbor->level() || + std::binary_search(inner_faces_at_proc_boundary[id2].shared_faces.begin(), + inner_faces_at_proc_boundary[id2].shared_faces.end(), + std::make_pair(id1id() : neighbor->id(), + id1id() : dcell->id()))) + { + face_is_owned[dcell->face(f)->index()] = FaceCategory::locally_active_done_here; + if (dcell->level() == neighbor->level()) + face_is_owned[neighbor->face(dcell->has_periodic_neighbor(f) ? + dcell->periodic_neighbor_face_no(f) : + dcell->neighbor_face_no(f))->index()] + = FaceCategory::locally_active_done_here; + + // If neighbor is a ghost element (i.e. dcell->subdomain_id + // ! dcell->neighbor(f)->subdomain_id()), we need to add its + // index into cell level list. + if (use_active_cells) + add_to_ghost = (dcell->subdomain_id() != + neighbor->subdomain_id()); + else + add_to_ghost = (dcell->level_subdomain_id() != + neighbor->level_subdomain_id()); + } + else if (additional_data.hold_all_faces_to_owned_cells == false) + { + // mark the cell to be close to the boundary + cells_close_to_boundary.emplace_back(i); + } + else + { + // add all cells to ghost layer... + face_is_owned[dcell->face(f)->index()] = FaceCategory::ghosted; + if (use_active_cells) + { + if (neighbor->has_children()) + for (unsigned int s=0; sface(f)->n_children(); ++s) + if (dcell->at_boundary(f)) + { + if (dcell->periodic_neighbor_child_on_subface(f,s)->subdomain_id() != + dcell->subdomain_id()) + add_to_ghost = true; + } + else + { + if (dcell->neighbor_child_on_subface(f,s)->subdomain_id() != + dcell->subdomain_id()) + add_to_ghost = true; + } + else + add_to_ghost = (dcell->subdomain_id() != + neighbor->subdomain_id()); + } + else + add_to_ghost = (dcell->level_subdomain_id() != + neighbor->level_subdomain_id()); + } + + if (add_to_ghost) + { + ghost_cells.insert + (std::pair + (neighbor->level(), neighbor->index())); + at_processor_boundary[i] = true; + } + } + } + } + + // step 2: append the ghost cells at the end of the locally owned + // cells + for (std::set >::iterator + it=ghost_cells.begin(); it!=ghost_cells.end(); ++it) + cell_levels.push_back(*it); + + // step 3: clean up the cells close to the boundary + std::sort(cells_close_to_boundary.begin(), cells_close_to_boundary.end()); + cells_close_to_boundary.erase(std::unique(cells_close_to_boundary.begin(), + cells_close_to_boundary.end()), + cells_close_to_boundary.end()); + std::vector final_cells; + final_cells.reserve(cells_close_to_boundary.size()); + for (unsigned int i=0; i + void + FaceSetup + ::generate_faces(const dealii::Triangulation &triangulation, + const std::vector > &cell_levels, + TaskInfo &task_info) + { + // step 1: create the inverse map between cell iterators and the + // cell_level_index field + std::map, unsigned int> + map_to_vectorized; + for (unsigned int cell=0; cell::cell_iterator dcell + (&triangulation, cell_levels[cell].first, cell_levels[cell].second); + std::pair level_index(dcell->level(), + dcell->index()); + map_to_vectorized[level_index] = cell; + } + + // step 2: fill the information about inner faces and boundary faces + const unsigned int vectorization_length = task_info.vectorization_length; + task_info.face_partition_data.resize(task_info.cell_partition_data.size()-1, 0); + task_info.boundary_partition_data.resize(task_info.cell_partition_data.size()-1, 0); + std::vector face_visited(face_is_owned.size(), 0); + for (unsigned int partition=0; partition::cell_iterator dcell + (&triangulation, cell_levels[cell].first, cell_levels[cell].second); + for (unsigned int f=0; f::faces_per_cell; ++f) + { + // boundary face + if (face_is_owned[dcell->face(f)->index()] == FaceCategory::locally_active_at_boundary) + { + Assert(dcell->at_boundary(f), ExcInternalError()); + ++boundary_counter; + FaceToCellTopology<1> info; + info.cells_interior[0] = cell; + info.cells_exterior[0] = numbers::invalid_unsigned_int; + info.interior_face_no = f; + info.exterior_face_no = dcell->face(f)->boundary_id(); + info.subface_index = GeometryInfo::max_children_per_cell; + info.face_orientation = 0; + boundary_faces.push_back(info); + + face_visited[dcell->face(f)->index()]++; + } + // interior face, including faces over periodic boundaries + else + { + typename dealii::Triangulation::cell_iterator neighbor + = dcell->neighbor_or_periodic_neighbor(f); + if (use_active_cells && neighbor->has_children()) + { + for (unsigned int c=0; cface(f)->n_children(); ++c) + { + typename dealii::Triangulation::cell_iterator + neighbor_c = dcell->at_boundary(f) ? + dcell->periodic_neighbor_child_on_subface(f, c) : + dcell->neighbor_child_on_subface(f, c); + const types::subdomain_id neigh_domain = neighbor_c->subdomain_id(); + const unsigned int neighbor_face_no = + dcell->has_periodic_neighbor(f) ? + dcell->periodic_neighbor_face_no(f) : + dcell->neighbor_face_no(f); + if (neigh_domain != dcell->subdomain_id() || + face_visited[dcell->face(f)->child(c)->index()] == 1) + { + std::pair + level_index(neighbor_c->level(), neighbor_c->index()); + if (face_is_owned[dcell->face(f)->child(c)->index()]== + FaceCategory::locally_active_done_here) + { + ++inner_counter; + inner_faces. + push_back(create_face(neighbor_face_no, neighbor_c, + map_to_vectorized[level_index], + dcell, cell)); + + } + else if (face_is_owned[dcell->face(f)->child(c)->index()] == + FaceCategory::ghosted) + { + inner_ghost_faces. + push_back(create_face(neighbor_face_no, neighbor_c, + map_to_vectorized[level_index], + dcell, cell)); + } + else + Assert(face_is_owned[dcell->face(f)->index()] == + FaceCategory::locally_active_done_elsewhere, + ExcInternalError()); + } + else + { + face_visited[dcell->face(f)->child(c)->index()] = 1; + } + } + } + else + { + const types::subdomain_id my_domain = use_active_cells ? + dcell->subdomain_id() : dcell->level_subdomain_id(); + const types::subdomain_id neigh_domain = use_active_cells ? + neighbor->subdomain_id() : neighbor->level_subdomain_id(); + if (neigh_domain != my_domain || + face_visited[dcell->face(f)->index()] == 1) + { + std::pair + level_index(neighbor->level(), neighbor->index()); + if (face_is_owned[dcell->face(f)->index()]== + FaceCategory::locally_active_done_here) + { + ++inner_counter; + inner_faces. + push_back(create_face(f, dcell, cell, neighbor, + map_to_vectorized[level_index])); + } + else if (face_is_owned[dcell->face(f)->index()] == + FaceCategory::ghosted) + { + inner_ghost_faces. + push_back(create_face(f, dcell, cell, neighbor, + map_to_vectorized[level_index])); + } + } + else + { + face_visited[dcell->face(f)->index()] = 1; + if (dcell->has_periodic_neighbor(f)) + face_visited[neighbor->face(dcell->periodic_neighbor_face_no(f))->index()] = 1; + } + if (face_is_owned[dcell->face(f)->index()] == + FaceCategory::multigrid_refinement_edge) + { + refinement_edge_faces. + push_back(create_face(f, dcell, cell, neighbor, + refinement_edge_faces.size())); + } + } + } + } + } + task_info.face_partition_data[partition+1]= + task_info.face_partition_data[partition] + inner_counter; + task_info.boundary_partition_data[partition+1] = + task_info.boundary_partition_data[partition] + boundary_counter; + } + task_info.ghost_face_partition_data.resize(2); + task_info.ghost_face_partition_data[0] = 0; + task_info.ghost_face_partition_data[1] = inner_ghost_faces.size(); + task_info.refinement_edge_face_partition_data.resize(2); + task_info.refinement_edge_face_partition_data[0] = 0; + task_info.refinement_edge_face_partition_data[1] = refinement_edge_faces.size(); + } + + + + template + FaceToCellTopology<1> + FaceSetup + ::create_face(const unsigned int face_no, + const typename dealii::Triangulation::cell_iterator &cell, + const unsigned int number_cell_interior, + const typename dealii::Triangulation::cell_iterator &neighbor, + const unsigned int number_cell_exterior) + { + FaceToCellTopology<1> info; + info.cells_interior[0] = number_cell_interior; + info.cells_exterior[0] = number_cell_exterior; + info.interior_face_no = face_no; + if (cell->has_periodic_neighbor(face_no)) + info.exterior_face_no = cell->periodic_neighbor_face_no(face_no); + else + info.exterior_face_no = cell->neighbor_face_no(face_no); + + info.subface_index = GeometryInfo::max_children_per_cell; + Assert(neighbor->level() <= cell->level(), + ExcInternalError()); + if (cell->level() > neighbor->level()) + { + if (cell->has_periodic_neighbor(face_no)) + info.subface_index = + cell->periodic_neighbor_of_coarser_periodic_neighbor(face_no).second; + else + info.subface_index = + cell->neighbor_of_coarser_neighbor(face_no).second; + } + + info.face_orientation = 0; + unsigned int left_face_orientation = + !cell->face_orientation(face_no) + 2 * cell->face_flip(face_no) + + 4 * cell->face_rotation(face_no); + unsigned int right_face_orientation = + !neighbor->face_orientation(info.exterior_face_no) + + 2 * neighbor->face_flip(info.exterior_face_no) + + 4 * neighbor->face_rotation(info.exterior_face_no); + if (left_face_orientation != 0) + { + info.face_orientation = 8 + left_face_orientation; + Assert(right_face_orientation == 0, + ExcMessage("Face seems to be wrongly oriented from both sides")); + } + else + info.face_orientation = right_face_orientation; + return info; + } + + + + /** + * This simple comparison for collect_faces_vectorization() identifies + * faces of the same type, i.e., where all of the interior and exterior + * face number, subface index and orientation are the same. This is used + * to batch similar faces together for vectorization. + */ + bool compare_faces_for_vectorization + (const FaceToCellTopology<1> &face1, + const FaceToCellTopology<1> &face2) + { + if (face1.interior_face_no != face2.interior_face_no) + return false; + if (face1.exterior_face_no != face2.exterior_face_no) + return false; + if (face1.subface_index != face2.subface_index) + return false; + if (face1.face_orientation != face2.face_orientation) + return false; + return true; + } + + + + /** + * This comparator is used within collect_faces_vectorization() to create + * a sorting of FaceToCellTopology objects based on their + * identifiers. This is used to obtain a good data locality when + * processing the face integrals. + */ + template + struct FaceComparator + { + bool operator() (const FaceToCellTopology &face1, + const FaceToCellTopology &face2) + { + for (unsigned int i=0; i face2.cells_interior[i]) + return false; + for (unsigned int i=0; i face2.cells_exterior[i]) + return false; + if (face1.interior_face_no < face2.interior_face_no) + return true; + else if (face1.interior_face_no > face2.interior_face_no) + return false; + if (face1.exterior_face_no < face2.exterior_face_no) + return true; + else if (face1.exterior_face_no > face2.exterior_face_no) + return false; + + // we do not need to check for subface_index and orientation because + // those cannot be different if when all the other values are the + // same. + AssertDimension(face1.subface_index, face2.subface_index); + AssertDimension(face1.face_orientation, face2.face_orientation); + + return false; + } + }; + + + + template + void + collect_faces_vectorization + (const std::vector > &faces_in, + const std::vector &hard_vectorization_boundary, + std::vector &face_partition_data, + std::vector > &faces_out) + { + FaceToCellTopology macro_face; + std::vector > faces_type; + + unsigned int face_start = face_partition_data[0], + face_end = face_partition_data[0]; + + face_partition_data[0] = faces_out.size(); + for ( unsigned int partition=0; partition > new_faces_type; + + // start with the end point for the last partition + face_start = face_end; + face_end = face_partition_data[partition+1]; + + // set the partitioner to the new vectorized lengths + face_partition_data[partition+1] = face_partition_data[partition]; + + // loop over the faces in the current partition and reorder according to + // the face type + for (unsigned int face=face_start; face(1, face)); +face_found : + {} + } + + // insert new faces in sorted list to get good data locality + std::set, + FaceComparator > new_faces; + for (unsigned int type=0; type touched(no_faces, 0); + + // do two passes through the data. The first is to identify similar + // faces within the same index range as the cells which will allow + // for vectorized read operations, the second picks up all the rest + unsigned int n_vectorized = 0; + for (unsigned int f=0; f no_faces) + is_contiguous = false; + else + for (unsigned int v=1; v untouched; + untouched.reserve(no_faces - n_vectorized); + for (unsigned int f=0; f 0 && v < vectorization_width) + { + // must add non-filled face + if (hard_vectorization_boundary[partition+1] || + partition == face_partition_data.size()-2) + { + for ( ; v < vectorization_width; ++v) + { + // Dummy cell, not used + macro_face.cells_interior[v] = numbers::invalid_unsigned_int; + macro_face.cells_exterior[v] = numbers::invalid_unsigned_int; + } + new_faces.insert(macro_face); + } + else + { + // postpone to the next partition + std::vector untreated(v); + for (unsigned int f=0; f> in_faces, out_faces; + for (unsigned int i=0; i