From: David Wells Date: Sun, 28 Jan 2018 23:49:55 +0000 (-0500) Subject: Split the grid tools instantiations. X-Git-Tag: v9.0.0-rc1~464^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e28beee3c9ecf26e97399b2d64a0a5402ec5a4dc;p=dealii.git Split the grid tools instantiations. This commit splits the grid tools instantiation files in two: one file deals with instantiations for DoFHandlers and Triangulations and the other deals with instantiations of just Triangulations. Here is the timing information: before this patch: debug: grid tools: 74s, 3.6 GB release: grid tools: 97s 3.9 GB after this patch: debug: grid tools: 31s, 1.8 GB grid tools dof handlers: 48s, 2.4 GB release: grid tools: 46s, 1.9 GB grid tools dof handlers: 65s, 2.4 GB This lowers the high water mark for memory usage with only a small impact on the total compilation time. --- diff --git a/source/grid/CMakeLists.txt b/source/grid/CMakeLists.txt index a5888e4730..a1e4c1e165 100644 --- a/source/grid/CMakeLists.txt +++ b/source/grid/CMakeLists.txt @@ -37,6 +37,7 @@ SET(_unity_include_src SET(_separate_src grid_reordering.cc grid_tools.cc + grid_tools_dof_handlers.cc tria.cc ) @@ -56,6 +57,7 @@ SET(_inst grid_out.inst.in grid_refinement.inst.in grid_tools.inst.in + grid_tools_dof_handlers.inst.in grid_tools_cache.inst.in intergrid_map.inst.in manifold.inst.in diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 43e7f53cbf..3dd3083319 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -44,7 +44,6 @@ #include #include #include -#include #include #include @@ -1068,477 +1067,6 @@ namespace GridTools - template class MeshType, int spacedim> - unsigned int - find_closest_vertex (const MeshType &mesh, - const Point &p, - const std::vector &marked_vertices) - { - // first get the underlying - // triangulation from the - // mesh and determine vertices - // and used vertices - const Triangulation &tria = mesh.get_triangulation(); - - const std::vector< Point > &vertices = tria.get_vertices(); - - Assert ( tria.get_vertices().size() == marked_vertices.size() || marked_vertices.size() ==0, - ExcDimensionMismatch(tria.get_vertices().size(), marked_vertices.size())); - - // If p is an element of marked_vertices, - // and q is that of used_Vertices, - // the vector marked_vertices does NOT - // contain unused vertices if p implies q. - // I.e., if p is true q must be true - // (if p is false, q could be false or true). - // p implies q logic is encapsulated in ~p|q. - Assert( marked_vertices.size()==0 - || - std::equal( marked_vertices.begin(), - marked_vertices.end(), - tria.get_used_vertices().begin(), - [](bool p, bool q) - { - return !p || q; - }), - ExcMessage("marked_vertices should be a subset of used vertices in the triangulation " - "but marked_vertices contains one or more vertices that are not used vertices!") ); - - // In addition, if a vector bools - // is specified (marked_vertices) - // marking all the vertices which - // could be the potentially closest - // vertex to the point, use it instead - // of used vertices - const std::vector &used = - (marked_vertices.size()==0) ? tria.get_used_vertices() : marked_vertices; - - // At the beginning, the first - // used vertex is the closest one - std::vector::const_iterator first = - std::find(used.begin(), used.end(), true); - - // Assert that at least one vertex - // is actually used - Assert(first != used.end(), ExcInternalError()); - - unsigned int best_vertex = std::distance(used.begin(), first); - double best_dist = (p - vertices[best_vertex]).norm_square(); - - // For all remaining vertices, test - // whether they are any closer - for (unsigned int j = best_vertex+1; j < vertices.size(); j++) - if (used[j]) - { - double dist = (p - vertices[j]).norm_square(); - if (dist < best_dist) - { - best_vertex = j; - best_dist = dist; - } - } - - return best_vertex; - } - - - - template class MeshType, int spacedim> - unsigned int - find_closest_vertex (const Mapping &mapping, - const MeshType &mesh, - const Point &p, - const std::vector &marked_vertices) - { - // Take a shortcut in the simple case. - if (mapping.preserves_vertex_locations() == true) - return find_closest_vertex(mesh, p, marked_vertices); - - // first get the underlying - // triangulation from the - // mesh and determine vertices - // and used vertices - const Triangulation &tria = mesh.get_triangulation(); - - auto vertices = extract_used_vertices(tria, mapping); - - Assert ( tria.get_vertices().size() == marked_vertices.size() || marked_vertices.size() ==0, - ExcDimensionMismatch(tria.get_vertices().size(), marked_vertices.size())); - - // If p is an element of marked_vertices, - // and q is that of used_Vertices, - // the vector marked_vertices does NOT - // contain unused vertices if p implies q. - // I.e., if p is true q must be true - // (if p is false, q could be false or true). - // p implies q logic is encapsulated in ~p|q. - Assert( marked_vertices.size()==0 - || - std::equal( marked_vertices.begin(), - marked_vertices.end(), - tria.get_used_vertices().begin(), - [](bool p, bool q) - { - return !p || q; - }), - ExcMessage("marked_vertices should be a subset of used vertices in the triangulation " - "but marked_vertices contains one or more vertices that are not used vertices!") ); - - // Remove from the map unwanted elements. - if (marked_vertices.size()) - for (auto it = vertices.begin(); it != vertices.end(); ) - { - if (marked_vertices[it->first] == false) - { - vertices.erase(it++); - } - else - { - ++it; - } - } - - return find_closest_vertex(vertices, p); - } - - - - template class MeshType, int spacedim> -#ifndef _MSC_VER - std::vector::active_cell_iterator> -#else - std::vector >::type> -#endif - find_cells_adjacent_to_vertex(const MeshType &mesh, - const unsigned int vertex) - { - // make sure that the given vertex is - // an active vertex of the underlying - // triangulation - Assert(vertex < mesh.get_triangulation().n_vertices(), - ExcIndexRange(0,mesh.get_triangulation().n_vertices(),vertex)); - Assert(mesh.get_triangulation().get_used_vertices()[vertex], - ExcVertexNotUsed(vertex)); - - // use a set instead of a vector - // to ensure that cells are inserted only - // once - std::set >::type> adjacent_cells; - - typename dealii::internal::ActiveCellIterator >::type - cell = mesh.begin_active(), - endc = mesh.end(); - - // go through all active cells and look if the vertex is part of that cell - // - // in 1d, this is all we need to care about. in 2d/3d we also need to worry - // that the vertex might be a hanging node on a face or edge of a cell; in - // this case, we would want to add those cells as well on whose faces the - // vertex is located but for which it is not a vertex itself. - // - // getting this right is a lot simpler in 2d than in 3d. in 2d, a hanging - // node can only be in the middle of a face and we can query the neighboring - // cell from the current cell. on the other hand, in 3d a hanging node - // vertex can also be on an edge but there can be many other cells on - // this edge and we can not access them from the cell we are currently - // on. - // - // so, in the 3d case, if we run the algorithm as in 2d, we catch all - // those cells for which the vertex we seek is on a *subface*, but we - // miss the case of cells for which the vertex we seek is on a - // sub-edge for which there is no corresponding sub-face (because the - // immediate neighbor behind this face is not refined), see for example - // the bits/find_cells_adjacent_to_vertex_6 testcase. thus, if we - // haven't yet found the vertex for the current cell we also need to - // look at the mid-points of edges - // - // as a final note, deciding whether a neighbor is actually coarser is - // simple in the case of isotropic refinement (we just need to look at - // the level of the current and the neighboring cell). however, this - // isn't so simple if we have used anisotropic refinement since then - // the level of a cell is not indicative of whether it is coarser or - // not than the current cell. ultimately, we want to add all cells on - // which the vertex is, independent of whether they are coarser or - // finer and so in the 2d case below we simply add *any* *active* neighbor. - // in the worst case, we add cells multiple times to the adjacent_cells - // list, but std::set throws out those cells already entered - for (; cell != endc; ++cell) - { - for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; v++) - if (cell->vertex_index(v) == vertex) - { - // OK, we found a cell that contains - // the given vertex. We add it - // to the list. - adjacent_cells.insert(cell); - - // as explained above, in 2+d we need to check whether - // this vertex is on a face behind which there is a - // (possibly) coarser neighbor. if this is the case, - // then we need to also add this neighbor - if (dim >= 2) - for (unsigned int vface = 0; vface < dim; vface++) - { - const unsigned int face = - GeometryInfo::vertex_to_face[v][vface]; - - if (!cell->at_boundary(face) - && - cell->neighbor(face)->active()) - { - // there is a (possibly) coarser cell behind a - // face to which the vertex belongs. the - // vertex we are looking at is then either a - // vertex of that coarser neighbor, or it is a - // hanging node on one of the faces of that - // cell. in either case, it is adjacent to the - // vertex, so add it to the list as well (if - // the cell was already in the list then the - // std::set makes sure that we get it only - // once) - adjacent_cells.insert (cell->neighbor(face)); - } - } - - // in any case, we have found a cell, so go to the next cell - goto next_cell; - } - - // in 3d also loop over the edges - if (dim >= 3) - { - for (unsigned int e=0; e::lines_per_cell; ++e) - if (cell->line(e)->has_children()) - // the only place where this vertex could have been - // hiding is on the mid-edge point of the edge we - // are looking at - if (cell->line(e)->child(0)->vertex_index(1) == vertex) - { - adjacent_cells.insert(cell); - - // jump out of this tangle of nested loops - goto next_cell; - } - } - - // in more than 3d we would probably have to do the same as - // above also for even lower-dimensional objects - Assert (dim <= 3, ExcNotImplemented()); - - // move on to the next cell if we have found the - // vertex on the current one -next_cell: - ; - } - - // if this was an active vertex then there needs to have been - // at least one cell to which it is adjacent! - Assert (adjacent_cells.size() > 0, ExcInternalError()); - - // return the result as a vector, rather than the set we built above - return - std::vector >::type> - (adjacent_cells.begin(), adjacent_cells.end()); - } - - - - namespace - { - template class MeshType, int spacedim> - void find_active_cell_around_point_internal - (const MeshType &mesh, -#ifndef _MSC_VER - std::set::active_cell_iterator> &searched_cells, - std::set::active_cell_iterator> &adjacent_cells) -#else - std::set >::type> &searched_cells, - std::set >::type> &adjacent_cells) -#endif - { -#ifndef _MSC_VER - typedef typename MeshType::active_cell_iterator cell_iterator; -#else - typedef typename dealii::internal::ActiveCellIterator >::type cell_iterator; -#endif - - // update the searched cells - searched_cells.insert(adjacent_cells.begin(), adjacent_cells.end()); - // now we to collect all neighbors - // of the cells in adjacent_cells we - // have not yet searched. - std::set adjacent_cells_new; - - typename std::set::const_iterator - cell = adjacent_cells.begin(), - endc = adjacent_cells.end(); - for (; cell != endc; ++cell) - { - std::vector active_neighbors; - get_active_neighbors >(*cell, active_neighbors); - for (unsigned int i=0; i class MeshType, int spacedim> -#ifndef _MSC_VER - typename MeshType::active_cell_iterator -#else - typename dealii::internal::ActiveCellIterator >::type -#endif - find_active_cell_around_point (const MeshType &mesh, - const Point &p, - const std::vector &marked_vertices) - { - return - find_active_cell_around_point - (StaticMappingQ1::mapping, - mesh, p, marked_vertices).first; - } - - - template class MeshType, int spacedim> -#ifndef _MSC_VER - std::pair::active_cell_iterator, Point > -#else - std::pair >::type, Point > -#endif - find_active_cell_around_point (const Mapping &mapping, - const MeshType &mesh, - const Point &p, - const std::vector &marked_vertices) - { - typedef typename dealii::internal::ActiveCellIterator >::type active_cell_iterator; - - // The best distance is set to the - // maximum allowable distance from - // the unit cell; we assume a - // max. deviation of 1e-10 - double best_distance = 1e-10; - int best_level = -1; - std::pair > best_cell; - - // Find closest vertex and determine - // all adjacent cells - std::vector adjacent_cells_tmp - = find_cells_adjacent_to_vertex(mesh, - find_closest_vertex(mapping, mesh, p, marked_vertices)); - - // Make sure that we have found - // at least one cell adjacent to vertex. - Assert(adjacent_cells_tmp.size()>0, ExcInternalError()); - - // Copy all the cells into a std::set - std::set adjacent_cells (adjacent_cells_tmp.begin(), - adjacent_cells_tmp.end()); - std::set searched_cells; - - // Determine the maximal number of cells - // in the grid. - // As long as we have not found - // the cell and have not searched - // every cell in the triangulation, - // we keep on looking. - const unsigned int n_active_cells = mesh.get_triangulation().n_active_cells(); - bool found = false; - unsigned int cells_searched = 0; - while (!found && cells_searched < n_active_cells) - { - typename std::set::const_iterator - cell = adjacent_cells.begin(), - endc = adjacent_cells.end(); - for (; cell != endc; ++cell) - { - try - { - const Point p_cell = mapping.transform_real_to_unit_cell(*cell, p); - - // calculate the infinity norm of - // the distance vector to the unit cell. - const double dist = GeometryInfo::distance_to_unit_cell(p_cell); - - // We compare if the point is inside the - // unit cell (or at least not too far - // outside). If it is, it is also checked - // that the cell has a more refined state - if ((dist < best_distance) - || - ((dist == best_distance) - && - ((*cell)->level() > best_level))) - { - found = true; - best_distance = dist; - best_level = (*cell)->level(); - best_cell = std::make_pair(*cell, p_cell); - } - } - catch (typename MappingQGeneric::ExcTransformationFailed &) - { - // ok, the transformation - // failed presumably - // because the point we - // are looking for lies - // outside the current - // cell. this means that - // the current cell can't - // be the cell around the - // point, so just ignore - // this cell and move on - // to the next - } - } - - // update the number of cells searched - cells_searched += adjacent_cells.size(); - - // if the user provided a custom mask for vertices, - // terminate the search without trying to expand the search - // to all cells of the triangulation, as done below. - if (marked_vertices.size() > 0) - cells_searched = n_active_cells; - - // if we have not found the cell in - // question and have not yet searched every - // cell, we expand our search to - // all the not already searched neighbors of - // the cells in adjacent_cells. This is - // what find_active_cell_around_point_internal - // is for. - if (!found && cells_searched < n_active_cells) - { - find_active_cell_around_point_internal - (mesh, searched_cells, adjacent_cells); - } - } - - AssertThrow (best_cell.first.state() == IteratorState::valid, - ExcPointNotFound(p)); - - return best_cell; - } - template std::vector > > vertex_to_cell_centers_directions(const Triangulation &mesh, @@ -1735,465 +1263,6 @@ next_cell: - template - std::pair::active_cell_iterator, Point > - find_active_cell_around_point (const hp::MappingCollection &mapping, - const hp::DoFHandler &mesh, - const Point &p) - { - Assert ((mapping.size() == 1) || - (mapping.size() == mesh.get_fe_collection().size()), - ExcMessage ("Mapping collection needs to have either size 1 " - "or size equal to the number of elements in " - "the FECollection.")); - - typedef typename hp::DoFHandler::active_cell_iterator cell_iterator; - - std::pair > best_cell; - //If we have only one element in the MappingCollection, - //we use find_active_cell_around_point using only one - //mapping. - if (mapping.size() == 1) - best_cell = find_active_cell_around_point(mapping[0], mesh, p); - else - { - - - // The best distance is set to the - // maximum allowable distance from - // the unit cell; we assume a - // max. deviation of 1e-10 - double best_distance = 1e-10; - int best_level = -1; - - - // Find closest vertex and determine - // all adjacent cells - unsigned int vertex = find_closest_vertex(mesh, p); - - std::vector adjacent_cells_tmp = - find_cells_adjacent_to_vertex(mesh, vertex); - - // Make sure that we have found - // at least one cell adjacent to vertex. - Assert(adjacent_cells_tmp.size()>0, ExcInternalError()); - - // Copy all the cells into a std::set - std::set adjacent_cells(adjacent_cells_tmp.begin(), adjacent_cells_tmp.end()); - std::set searched_cells; - - // Determine the maximal number of cells - // in the grid. - // As long as we have not found - // the cell and have not searched - // every cell in the triangulation, - // we keep on looking. - const unsigned int n_cells = mesh.get_triangulation().n_cells(); - bool found = false; - unsigned int cells_searched = 0; - while (!found && cells_searched < n_cells) - { - typename std::set::const_iterator - cell = adjacent_cells.begin(), - endc = adjacent_cells.end(); - for (; cell != endc; ++cell) - { - try - { - const Point p_cell = mapping[(*cell)->active_fe_index()].transform_real_to_unit_cell(*cell, p); - - - // calculate the infinity norm of - // the distance vector to the unit cell. - const double dist = GeometryInfo::distance_to_unit_cell(p_cell); - - // We compare if the point is inside the - // unit cell (or at least not too far - // outside). If it is, it is also checked - // that the cell has a more refined state - if (dist < best_distance || - (dist == best_distance && (*cell)->level() > best_level)) - { - found = true; - best_distance = dist; - best_level = (*cell)->level(); - best_cell = std::make_pair(*cell, p_cell); - } - } - catch (typename MappingQGeneric::ExcTransformationFailed &) - { - // ok, the transformation - // failed presumably - // because the point we - // are looking for lies - // outside the current - // cell. this means that - // the current cell can't - // be the cell around the - // point, so just ignore - // this cell and move on - // to the next - } - } - //udpate the number of cells searched - cells_searched += adjacent_cells.size(); - // if we have not found the cell in - // question and have not yet searched every - // cell, we expand our search to - // all the not already searched neighbors of - // the cells in adjacent_cells. - if (!found && cells_searched < n_cells) - { - find_active_cell_around_point_internal - (mesh, searched_cells, adjacent_cells); - } - - } - } - - AssertThrow (best_cell.first.state() == IteratorState::valid, - ExcPointNotFound(p)); - - return best_cell; - } - - - namespace - { - - template - bool - contains_locally_owned_cells (const std::vector &cells) - { - for (typename std::vector::const_iterator - it = cells.begin(); it != cells.end(); ++it) - { - if ((*it)->is_locally_owned()) - return true; - } - return false; - } - - template - bool - contains_artificial_cells (const std::vector &cells) - { - for (typename std::vector::const_iterator - it = cells.begin(); it != cells.end(); ++it) - { - if ((*it)->is_artificial()) - return true; - } - return false; - } - - } - - - - template - std::vector - compute_active_cell_halo_layer - (const MeshType &mesh, - const std::function &predicate) - { - std::vector active_halo_layer; - std::vector locally_active_vertices_on_subdomain (mesh.get_triangulation().n_vertices(), - false); - - // Find the cells for which the predicate is true - // These are the cells around which we wish to construct - // the halo layer - for (typename MeshType::active_cell_iterator - cell = mesh.begin_active(); - cell != mesh.end(); ++cell) - if (predicate(cell)) // True predicate --> Part of subdomain - for (unsigned int v=0; v::vertices_per_cell; ++v) - locally_active_vertices_on_subdomain[cell->vertex_index(v)] = true; - - // Find the cells that do not conform to the predicate - // but share a vertex with the selected subdomain - // These comprise the halo layer - for (typename MeshType::active_cell_iterator - cell = mesh.begin_active(); - cell != mesh.end(); ++cell) - if (!predicate(cell)) // False predicate --> Potential halo cell - for (unsigned int v=0; v::vertices_per_cell; ++v) - if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] == true) - { - active_halo_layer.push_back(cell); - break; - } - - return active_halo_layer; - } - - - - template - std::vector - compute_cell_halo_layer_on_level - (const MeshType &mesh, - const std::function &predicate, - const unsigned int level) - { - std::vector level_halo_layer; - std::vector locally_active_vertices_on_level_subdomain (mesh.get_triangulation().n_vertices(), - false); - - // Find the cells for which the predicate is true - // These are the cells around which we wish to construct - // the halo layer - for (typename MeshType::cell_iterator - cell = mesh.begin(level); - cell != mesh.end(level); ++cell) - if (predicate(cell)) // True predicate --> Part of subdomain - for (unsigned int v=0; v::vertices_per_cell; ++v) - locally_active_vertices_on_level_subdomain[cell->vertex_index(v)] = true; - - // Find the cells that do not conform to the predicate - // but share a vertex with the selected subdomain on that level - // These comprise the halo layer - for (typename MeshType::cell_iterator - cell = mesh.begin(level); - cell != mesh.end(level); ++cell) - if (!predicate(cell)) // False predicate --> Potential halo cell - for (unsigned int v=0; v::vertices_per_cell; ++v) - if (locally_active_vertices_on_level_subdomain[cell->vertex_index(v)] == true) - { - level_halo_layer.push_back(cell); - break; - } - - return level_halo_layer; - } - - - - template - std::vector - compute_ghost_cell_halo_layer (const MeshType &mesh) - { - std::function predicate - = IteratorFilters::LocallyOwnedCell(); - - const std::vector - active_halo_layer = compute_active_cell_halo_layer (mesh, predicate); - - // Check that we never return locally owned or artificial cells - // What is left should only be the ghost cells - Assert(contains_locally_owned_cells(active_halo_layer) == false, - ExcMessage("Halo layer contains locally owned cells")); - Assert(contains_artificial_cells(active_halo_layer) == false, - ExcMessage("Halo layer contains artificial cells")); - - return active_halo_layer; - } - - - - template - std::vector - compute_active_cell_layer_within_distance - (const MeshType &mesh, - const std::function &predicate, - const double layer_thickness) - { - std::vector subdomain_boundary_cells, active_cell_layer_within_distance; - std::vector vertices_outside_subdomain ( mesh.get_triangulation().n_vertices(), - false); - - const unsigned int spacedim = MeshType::space_dimension; - - unsigned int n_non_predicate_cells = 0; // Number of non predicate cells - - // Find the layer of cells for which predicate is true and that - // are on the boundary with other cells. These are - // subdomain boundary cells. - - // Find the cells for which the predicate is false - // These are the cells which are around the predicate subdomain - for ( typename MeshType::active_cell_iterator - cell = mesh.begin_active(); - cell != mesh.end(); ++cell) - if ( !predicate(cell)) // Negation of predicate --> Not Part of subdomain - { - for (unsigned int v=0; v::vertices_per_cell; ++v) - vertices_outside_subdomain[cell->vertex_index(v)] = true; - n_non_predicate_cells++; - } - - // If all the active cells conform to the predicate - // or if none of the active cells conform to the predicate - // there is no active cell layer around the predicate - // subdomain (within any distance) - if ( n_non_predicate_cells == 0 || n_non_predicate_cells == mesh.get_triangulation().n_active_cells() ) - return std::vector(); - - // Find the cells that conform to the predicate - // but share a vertex with the cell not in the predicate subdomain - for ( typename MeshType::active_cell_iterator - cell = mesh.begin_active(); - cell != mesh.end(); ++cell) - if ( predicate(cell)) // True predicate --> Potential boundary cell of the subdomain - for (unsigned int v=0; v::vertices_per_cell; ++v) - if (vertices_outside_subdomain[cell->vertex_index(v)] == true) - { - subdomain_boundary_cells.push_back(cell); - break; // No need to go through remaining vertices - } - - // To cheaply filter out some cells located far away from the predicate subdomain, - // get the bounding box of the predicate subdomain. - std::pair< Point, Point > bounding_box = compute_bounding_box( mesh, - predicate ); - - // DOUBLE_EPSILON to compare really close double values - const double &DOUBLE_EPSILON = 100.*std::numeric_limits::epsilon(); - - // Add layer_thickness to the bounding box - for ( unsigned int d=0; d > subdomain_boundary_cells_centers; // cache all the subdomain boundary cells centers here - std::vector subdomain_boundary_cells_radii; // cache all the subdomain boundary cells radii - subdomain_boundary_cells_centers.reserve (subdomain_boundary_cells.size()); - subdomain_boundary_cells_radii.reserve (subdomain_boundary_cells.size()); - // compute cell radius for each boundary cell of the predicate subdomain - for ( typename std::vector::const_iterator - subdomain_boundary_cell_iterator = subdomain_boundary_cells.begin(); - subdomain_boundary_cell_iterator != subdomain_boundary_cells.end(); ++subdomain_boundary_cell_iterator ) - { - const std::pair, double> & - subdomain_boundary_cell_enclosing_ball = (*subdomain_boundary_cell_iterator)->enclosing_ball(); - - subdomain_boundary_cells_centers.push_back( subdomain_boundary_cell_enclosing_ball.first); - subdomain_boundary_cells_radii.push_back( subdomain_boundary_cell_enclosing_ball.second); - } - AssertThrow( subdomain_boundary_cells_radii.size() == subdomain_boundary_cells_centers.size(), - ExcInternalError()); - - // Find the cells that are within layer_thickness of predicate subdomain boundary - // distance but are inside the extended bounding box. - // Most cells might be outside the extended bounding box, so we could skip them. - // Those cells that are inside the extended bounding box but are not part of the - // predicate subdomain are possible candidates to be within the distance to the - // boundary cells of the predicate subdomain. - for ( typename MeshType::active_cell_iterator - cell = mesh.begin_active(); - cell != mesh.end(); ++cell) - { - // Ignore all the cells that are in the predicate subdomain - if ( predicate(cell)) - continue; - - const std::pair, double> &cell_enclosing_ball - = cell->enclosing_ball(); - - const Point &cell_enclosing_ball_center = cell_enclosing_ball.first; - const double &cell_enclosing_ball_radius = cell_enclosing_ball.second; - - bool cell_inside = true; // reset for each cell - - for (unsigned int d = 0; d < spacedim; ++d) - cell_inside &= (cell_enclosing_ball_center[d] + cell_enclosing_ball_radius > bounding_box.first[d]) - && (cell_enclosing_ball_center[d] - cell_enclosing_ball_radius < bounding_box.second[d]); - // cell_inside is true if its enclosing ball intersects the extended bounding box - - // Ignore all the cells that are outside the extended bounding box - if (cell_inside) - for (unsigned int i =0; i< subdomain_boundary_cells_radii.size(); ++i) - if ( cell_enclosing_ball_center.distance_square(subdomain_boundary_cells_centers[i]) - < Utilities::fixed_power<2>( cell_enclosing_ball_radius + - subdomain_boundary_cells_radii[i] + - layer_thickness + DOUBLE_EPSILON )) - { - active_cell_layer_within_distance.push_back(cell); - break; // Exit the loop checking all the remaining subdomain boundary cells - } - - } - return active_cell_layer_within_distance; - } - - - - template - std::vector - compute_ghost_cell_layer_within_distance ( const MeshType &mesh, const double layer_thickness) - { - IteratorFilters::LocallyOwnedCell locally_owned_cell_predicate; - std::function predicate (locally_owned_cell_predicate); - - const std::vector - ghost_cell_layer_within_distance = compute_active_cell_layer_within_distance (mesh, predicate, layer_thickness); - - // Check that we never return locally owned or artificial cells - // What is left should only be the ghost cells - Assert(contains_locally_owned_cells(ghost_cell_layer_within_distance) == false, - ExcMessage("Ghost cells within layer_thickness contains locally owned cells.")); - Assert(contains_artificial_cells(ghost_cell_layer_within_distance) == false, - ExcMessage("Ghost cells within layer_thickness contains artificial cells." - "The function compute_ghost_cell_layer_within_distance " - "is probably called while using parallel::distributed::Triangulation. " - "In such case please refer to the description of this function.")); - - return ghost_cell_layer_within_distance; - } - - - - template < class MeshType> - std::pair< Point, Point > - compute_bounding_box - ( const MeshType &mesh, - const std::function &predicate ) - { - std::vector locally_active_vertices_on_subdomain (mesh.get_triangulation().n_vertices(), - false); - - const unsigned int spacedim = MeshType::space_dimension; - - // Two extreme points can define the bounding box - // around the active cells that conform to the given predicate. - Point maxp, minp; - - // initialize minp and maxp with the first predicate cell center - for ( typename MeshType::active_cell_iterator - cell = mesh.begin_active(); - cell != mesh.end(); ++cell) - if ( predicate(cell)) - { - minp = cell->center(); - maxp = cell->center(); - break; - } - - // Run through all the cells to check if it belongs to predicate domain, - // if it belongs to the predicate domain, extend the bounding box. - for ( typename MeshType::active_cell_iterator - cell = mesh.begin_active(); - cell != mesh.end(); ++cell) - if (predicate(cell)) // True predicate --> Part of subdomain - for (unsigned int v=0; v::vertices_per_cell; ++v) - if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] == false) - { - locally_active_vertices_on_subdomain[cell->vertex_index(v)] = true; - for ( unsigned int d=0; dvertex(v)[d]); - maxp[d] = std::max( maxp[d], cell->vertex(v)[d]); - } - } - - return std::make_pair( minp, maxp ); - } - - - namespace internal { namespace BoundingBoxPredicate @@ -3150,143 +2219,9 @@ next_cell: - template - std::list > - get_finest_common_cells (const MeshType &mesh_1, - const MeshType &mesh_2) - { - Assert (have_same_coarse_mesh (mesh_1, mesh_2), - ExcMessage ("The two meshes must be represent triangulations that " - "have the same coarse meshes")); - - // the algorithm goes as follows: - // first, we fill a list with pairs - // of iterators common to the two - // meshes on the coarsest - // level. then we traverse the - // list; each time, we find a pair - // of iterators for which both - // correspond to non-active cells, - // we delete this item and push the - // pairs of iterators to their - // children to the back. if these - // again both correspond to - // non-active cells, we will get to - // the later on for further - // consideration - typedef - std::list > - CellList; - - CellList cell_list; - - // first push the coarse level cells - typename MeshType::cell_iterator - cell_1 = mesh_1.begin(0), - cell_2 = mesh_2.begin(0); - for (; cell_1 != mesh_1.end(0); ++cell_1, ++cell_2) - cell_list.emplace_back (cell_1, cell_2); - - // then traverse list as described - // above - typename CellList::iterator cell_pair = cell_list.begin(); - while (cell_pair != cell_list.end()) - { - // if both cells in this pair - // have children, then erase - // this element and push their - // children instead - if (cell_pair->first->has_children() - && - cell_pair->second->has_children()) - { - Assert(cell_pair->first->refinement_case()== - cell_pair->second->refinement_case(), ExcNotImplemented()); - for (unsigned int c=0; cfirst->n_children(); ++c) - cell_list.emplace_back (cell_pair->first->child(c), - cell_pair->second->child(c)); - - // erasing an iterator - // keeps other iterators - // valid, so already - // advance the present - // iterator by one and then - // delete the element we've - // visited before - const typename CellList::iterator previous_cell_pair = cell_pair; - ++cell_pair; - - cell_list.erase (previous_cell_pair); - } - else - // both cells are active, do - // nothing - ++cell_pair; - } - - // just to make sure everything is ok, - // validate that all pairs have at least one - // active iterator or have different - // refinement_cases - for (cell_pair = cell_list.begin(); cell_pair != cell_list.end(); ++cell_pair) - Assert (cell_pair->first->active() - || - cell_pair->second->active() - || - (cell_pair->first->refinement_case() - != cell_pair->second->refinement_case()), - ExcInternalError()); - - return cell_list; - } - - template - bool - have_same_coarse_mesh (const Triangulation &mesh_1, - const Triangulation &mesh_2) - { - // make sure the two meshes have - // the same number of coarse cells - if (mesh_1.n_cells (0) != mesh_2.n_cells (0)) - return false; - - // if so, also make sure they have - // the same vertices on the cells - // of the coarse mesh - typename Triangulation::cell_iterator - cell_1 = mesh_1.begin(0), - cell_2 = mesh_2.begin(0), - endc = mesh_1.end(0); - for (; cell_1!=endc; ++cell_1, ++cell_2) - for (unsigned int v=0; v::vertices_per_cell; ++v) - if (cell_1->vertex(v) != cell_2->vertex(v)) - return false; - - // if we've gotten through all - // this, then the meshes really - // seem to have a common coarse - // mesh - return true; - } - - - - template - bool - have_same_coarse_mesh (const MeshType &mesh_1, - const MeshType &mesh_2) - { - return have_same_coarse_mesh (mesh_1.get_triangulation(), - mesh_2.get_triangulation()); - } - - - - template - double - minimal_cell_diameter (const Triangulation &triangulation) + template + double + minimal_cell_diameter (const Triangulation &triangulation) { double min_diameter = triangulation.begin_active()->diameter(); for (typename Triangulation::active_cell_iterator @@ -3790,920 +2725,6 @@ next_cell: - template - std::vector - get_patch_around_cell(const typename MeshType::active_cell_iterator &cell) - { - Assert (cell->is_locally_owned(), - ExcMessage ("This function only makes sense if the cell for " - "which you are asking for a patch, is locally " - "owned.")); - - std::vector patch; - patch.push_back (cell); - for (unsigned int face_number=0; face_number::faces_per_cell; ++face_number) - if (cell->face(face_number)->at_boundary()==false) - { - if (cell->neighbor(face_number)->has_children() == false) - patch.push_back (cell->neighbor(face_number)); - else - // the neighbor is refined. in 2d/3d, we can simply ask for the children - // of the neighbor because they can not be further refined and, - // consequently, the children is active - if (MeshType::dimension > 1) - { - for (unsigned int subface=0; subfaceface(face_number)->n_children(); ++subface) - patch.push_back (cell->neighbor_child_on_subface (face_number, subface)); - } - else - { - // in 1d, we need to work a bit harder: iterate until we find - // the child by going from cell to child to child etc - typename MeshType::cell_iterator neighbor - = cell->neighbor (face_number); - while (neighbor->has_children()) - neighbor = neighbor->child(1-face_number); - - Assert (neighbor->neighbor(1-face_number) == cell, ExcInternalError()); - patch.push_back (neighbor); - } - } - return patch; - } - - - - template - std::vector - get_cells_at_coarsest_common_level ( - const std::vector &patch) - { - Assert (patch.size() > 0, ExcMessage("Vector containing patch cells should not be an empty vector!")); - // In order to extract the set of cells with the coarsest common level from the give vector of cells: - // First it finds the number associated with the minimum level of refinmenet, namely "min_level" - int min_level = patch[0]->level(); - - for (unsigned int i=0; ilevel() ); - std::set uniform_cells; - typename std::vector::const_iterator patch_cell; - // it loops through all cells of the input vector - for (patch_cell=patch.begin(); patch_cell!=patch.end () ; ++patch_cell) - { - // If the refinement level of each cell i the loop be equal to the min_level, so that - // that cell inserted into the set of uniform_cells, as the set of cells with the coarsest common refinement level - if ((*patch_cell)->level() == min_level) - uniform_cells.insert (*patch_cell); - else - // If not, it asks for the parent of the cell, until it finds the parent cell - // with the refinement level equal to the min_level and inserts that parent cell into the - // the set of uniform_cells, as the set of cells with the coarsest common refinement level. - { - typename Container::cell_iterator parent = *patch_cell; - - while (parent->level() > min_level) - parent = parent-> parent(); - uniform_cells.insert (parent); - } - } - - return std::vector (uniform_cells.begin(), - uniform_cells.end()); - } - - - - template - void build_triangulation_from_patch(const std::vector &patch, - Triangulation &local_triangulation, - std::map::active_cell_iterator, - typename Container::active_cell_iterator> &patch_to_global_tria_map) - - { - const std::vector uniform_cells = - get_cells_at_coarsest_common_level (patch); - // First it creates triangulation from the vector of "uniform_cells" - local_triangulation.clear(); - std::vector > vertices; - const unsigned int n_uniform_cells=uniform_cells.size(); - std::vector > cells(n_uniform_cells); - unsigned int k=0;// for enumerating cells - unsigned int i=0;// for enumerating vertices - typename std::vector::const_iterator uniform_cell; - for (uniform_cell=uniform_cells.begin(); uniform_cell!=uniform_cells.end(); ++uniform_cell) - { - for (unsigned int v=0; v::vertices_per_cell; ++v) - { - Point position=(*uniform_cell)->vertex (v); - bool repeat_vertex=false; - - for (unsigned int m=0; m::cell_iterator, - typename Container::cell_iterator> patch_to_global_tria_map_tmp; - for (typename Triangulation::cell_iterator coarse_cell = local_triangulation.begin(); - coarse_cell != local_triangulation.end(); ++coarse_cell, ++index) - { - patch_to_global_tria_map_tmp.insert (std::make_pair(coarse_cell, uniform_cells[index])); - // To ensure that the cells with the same coordinates (here, we compare their centers) are mapped into each other. - - Assert(coarse_cell->center().distance( uniform_cells[index]->center())<=1e-15*coarse_cell->diameter(), - ExcInternalError()); - } - bool refinement_necessary; - // In this loop we start to do refinement on the above coarse triangulation to reach - // to the same level of refinement as the patch cells are really on - do - { - refinement_necessary = false; - for (typename Triangulation::active_cell_iterator - active_tria_cell = local_triangulation.begin_active(); - active_tria_cell != local_triangulation.end(); ++active_tria_cell) - { - if (patch_to_global_tria_map_tmp[active_tria_cell]->has_children()) - { - active_tria_cell -> set_refine_flag(); - refinement_necessary = true; - } - else for (unsigned int i=0; i::vertices_per_cell; ++v) - active_tria_cell->vertex(v) = patch[i]->vertex(v); - - Assert(active_tria_cell->center().distance(patch_to_global_tria_map_tmp[active_tria_cell]->center()) - <=1e-15*active_tria_cell->diameter(), ExcInternalError()); - - active_tria_cell->set_user_flag(); - break; - } - } - } - - if (refinement_necessary) - { - local_triangulation.execute_coarsening_and_refinement (); - - for (typename Triangulation::cell_iterator - cell = local_triangulation.begin(); - cell != local_triangulation.end(); ++cell) - { - - if (patch_to_global_tria_map_tmp.find(cell)!=patch_to_global_tria_map_tmp.end()) - { - if (cell-> has_children()) - { - // Note: Since the cell got children, then it should not be in the map anymore - // children may be added into the map, instead - - // these children may not yet be in the map - for (unsigned int c=0; cn_children(); ++c) - { - if (patch_to_global_tria_map_tmp.find(cell->child(c)) == - patch_to_global_tria_map_tmp.end()) - { - patch_to_global_tria_map_tmp.insert (std::make_pair(cell->child(c), - patch_to_global_tria_map_tmp[cell]->child(c))); - - // One might be tempted to assert that the cell - // being added here has the same center as the - // equivalent cell in the global triangulation, - // but it may not be the case. For triangulations - // that have been perturbed or smoothed, the cell - // indices and levels may be the same, but the - // vertex locations may not. We adjust - // the vertices of the cells that have no - // children (ie the active cells) to be - // consistent with the global triangulation - // later on and add assertions at that time - // to guarantee the cells in the - // local_triangulation are physically at the same - // locations of the cells in the patch of the - // global triangulation. - - } - } - // The parent cell whose children were added - // into the map should be deleted from the map - patch_to_global_tria_map_tmp.erase(cell); - } - } - } - } - - } - while (refinement_necessary); - - - // Last assertion check to make sure we have the right cells and centers - // in the map, and hence the correct vertices of the triangulation - for (typename Triangulation::cell_iterator - cell = local_triangulation.begin(); - cell != local_triangulation.end(); ++cell) - { - if (cell->user_flag_set() ) - { - Assert(patch_to_global_tria_map_tmp.find(cell) != patch_to_global_tria_map_tmp.end(), - ExcInternalError() ); - - Assert(cell->center().distance( patch_to_global_tria_map_tmp[cell]->center())<=1e-15*cell->diameter(), - ExcInternalError()); - } - } - - - typename std::map::cell_iterator, - typename Container::cell_iterator>::iterator map_tmp_it = - patch_to_global_tria_map_tmp.begin(),map_tmp_end = patch_to_global_tria_map_tmp.end(); - // Now we just need to take the temporary map of pairs of type cell_iterator "patch_to_global_tria_map_tmp" - // making pair of active_cell_iterators so that filling out the final map "patch_to_global_tria_map" - for (; map_tmp_it!=map_tmp_end; ++map_tmp_it) - patch_to_global_tria_map[map_tmp_it->first] = map_tmp_it->second; - } - - - - - template - std::map< types::global_dof_index,std::vector > - get_dof_to_support_patch_map(DoFHandlerType &dof_handler) - { - - // This is the map from global_dof_index to - // a set of cells on patch. We first map into - // a set because it is very likely that we - // will attempt to add a cell more than once - // to a particular patch and we want to preserve - // uniqueness of cell iterators. std::set does this - // automatically for us. Later after it is all - // constructed, we will copy to a map of vectors - // since that is the prefered output for other - // functions. - std::map< types::global_dof_index,std::set > dof_to_set_of_cells_map; - - std::vector local_dof_indices; - std::vector local_face_dof_indices; - std::vector local_line_dof_indices; - - // a place to save the dof_handler user flags and restore them later - // to maintain const of dof_handler. - std::vector user_flags; - - - // in 3d, we need pointers from active lines to the - // active parent lines, so we construct it as needed. - std::map lines_to_parent_lines_map; - if (DoFHandlerType::dimension == 3) - { - - // save user flags as they will be modified and then later restored - dof_handler.get_triangulation().save_user_flags(user_flags); - const_cast &>(dof_handler.get_triangulation()).clear_user_flags (); - - - typename DoFHandlerType::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - { - // We only want lines that are locally_relevant - // although it doesn't hurt to have lines that - // are children of ghost cells since there are - // few and we don't have to use them. - if (cell->is_artificial() == false) - { - for (unsigned int l=0; l::lines_per_cell; ++l) - if (cell->line(l)->has_children()) - for (unsigned int c=0; cline(l)->n_children(); ++c) - { - lines_to_parent_lines_map[cell->line(l)->child(c)] = cell->line(l); - // set flags to know that child - // line has an active parent. - cell->line(l)->child(c)->set_user_flag(); - } - } - } - } - - - // We loop through all cells and add cell to the - // map for the dofs that it immediately touches - // and then account for all the other dofs of - // which it is a part, mainly the ones that must - // be added on account of adaptivity hanging node - // constraints. - typename DoFHandlerType::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - { - // Need to loop through all cells that could - // be in the patch of dofs on locally_owned - // cells including ghost cells - if (cell->is_artificial() == false) - { - const unsigned int n_dofs_per_cell = cell->get_fe().dofs_per_cell; - local_dof_indices.resize(n_dofs_per_cell); - - // Take care of adding cell pointer to each - // dofs that exists on cell. - cell->get_dof_indices(local_dof_indices); - for (unsigned int i=0; i< n_dofs_per_cell; ++i ) - dof_to_set_of_cells_map[local_dof_indices[i]].insert(cell); - - // In the case of the adjacent cell (over - // faces or edges) being more refined, we - // want to add all of the children to the - // patch since the support function at that - // dof could be non-zero along that entire - // face (or line). - - // Take care of dofs on neighbor faces - for (unsigned int f=0; f::faces_per_cell; ++f) - { - if (cell->face(f)->has_children()) - { - for (unsigned int c=0; cface(f)->n_children(); ++c) - { - // Add cell to dofs of all subfaces - // - // *-------------------*----------*---------* - // | | add cell | | - // | |<- to dofs| | - // | |of subface| | - // | cell *----------*---------* - // | | add cell | | - // | |<- to dofs| | - // | |of subface| | - // *-------------------*----------*---------* - // - Assert (cell->face(f)->child(c)->has_children() == false, ExcInternalError()); - - const unsigned int n_dofs_per_face = cell->get_fe().dofs_per_face; - local_face_dof_indices.resize(n_dofs_per_face); - - cell->face(f)->child(c)->get_dof_indices(local_face_dof_indices); - for (unsigned int i=0; i< n_dofs_per_face; ++i ) - dof_to_set_of_cells_map[local_face_dof_indices[i]].insert(cell); - } - } - else if ((cell->face(f)->at_boundary() == false) && (cell->neighbor_is_coarser(f))) - { - - // Add cell to dofs of parent face and all - // child faces of parent face - // - // *-------------------*----------*---------* - // | | | | - // | | cell | | - // | add cell | | | - // | to dofs -> *----------*---------* - // | of parent | add cell | | - // | face |<- to dofs| | - // | |of subface| | - // *-------------------*----------*---------* - // - - // Add cell to all dofs of parent face - std::pair neighbor_face_no_subface_no = cell->neighbor_of_coarser_neighbor(f); - unsigned int face_no = neighbor_face_no_subface_no.first; - unsigned int subface = neighbor_face_no_subface_no.second; - - const unsigned int n_dofs_per_face = cell->get_fe().dofs_per_face; - local_face_dof_indices.resize(n_dofs_per_face); - - cell->neighbor(f)->face(face_no)->get_dof_indices(local_face_dof_indices); - for (unsigned int i=0; i< n_dofs_per_face; ++i ) - dof_to_set_of_cells_map[local_face_dof_indices[i]].insert(cell); - - // Add cell to all dofs of children of - // parent face - for (unsigned int c=0; cneighbor(f)->face(face_no)->n_children(); ++c) - { - if (c != subface) // don't repeat work on dofs of original cell - { - const unsigned int n_dofs_per_face = cell->get_fe().dofs_per_face; - local_face_dof_indices.resize(n_dofs_per_face); - - Assert (cell->neighbor(f)->face(face_no)->child(c)->has_children() == false, ExcInternalError()); - cell->neighbor(f)->face(face_no)->child(c)->get_dof_indices(local_face_dof_indices); - for (unsigned int i=0; i::lines_per_cell; ++l) - { - if (cell->line(l)->has_children()) - { - for (unsigned int c=0; cline(l)->n_children(); ++c) - { - Assert (cell->line(l)->child(c)->has_children() == false, ExcInternalError()); - - // dofs_per_line returns number of dofs - // on line not including the vertices of the line. - const unsigned int n_dofs_per_line = 2*cell->get_fe().dofs_per_vertex - + cell->get_fe().dofs_per_line; - local_line_dof_indices.resize(n_dofs_per_line); - - cell->line(l)->child(c)->get_dof_indices(local_line_dof_indices); - for (unsigned int i=0; iline(l)->user_flag_set() == true) - { - typename DoFHandlerType::line_iterator parent_line = lines_to_parent_lines_map[cell->line(l)]; - Assert (parent_line->has_children(), ExcInternalError() ); - - // dofs_per_line returns number of dofs - // on line not including the vertices of the line. - const unsigned int n_dofs_per_line = 2*cell->get_fe().dofs_per_vertex - + cell->get_fe().dofs_per_line; - local_line_dof_indices.resize(n_dofs_per_line); - - parent_line->get_dof_indices(local_line_dof_indices); - for (unsigned int i=0; in_children(); ++c) - { - Assert (parent_line->child(c)->has_children() == false, ExcInternalError()); - - const unsigned int n_dofs_per_line = 2*cell->get_fe().dofs_per_vertex - + cell->get_fe().dofs_per_line; - local_line_dof_indices.resize(n_dofs_per_line); - - parent_line->child(c)->get_dof_indices(local_line_dof_indices); - for (unsigned int i=0; iis_locally_owned() - }// for cells - - - if (DoFHandlerType::dimension == 3) - { - // finally, restore user flags that were changed above - // to when we constructed the pointers to parent of lines - // Since dof_handler is const, we must leave it unchanged. - const_cast &>(dof_handler.get_triangulation()).load_user_flags (user_flags); - } - - // Finally, we copy map of sets to - // map of vectors using the std::vector::assign() function - std::map< types::global_dof_index, std::vector > dof_to_cell_patches; - - typename std::map >::iterator - it = dof_to_set_of_cells_map.begin(), - it_end = dof_to_set_of_cells_map.end(); - for ( ; it!=it_end; ++it) - dof_to_cell_patches[it->first].assign( it->second.begin(), it->second.end() ); - - return dof_to_cell_patches; - } - - - - /* - * Internally used in orthogonal_equality - * - * An orthogonal equality test for points: - * - * point1 and point2 are considered equal, if - * matrix.point1 + offset - point2 - * is parallel to the unit vector in - */ - template - inline bool orthogonal_equality (const Point &point1, - const Point &point2, - const int direction, - const Tensor<1,spacedim> &offset, - const FullMatrix &matrix) - { - Assert (0<=direction && direction distance; - - if (matrix.m() == spacedim) - for (int i = 0; i < spacedim; ++i) - for (int j = 0; j < spacedim; ++j) - distance(i) += matrix(i,j) * point1(j); - else - distance = point1; - - distance += offset - point2; - - for (int i = 0; i < spacedim; ++i) - { - // Only compare coordinate-components != direction: - if (i == direction) - continue; - - if (fabs(distance(i)) > 1.e-10) - return false; - } - - return true; - } - - - /* - * Internally used in orthogonal_equality - * - * A lookup table to transform vertex matchings to orientation flags of - * the form (face_orientation, face_flip, face_rotation) - * - * See the comment on the next function as well as the detailed - * documentation of make_periodicity_constraints and - * collect_periodic_faces for details - */ - template struct OrientationLookupTable {}; - - template <> struct OrientationLookupTable<1> - { - typedef std::array::vertices_per_face> MATCH_T; - static inline std::bitset<3> lookup (const MATCH_T &) - { - // The 1D case is trivial - return 1; // [true ,false,false] - } - }; - - template <> struct OrientationLookupTable<2> - { - typedef std::array::vertices_per_face> MATCH_T; - static inline std::bitset<3> lookup (const MATCH_T &matching) - { - // In 2D matching faces (=lines) results in two cases: Either - // they are aligned or flipped. We store this "line_flip" - // property somewhat sloppy as "face_flip" - // (always: face_orientation = true, face_rotation = false) - - static const MATCH_T m_tff = {{ 0, 1 }}; - if (matching == m_tff) return 1; // [true ,false,false] - static const MATCH_T m_ttf = {{ 1, 0 }}; - if (matching == m_ttf) return 3; // [true ,true ,false] - Assert(false, ExcInternalError()); - // what follows is dead code, but it avoids warnings about the lack - // of a return value - return 0; - } - }; - - template <> struct OrientationLookupTable<3> - { - typedef std::array::vertices_per_face> MATCH_T; - static inline std::bitset<3> lookup (const MATCH_T &matching) - { - // The full fledged 3D case. *Yay* - // See the documentation in include/deal.II/base/geometry_info.h - // as well as the actual implementation in source/grid/tria.cc - // for more details... - - static const MATCH_T m_tff = {{ 0, 1, 2, 3 }}; - if (matching == m_tff) return 1; // [true ,false,false] - static const MATCH_T m_tft = {{ 1, 3, 0, 2 }}; - if (matching == m_tft) return 5; // [true ,false,true ] - static const MATCH_T m_ttf = {{ 3, 2, 1, 0 }}; - if (matching == m_ttf) return 3; // [true ,true ,false] - static const MATCH_T m_ttt = {{ 2, 0, 3, 1 }}; - if (matching == m_ttt) return 7; // [true ,true ,true ] - static const MATCH_T m_fff = {{ 0, 2, 1, 3 }}; - if (matching == m_fff) return 0; // [false,false,false] - static const MATCH_T m_fft = {{ 2, 3, 0, 1 }}; - if (matching == m_fft) return 4; // [false,false,true ] - static const MATCH_T m_ftf = {{ 3, 1, 2, 0 }}; - if (matching == m_ftf) return 2; // [false,true ,false] - static const MATCH_T m_ftt = {{ 1, 0, 3, 2 }}; - if (matching == m_ftt) return 6; // [false,true ,true ] - Assert(false, ExcInternalError()); - // what follows is dead code, but it avoids warnings about the lack - // of a return value - return 0; - } - }; - - - - template - inline bool - orthogonal_equality (std::bitset<3> &orientation, - const FaceIterator &face1, - const FaceIterator &face2, - const int direction, - const Tensor<1,FaceIterator::AccessorType::space_dimension> &offset, - const FullMatrix &matrix) - { - Assert(matrix.m() == matrix.n(), - ExcMessage("The supplied matrix must be a square matrix")); - - static const int dim = FaceIterator::AccessorType::dimension; - - // Do a full matching of the face vertices: - - std:: - array::vertices_per_face> matching; - - std::set face2_vertices; - for (unsigned int i = 0; i < GeometryInfo::vertices_per_face; ++i) - face2_vertices.insert(i); - - for (unsigned int i = 0; i < GeometryInfo::vertices_per_face; ++i) - { - for (std::set::iterator it = face2_vertices.begin(); - it != face2_vertices.end(); - ++it) - { - if (orthogonal_equality(face1->vertex(i),face2->vertex(*it), - direction, offset, matrix)) - { - matching[i] = *it; - face2_vertices.erase(it); - break; // jump out of the innermost loop - } - } - } - - // And finally, a lookup to determine the ordering bitmask: - if (face2_vertices.empty()) - orientation = OrientationLookupTable::lookup(matching); - - return face2_vertices.empty(); - } - - - - template - inline bool - orthogonal_equality (const FaceIterator &face1, - const FaceIterator &face2, - const int direction, - const Tensor<1,FaceIterator::AccessorType::space_dimension> &offset, - const FullMatrix &matrix) - { - // Call the function above with a dummy orientation array - std::bitset<3> dummy; - return orthogonal_equality (dummy, face1, face2, direction, offset, matrix); - } - - - - /* - * Internally used in collect_periodic_faces - */ - template - void - match_periodic_face_pairs - (std::set > &pairs1, - std::set::type, unsigned int> > &pairs2, - const int direction, - std::vector > &matched_pairs, - const dealii::Tensor<1,CellIterator::AccessorType::space_dimension> &offset, - const FullMatrix &matrix) - { - static const int space_dim = CellIterator::AccessorType::space_dimension; - (void)space_dim; - Assert (0<=direction && direction orientation; - typedef typename std::set - >::const_iterator PairIterator; - for (PairIterator it1 = pairs1.begin(); it1 != pairs1.end(); ++it1) - { - for (PairIterator it2 = pairs2.begin(); it2 != pairs2.end(); ++it2) - { - const CellIterator cell1 = it1->first; - const CellIterator cell2 = it2->first; - const unsigned int face_idx1 = it1->second; - const unsigned int face_idx2 = it2->second; - if (GridTools::orthogonal_equality(orientation, - cell1->face(face_idx1), - cell2->face(face_idx2), - direction, offset, - matrix)) - { - // We have a match, so insert the matching pairs and - // remove the matched cell in pairs2 to speed up the - // matching: - const PeriodicFacePair matched_face = - { - {cell1, cell2}, - {face_idx1, face_idx2}, - orientation, - matrix - }; - matched_pairs.push_back(matched_face); - pairs2.erase(it2); - ++n_matches; - break; - } - } - } - - //Assure that all faces are matched - AssertThrow (n_matches == pairs1.size() && pairs2.size() == 0, - ExcMessage ("Unmatched faces on periodic boundaries")); - } - - - - template - void - collect_periodic_faces - (const MeshType &mesh, - const types::boundary_id b_id1, - const types::boundary_id b_id2, - const int direction, - std::vector > &matched_pairs, - const Tensor<1,MeshType::space_dimension> &offset, - const FullMatrix &matrix) - { - static const int dim = MeshType::dimension; - static const int space_dim = MeshType::space_dimension; - (void)dim; - (void)space_dim; - Assert (0<=direction && direction > pairs1; - std::set > pairs2; - - for (typename MeshType::cell_iterator cell = mesh.begin(0); - cell != mesh.end(0); ++cell) - { - for (unsigned int i = 0; i < GeometryInfo::faces_per_cell; ++i) - { - const typename MeshType::face_iterator face = cell->face(i); - if (face->at_boundary() && face->boundary_id() == b_id1) - { - const std::pair pair1 - = std::make_pair(cell, i); - pairs1.insert(pair1); - } - - if (face->at_boundary() && face->boundary_id() == b_id2) - { - const std::pair pair2 - = std::make_pair(cell, i); - pairs2.insert(pair2); - } - } - } - - Assert (pairs1.size() == pairs2.size(), - ExcMessage ("Unmatched faces on periodic boundaries")); - - Assert (pairs1.size() > 0, - ExcMessage("No new periodic face pairs have been found. " - "Are you sure that you've selected the correct boundary " - "id's and that the coarsest level mesh is colorized?")); - - // and call match_periodic_face_pairs that does the actual matching: - match_periodic_face_pairs(pairs1, pairs2, direction, matched_pairs, offset, - matrix); - } - - - - template - void - collect_periodic_faces - (const MeshType &mesh, - const types::boundary_id b_id, - const int direction, - std::vector > &matched_pairs, - const Tensor<1,MeshType::space_dimension> &offset, - const FullMatrix &matrix) - { - static const int dim = MeshType::dimension; - static const int space_dim = MeshType::space_dimension; - (void)dim; - (void)space_dim; - Assert (0<=direction && direction > pairs1; - std::set > pairs2; - - for (typename MeshType::cell_iterator cell = mesh.begin(0); - cell != mesh.end(0); ++cell) - { - const typename MeshType::face_iterator face_1 = cell->face(2*direction); - const typename MeshType::face_iterator face_2 = cell->face(2*direction+1); - - if (face_1->at_boundary() && face_1->boundary_id() == b_id) - { - const std::pair pair1 - = std::make_pair(cell, 2*direction); - pairs1.insert(pair1); - } - - if (face_2->at_boundary() && face_2->boundary_id() == b_id) - { - const std::pair pair2 - = std::make_pair(cell, 2*direction+1); - pairs2.insert(pair2); - } - } - - Assert (pairs1.size() == pairs2.size(), - ExcMessage ("Unmatched faces on periodic boundaries")); - - Assert (pairs1.size() > 0, - ExcMessage("No new periodic face pairs have been found. " - "Are you sure that you've selected the correct boundary " - "id's and that the coarsest level mesh is colorized?")); - -#ifdef DEBUG - const unsigned int size_old = matched_pairs.size(); -#endif - - // and call match_periodic_face_pairs that does the actual matching: - match_periodic_face_pairs(pairs1, pairs2, direction, matched_pairs, offset, - matrix); - -#ifdef DEBUG - //check for standard orientation - const unsigned int size_new = matched_pairs.size(); - for (unsigned int i = size_old; i < size_new; ++i) - { - Assert(matched_pairs[i].orientation == 1, - ExcMessage("Found a face match with non standard orientation. " - "This function is only suitable for meshes with cells " - "in default orientation")); - } -#endif - } - - - template void copy_boundary_to_manifold_id(Triangulation &tria, const bool reset_boundary_ids) diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index 51d4df9a32..22d8a2afdb 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -39,7 +39,8 @@ for (X : TRIANGULATIONS; deal_II_dimension : DIMENSIONS ; deal_II_space_dimensio } -for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS) +// now also instantiate a few additional functions for parallel::distributed::Triangulation +for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS) { #if deal_II_dimension <= deal_II_space_dimension @@ -47,122 +48,38 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II template unsigned int - find_closest_vertex (const X &, - const Point &, - const std::vector &); - - template - unsigned int - find_closest_vertex (const Mapping&, - const X &, - const Point &, - const std::vector &); - - template - std::vector::type> - find_cells_adjacent_to_vertex(const X &, const unsigned int); - - template - dealii::internal::ActiveCellIterator::type - find_active_cell_around_point (const X &, const Point &, const std::vector &); - - template - std::pair::type, Point > - find_active_cell_around_point (const Mapping &, - const X &, - const Point &, - const std::vector &); - - template - std::vector::type> - compute_active_cell_halo_layer (const X &, - const std::function::type&)> &); - - template - std::vector - compute_cell_halo_layer_on_level (const X &, - const std::function &, - const unsigned int); - - template - std::vector::type> - compute_ghost_cell_halo_layer (const X &); - - - template - std::vector::type> - compute_active_cell_layer_within_distance (const X &, - const std::function::type&)> &, - const double); - + find_closest_vertex_of_cell + (const typename Triangulation::active_cell_iterator &, + const Point &); template - std::vector::type> - compute_ghost_cell_layer_within_distance (const X &, const double); - + std::map + compute_local_to_global_vertex_index_map(const parallel::distributed::Triangulation &triangulation); template - std::pair< Point, Point > - compute_bounding_box (const X &, - const std::function::type&)> &); - + std::map > + extract_used_vertices(const Triangulation&mesh, + const Mapping &mapping); template - std::list > - get_finest_common_cells (const X &mesh_1, - const X &mesh_2); - + std::pair::active_cell_iterator, + Point > + find_active_cell_around_point(const Cache&, + const Point &, + const typename Triangulation::active_cell_iterator &, + const std::vector &); template - bool - have_same_coarse_mesh (const X &mesh_1, - const X &mesh_2); + std::tuple< std::vector< typename Triangulation< deal_II_dimension, deal_II_space_dimension>::active_cell_iterator >, + std::vector< std::vector< Point< deal_II_dimension > > >, std::vector< std::vector< unsigned int > > > + compute_point_locations(const Cache< deal_II_dimension, deal_II_space_dimension > &, + const std::vector< Point< deal_II_space_dimension > > &, + const typename Triangulation< deal_II_dimension, deal_II_space_dimension>::active_cell_iterator &); \} #endif } -// now also instantiate a few additional functions for parallel::distributed::Triangulation -for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS) -{ - -#if deal_II_dimension <= deal_II_space_dimension - namespace GridTools \{ - - template - unsigned int - find_closest_vertex_of_cell - (const typename Triangulation::active_cell_iterator &, - const Point &); - - template - std::map - compute_local_to_global_vertex_index_map(const parallel::distributed::Triangulation &triangulation); - - template - std::map > - extract_used_vertices(const Triangulation&mesh, - const Mapping &mapping); - - template - std::pair::active_cell_iterator, - Point > - find_active_cell_around_point(const Cache&, - const Point &, - const typename Triangulation::active_cell_iterator &, - const std::vector &); - - template - std::tuple< std::vector< typename Triangulation< deal_II_dimension, deal_II_space_dimension>::active_cell_iterator >, - std::vector< std::vector< Point< deal_II_dimension > > >, std::vector< std::vector< unsigned int > > > - compute_point_locations(const Cache< deal_II_dimension, deal_II_space_dimension > &, - const std::vector< Point< deal_II_space_dimension > > &, - const typename Triangulation< deal_II_dimension, deal_II_space_dimension>::active_cell_iterator &); - \} - -#endif -} - for (deal_II_space_dimension : SPACE_DIMENSIONS) @@ -284,14 +201,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS template void partition_multigrid_levels (Triangulation &); - template - std::pair::active_cell_iterator, - Point > - find_active_cell_around_point - (const hp::MappingCollection &, - const hp::DoFHandler &, - const Point &); - template void get_subdomain_association (const Triangulation &, std::vector &); @@ -351,141 +260,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS -} - -for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; Container : DOFHANDLER_TEMPLATES) -{ -#if deal_II_dimension <= deal_II_space_dimension - namespace GridTools \{ - - template - std::map< types::global_dof_index,std::vector::active_cell_iterator> > - get_dof_to_support_patch_map > - (Container &dof_handler); - - \} -#endif -} - - -for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; Container : TRIANGULATION_AND_DOFHANDLER_TEMPLATES) -{ -#if deal_II_dimension <= deal_II_space_dimension - namespace GridTools \{ - - template - std::vector::active_cell_iterator> - get_patch_around_cell > - (const Container::active_cell_iterator &cell); - - template - std::vector< Container::cell_iterator> - get_cells_at_coarsest_common_level > ( - const std::vector< Container::active_cell_iterator> & patch_cells); - - template - void build_triangulation_from_patch > ( - const std::vector::active_cell_iterator> &patch, - Triangulation::dimension,Container::space_dimension> &local_triangulation, - std::map::active_cell_iterator, - Container::active_cell_iterator > &patch_to_global_tria_map); - - \} -#endif -} - - -// instantiate the following functions only for the "sequential" containers. this -// is a misnomer here, however: the point is simply that we only instantiate -// these functions for certain *iterator* types, and the iterator types are -// the same for sequential and parallel containers; consequently, we get duplicate -// instantiation errors if we instantiate for *all* container types, rather than -// only the sequential ones -for (X : SEQUENTIAL_TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS) -{ -#if deal_II_dimension <= deal_II_space_dimension - namespace GridTools \{ - - template - bool orthogonal_equality (std::bitset<3> &, - const X::active_face_iterator&, - const X::active_face_iterator&, - const int, - const Tensor<1,deal_II_space_dimension> &, - const FullMatrix &); - - template - bool orthogonal_equality (std::bitset<3> &, - const X::face_iterator&, - const X::face_iterator&, - const int, - const Tensor<1,deal_II_space_dimension> &, - const FullMatrix &); - - template - bool orthogonal_equality (const X::active_face_iterator&, - const X::active_face_iterator&, - const int, - const Tensor<1,deal_II_space_dimension> &, - const FullMatrix &); - - template - bool orthogonal_equality (const X::face_iterator&, - const X::face_iterator&, - const int, - const Tensor<1,deal_II_space_dimension> &, - const FullMatrix &); - - template - void collect_periodic_faces (const X &, - const types::boundary_id, - const types::boundary_id, - const int, - std::vector > &, - const Tensor<1,X::space_dimension> &, - const FullMatrix &); - - template - void collect_periodic_faces (const X &, - const types::boundary_id, - const int, - std::vector > &, - const Tensor<1,X::space_dimension> &, - const FullMatrix &); - - \} -#endif -} - -for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS) -{ -#if deal_II_dimension <= deal_II_space_dimension -#if deal_II_dimension >= 2 - - namespace GridTools \{ - template - void - collect_periodic_faces > - (const parallel::distributed::Triangulation &, - const types::boundary_id, - const types::boundary_id, - const int, - std::vector::cell_iterator> > &, - const Tensor<1,parallel::distributed::Triangulation::space_dimension> &, - const FullMatrix &); - - template - void - collect_periodic_faces > - (const parallel::distributed::Triangulation &, - const types::boundary_id, - const int, - std::vector::cell_iterator> > &, - const Tensor<1,parallel::distributed::Triangulation::space_dimension> &, - const FullMatrix &); - \} -#endif -#endif } for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS) diff --git a/source/grid/grid_tools_dof_handlers.cc b/source/grid/grid_tools_dof_handlers.cc new file mode 100644 index 0000000000..2caf255138 --- /dev/null +++ b/source/grid/grid_tools_dof_handlers.cc @@ -0,0 +1,2050 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2001 - 2017 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. +// +// --------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + +namespace GridTools +{ + template class MeshType, int spacedim> + unsigned int + find_closest_vertex (const MeshType &mesh, + const Point &p, + const std::vector &marked_vertices) + { + // first get the underlying + // triangulation from the + // mesh and determine vertices + // and used vertices + const Triangulation &tria = mesh.get_triangulation(); + + const std::vector< Point > &vertices = tria.get_vertices(); + + Assert ( tria.get_vertices().size() == marked_vertices.size() || marked_vertices.size() ==0, + ExcDimensionMismatch(tria.get_vertices().size(), marked_vertices.size())); + + // If p is an element of marked_vertices, + // and q is that of used_Vertices, + // the vector marked_vertices does NOT + // contain unused vertices if p implies q. + // I.e., if p is true q must be true + // (if p is false, q could be false or true). + // p implies q logic is encapsulated in ~p|q. + Assert( marked_vertices.size()==0 + || + std::equal( marked_vertices.begin(), + marked_vertices.end(), + tria.get_used_vertices().begin(), + [](bool p, bool q) + { + return !p || q; + }), + ExcMessage("marked_vertices should be a subset of used vertices in the triangulation " + "but marked_vertices contains one or more vertices that are not used vertices!") ); + + // In addition, if a vector bools + // is specified (marked_vertices) + // marking all the vertices which + // could be the potentially closest + // vertex to the point, use it instead + // of used vertices + const std::vector &used = + (marked_vertices.size()==0) ? tria.get_used_vertices() : marked_vertices; + + // At the beginning, the first + // used vertex is the closest one + std::vector::const_iterator first = + std::find(used.begin(), used.end(), true); + + // Assert that at least one vertex + // is actually used + Assert(first != used.end(), ExcInternalError()); + + unsigned int best_vertex = std::distance(used.begin(), first); + double best_dist = (p - vertices[best_vertex]).norm_square(); + + // For all remaining vertices, test + // whether they are any closer + for (unsigned int j = best_vertex+1; j < vertices.size(); j++) + if (used[j]) + { + double dist = (p - vertices[j]).norm_square(); + if (dist < best_dist) + { + best_vertex = j; + best_dist = dist; + } + } + + return best_vertex; + } + + + + template class MeshType, int spacedim> + unsigned int + find_closest_vertex (const Mapping &mapping, + const MeshType &mesh, + const Point &p, + const std::vector &marked_vertices) + { + // Take a shortcut in the simple case. + if (mapping.preserves_vertex_locations() == true) + return find_closest_vertex(mesh, p, marked_vertices); + + // first get the underlying + // triangulation from the + // mesh and determine vertices + // and used vertices + const Triangulation &tria = mesh.get_triangulation(); + + auto vertices = extract_used_vertices(tria, mapping); + + Assert ( tria.get_vertices().size() == marked_vertices.size() || marked_vertices.size() ==0, + ExcDimensionMismatch(tria.get_vertices().size(), marked_vertices.size())); + + // If p is an element of marked_vertices, + // and q is that of used_Vertices, + // the vector marked_vertices does NOT + // contain unused vertices if p implies q. + // I.e., if p is true q must be true + // (if p is false, q could be false or true). + // p implies q logic is encapsulated in ~p|q. + Assert( marked_vertices.size()==0 + || + std::equal( marked_vertices.begin(), + marked_vertices.end(), + tria.get_used_vertices().begin(), + [](bool p, bool q) + { + return !p || q; + }), + ExcMessage("marked_vertices should be a subset of used vertices in the triangulation " + "but marked_vertices contains one or more vertices that are not used vertices!") ); + + // Remove from the map unwanted elements. + if (marked_vertices.size()) + for (auto it = vertices.begin(); it != vertices.end(); ) + { + if (marked_vertices[it->first] == false) + { + vertices.erase(it++); + } + else + { + ++it; + } + } + + return find_closest_vertex(vertices, p); + } + + + + template class MeshType, int spacedim> +#ifndef _MSC_VER + std::vector::active_cell_iterator> +#else + std::vector >::type> +#endif + find_cells_adjacent_to_vertex(const MeshType &mesh, + const unsigned int vertex) + { + // make sure that the given vertex is + // an active vertex of the underlying + // triangulation + Assert(vertex < mesh.get_triangulation().n_vertices(), + ExcIndexRange(0,mesh.get_triangulation().n_vertices(),vertex)); + Assert(mesh.get_triangulation().get_used_vertices()[vertex], + ExcVertexNotUsed(vertex)); + + // use a set instead of a vector + // to ensure that cells are inserted only + // once + std::set >::type> adjacent_cells; + + typename dealii::internal::ActiveCellIterator >::type + cell = mesh.begin_active(), + endc = mesh.end(); + + // go through all active cells and look if the vertex is part of that cell + // + // in 1d, this is all we need to care about. in 2d/3d we also need to worry + // that the vertex might be a hanging node on a face or edge of a cell; in + // this case, we would want to add those cells as well on whose faces the + // vertex is located but for which it is not a vertex itself. + // + // getting this right is a lot simpler in 2d than in 3d. in 2d, a hanging + // node can only be in the middle of a face and we can query the neighboring + // cell from the current cell. on the other hand, in 3d a hanging node + // vertex can also be on an edge but there can be many other cells on + // this edge and we can not access them from the cell we are currently + // on. + // + // so, in the 3d case, if we run the algorithm as in 2d, we catch all + // those cells for which the vertex we seek is on a *subface*, but we + // miss the case of cells for which the vertex we seek is on a + // sub-edge for which there is no corresponding sub-face (because the + // immediate neighbor behind this face is not refined), see for example + // the bits/find_cells_adjacent_to_vertex_6 testcase. thus, if we + // haven't yet found the vertex for the current cell we also need to + // look at the mid-points of edges + // + // as a final note, deciding whether a neighbor is actually coarser is + // simple in the case of isotropic refinement (we just need to look at + // the level of the current and the neighboring cell). however, this + // isn't so simple if we have used anisotropic refinement since then + // the level of a cell is not indicative of whether it is coarser or + // not than the current cell. ultimately, we want to add all cells on + // which the vertex is, independent of whether they are coarser or + // finer and so in the 2d case below we simply add *any* *active* neighbor. + // in the worst case, we add cells multiple times to the adjacent_cells + // list, but std::set throws out those cells already entered + for (; cell != endc; ++cell) + { + for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; v++) + if (cell->vertex_index(v) == vertex) + { + // OK, we found a cell that contains + // the given vertex. We add it + // to the list. + adjacent_cells.insert(cell); + + // as explained above, in 2+d we need to check whether + // this vertex is on a face behind which there is a + // (possibly) coarser neighbor. if this is the case, + // then we need to also add this neighbor + if (dim >= 2) + for (unsigned int vface = 0; vface < dim; vface++) + { + const unsigned int face = + GeometryInfo::vertex_to_face[v][vface]; + + if (!cell->at_boundary(face) + && + cell->neighbor(face)->active()) + { + // there is a (possibly) coarser cell behind a + // face to which the vertex belongs. the + // vertex we are looking at is then either a + // vertex of that coarser neighbor, or it is a + // hanging node on one of the faces of that + // cell. in either case, it is adjacent to the + // vertex, so add it to the list as well (if + // the cell was already in the list then the + // std::set makes sure that we get it only + // once) + adjacent_cells.insert (cell->neighbor(face)); + } + } + + // in any case, we have found a cell, so go to the next cell + goto next_cell; + } + + // in 3d also loop over the edges + if (dim >= 3) + { + for (unsigned int e=0; e::lines_per_cell; ++e) + if (cell->line(e)->has_children()) + // the only place where this vertex could have been + // hiding is on the mid-edge point of the edge we + // are looking at + if (cell->line(e)->child(0)->vertex_index(1) == vertex) + { + adjacent_cells.insert(cell); + + // jump out of this tangle of nested loops + goto next_cell; + } + } + + // in more than 3d we would probably have to do the same as + // above also for even lower-dimensional objects + Assert (dim <= 3, ExcNotImplemented()); + + // move on to the next cell if we have found the + // vertex on the current one +next_cell: + ; + } + + // if this was an active vertex then there needs to have been + // at least one cell to which it is adjacent! + Assert (adjacent_cells.size() > 0, ExcInternalError()); + + // return the result as a vector, rather than the set we built above + return + std::vector >::type> + (adjacent_cells.begin(), adjacent_cells.end()); + } + + + + namespace + { + template class MeshType, int spacedim> + void find_active_cell_around_point_internal + (const MeshType &mesh, +#ifndef _MSC_VER + std::set::active_cell_iterator> &searched_cells, + std::set::active_cell_iterator> &adjacent_cells) +#else + std::set >::type> &searched_cells, + std::set >::type> &adjacent_cells) +#endif + { +#ifndef _MSC_VER + typedef typename MeshType::active_cell_iterator cell_iterator; +#else + typedef typename dealii::internal::ActiveCellIterator >::type cell_iterator; +#endif + + // update the searched cells + searched_cells.insert(adjacent_cells.begin(), adjacent_cells.end()); + // now we to collect all neighbors + // of the cells in adjacent_cells we + // have not yet searched. + std::set adjacent_cells_new; + + typename std::set::const_iterator + cell = adjacent_cells.begin(), + endc = adjacent_cells.end(); + for (; cell != endc; ++cell) + { + std::vector active_neighbors; + get_active_neighbors >(*cell, active_neighbors); + for (unsigned int i=0; i class MeshType, int spacedim> +#ifndef _MSC_VER + typename MeshType::active_cell_iterator +#else + typename dealii::internal::ActiveCellIterator >::type +#endif + find_active_cell_around_point (const MeshType &mesh, + const Point &p, + const std::vector &marked_vertices) + { + return + find_active_cell_around_point + (StaticMappingQ1::mapping, + mesh, p, marked_vertices).first; + } + + + template class MeshType, int spacedim> +#ifndef _MSC_VER + std::pair::active_cell_iterator, Point > +#else + std::pair >::type, Point > +#endif + find_active_cell_around_point (const Mapping &mapping, + const MeshType &mesh, + const Point &p, + const std::vector &marked_vertices) + { + typedef typename dealii::internal::ActiveCellIterator >::type active_cell_iterator; + + // The best distance is set to the + // maximum allowable distance from + // the unit cell; we assume a + // max. deviation of 1e-10 + double best_distance = 1e-10; + int best_level = -1; + std::pair > best_cell; + + // Find closest vertex and determine + // all adjacent cells + std::vector adjacent_cells_tmp + = find_cells_adjacent_to_vertex(mesh, + find_closest_vertex(mapping, mesh, p, marked_vertices)); + + // Make sure that we have found + // at least one cell adjacent to vertex. + Assert(adjacent_cells_tmp.size()>0, ExcInternalError()); + + // Copy all the cells into a std::set + std::set adjacent_cells (adjacent_cells_tmp.begin(), + adjacent_cells_tmp.end()); + std::set searched_cells; + + // Determine the maximal number of cells + // in the grid. + // As long as we have not found + // the cell and have not searched + // every cell in the triangulation, + // we keep on looking. + const unsigned int n_active_cells = mesh.get_triangulation().n_active_cells(); + bool found = false; + unsigned int cells_searched = 0; + while (!found && cells_searched < n_active_cells) + { + typename std::set::const_iterator + cell = adjacent_cells.begin(), + endc = adjacent_cells.end(); + for (; cell != endc; ++cell) + { + try + { + const Point p_cell = mapping.transform_real_to_unit_cell(*cell, p); + + // calculate the infinity norm of + // the distance vector to the unit cell. + const double dist = GeometryInfo::distance_to_unit_cell(p_cell); + + // We compare if the point is inside the + // unit cell (or at least not too far + // outside). If it is, it is also checked + // that the cell has a more refined state + if ((dist < best_distance) + || + ((dist == best_distance) + && + ((*cell)->level() > best_level))) + { + found = true; + best_distance = dist; + best_level = (*cell)->level(); + best_cell = std::make_pair(*cell, p_cell); + } + } + catch (typename MappingQGeneric::ExcTransformationFailed &) + { + // ok, the transformation + // failed presumably + // because the point we + // are looking for lies + // outside the current + // cell. this means that + // the current cell can't + // be the cell around the + // point, so just ignore + // this cell and move on + // to the next + } + } + + // update the number of cells searched + cells_searched += adjacent_cells.size(); + + // if the user provided a custom mask for vertices, + // terminate the search without trying to expand the search + // to all cells of the triangulation, as done below. + if (marked_vertices.size() > 0) + cells_searched = n_active_cells; + + // if we have not found the cell in + // question and have not yet searched every + // cell, we expand our search to + // all the not already searched neighbors of + // the cells in adjacent_cells. This is + // what find_active_cell_around_point_internal + // is for. + if (!found && cells_searched < n_active_cells) + { + find_active_cell_around_point_internal + (mesh, searched_cells, adjacent_cells); + } + } + + AssertThrow (best_cell.first.state() == IteratorState::valid, + ExcPointNotFound(p)); + + return best_cell; + } + + + + template + std::vector + compute_active_cell_halo_layer + (const MeshType &mesh, + const std::function &predicate) + { + std::vector active_halo_layer; + std::vector locally_active_vertices_on_subdomain (mesh.get_triangulation().n_vertices(), + false); + + // Find the cells for which the predicate is true + // These are the cells around which we wish to construct + // the halo layer + for (typename MeshType::active_cell_iterator + cell = mesh.begin_active(); + cell != mesh.end(); ++cell) + if (predicate(cell)) // True predicate --> Part of subdomain + for (unsigned int v=0; v::vertices_per_cell; ++v) + locally_active_vertices_on_subdomain[cell->vertex_index(v)] = true; + + // Find the cells that do not conform to the predicate + // but share a vertex with the selected subdomain + // These comprise the halo layer + for (typename MeshType::active_cell_iterator + cell = mesh.begin_active(); + cell != mesh.end(); ++cell) + if (!predicate(cell)) // False predicate --> Potential halo cell + for (unsigned int v=0; v::vertices_per_cell; ++v) + if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] == true) + { + active_halo_layer.push_back(cell); + break; + } + + return active_halo_layer; + } + + + + template + std::vector + compute_cell_halo_layer_on_level + (const MeshType &mesh, + const std::function &predicate, + const unsigned int level) + { + std::vector level_halo_layer; + std::vector locally_active_vertices_on_level_subdomain (mesh.get_triangulation().n_vertices(), + false); + + // Find the cells for which the predicate is true + // These are the cells around which we wish to construct + // the halo layer + for (typename MeshType::cell_iterator + cell = mesh.begin(level); + cell != mesh.end(level); ++cell) + if (predicate(cell)) // True predicate --> Part of subdomain + for (unsigned int v=0; v::vertices_per_cell; ++v) + locally_active_vertices_on_level_subdomain[cell->vertex_index(v)] = true; + + // Find the cells that do not conform to the predicate + // but share a vertex with the selected subdomain on that level + // These comprise the halo layer + for (typename MeshType::cell_iterator + cell = mesh.begin(level); + cell != mesh.end(level); ++cell) + if (!predicate(cell)) // False predicate --> Potential halo cell + for (unsigned int v=0; v::vertices_per_cell; ++v) + if (locally_active_vertices_on_level_subdomain[cell->vertex_index(v)] == true) + { + level_halo_layer.push_back(cell); + break; + } + + return level_halo_layer; + } + + + namespace + { + template + bool + contains_locally_owned_cells (const std::vector &cells) + { + for (typename std::vector::const_iterator + it = cells.begin(); it != cells.end(); ++it) + { + if ((*it)->is_locally_owned()) + return true; + } + return false; + } + + template + bool + contains_artificial_cells (const std::vector &cells) + { + for (typename std::vector::const_iterator + it = cells.begin(); it != cells.end(); ++it) + { + if ((*it)->is_artificial()) + return true; + } + return false; + } + } + + + + + template + std::vector + compute_ghost_cell_halo_layer (const MeshType &mesh) + { + std::function predicate + = IteratorFilters::LocallyOwnedCell(); + + const std::vector + active_halo_layer = compute_active_cell_halo_layer (mesh, predicate); + + // Check that we never return locally owned or artificial cells + // What is left should only be the ghost cells + Assert(contains_locally_owned_cells(active_halo_layer) == false, + ExcMessage("Halo layer contains locally owned cells")); + Assert(contains_artificial_cells(active_halo_layer) == false, + ExcMessage("Halo layer contains artificial cells")); + + return active_halo_layer; + } + + + + template + std::vector + compute_active_cell_layer_within_distance + (const MeshType &mesh, + const std::function &predicate, + const double layer_thickness) + { + std::vector subdomain_boundary_cells, active_cell_layer_within_distance; + std::vector vertices_outside_subdomain ( mesh.get_triangulation().n_vertices(), + false); + + const unsigned int spacedim = MeshType::space_dimension; + + unsigned int n_non_predicate_cells = 0; // Number of non predicate cells + + // Find the layer of cells for which predicate is true and that + // are on the boundary with other cells. These are + // subdomain boundary cells. + + // Find the cells for which the predicate is false + // These are the cells which are around the predicate subdomain + for ( typename MeshType::active_cell_iterator + cell = mesh.begin_active(); + cell != mesh.end(); ++cell) + if ( !predicate(cell)) // Negation of predicate --> Not Part of subdomain + { + for (unsigned int v=0; v::vertices_per_cell; ++v) + vertices_outside_subdomain[cell->vertex_index(v)] = true; + n_non_predicate_cells++; + } + + // If all the active cells conform to the predicate + // or if none of the active cells conform to the predicate + // there is no active cell layer around the predicate + // subdomain (within any distance) + if ( n_non_predicate_cells == 0 || n_non_predicate_cells == mesh.get_triangulation().n_active_cells() ) + return std::vector(); + + // Find the cells that conform to the predicate + // but share a vertex with the cell not in the predicate subdomain + for ( typename MeshType::active_cell_iterator + cell = mesh.begin_active(); + cell != mesh.end(); ++cell) + if ( predicate(cell)) // True predicate --> Potential boundary cell of the subdomain + for (unsigned int v=0; v::vertices_per_cell; ++v) + if (vertices_outside_subdomain[cell->vertex_index(v)] == true) + { + subdomain_boundary_cells.push_back(cell); + break; // No need to go through remaining vertices + } + + // To cheaply filter out some cells located far away from the predicate subdomain, + // get the bounding box of the predicate subdomain. + std::pair< Point, Point > bounding_box = compute_bounding_box( mesh, + predicate ); + + // DOUBLE_EPSILON to compare really close double values + const double &DOUBLE_EPSILON = 100.*std::numeric_limits::epsilon(); + + // Add layer_thickness to the bounding box + for ( unsigned int d=0; d > subdomain_boundary_cells_centers; // cache all the subdomain boundary cells centers here + std::vector subdomain_boundary_cells_radii; // cache all the subdomain boundary cells radii + subdomain_boundary_cells_centers.reserve (subdomain_boundary_cells.size()); + subdomain_boundary_cells_radii.reserve (subdomain_boundary_cells.size()); + // compute cell radius for each boundary cell of the predicate subdomain + for ( typename std::vector::const_iterator + subdomain_boundary_cell_iterator = subdomain_boundary_cells.begin(); + subdomain_boundary_cell_iterator != subdomain_boundary_cells.end(); ++subdomain_boundary_cell_iterator ) + { + const std::pair, double> & + subdomain_boundary_cell_enclosing_ball = (*subdomain_boundary_cell_iterator)->enclosing_ball(); + + subdomain_boundary_cells_centers.push_back( subdomain_boundary_cell_enclosing_ball.first); + subdomain_boundary_cells_radii.push_back( subdomain_boundary_cell_enclosing_ball.second); + } + AssertThrow( subdomain_boundary_cells_radii.size() == subdomain_boundary_cells_centers.size(), + ExcInternalError()); + + // Find the cells that are within layer_thickness of predicate subdomain boundary + // distance but are inside the extended bounding box. + // Most cells might be outside the extended bounding box, so we could skip them. + // Those cells that are inside the extended bounding box but are not part of the + // predicate subdomain are possible candidates to be within the distance to the + // boundary cells of the predicate subdomain. + for ( typename MeshType::active_cell_iterator + cell = mesh.begin_active(); + cell != mesh.end(); ++cell) + { + // Ignore all the cells that are in the predicate subdomain + if ( predicate(cell)) + continue; + + const std::pair, double> &cell_enclosing_ball + = cell->enclosing_ball(); + + const Point &cell_enclosing_ball_center = cell_enclosing_ball.first; + const double &cell_enclosing_ball_radius = cell_enclosing_ball.second; + + bool cell_inside = true; // reset for each cell + + for (unsigned int d = 0; d < spacedim; ++d) + cell_inside &= (cell_enclosing_ball_center[d] + cell_enclosing_ball_radius > bounding_box.first[d]) + && (cell_enclosing_ball_center[d] - cell_enclosing_ball_radius < bounding_box.second[d]); + // cell_inside is true if its enclosing ball intersects the extended bounding box + + // Ignore all the cells that are outside the extended bounding box + if (cell_inside) + for (unsigned int i =0; i< subdomain_boundary_cells_radii.size(); ++i) + if ( cell_enclosing_ball_center.distance_square(subdomain_boundary_cells_centers[i]) + < Utilities::fixed_power<2>( cell_enclosing_ball_radius + + subdomain_boundary_cells_radii[i] + + layer_thickness + DOUBLE_EPSILON )) + { + active_cell_layer_within_distance.push_back(cell); + break; // Exit the loop checking all the remaining subdomain boundary cells + } + + } + return active_cell_layer_within_distance; + } + + + + template + std::vector + compute_ghost_cell_layer_within_distance ( const MeshType &mesh, const double layer_thickness) + { + IteratorFilters::LocallyOwnedCell locally_owned_cell_predicate; + std::function predicate (locally_owned_cell_predicate); + + const std::vector + ghost_cell_layer_within_distance = compute_active_cell_layer_within_distance (mesh, predicate, layer_thickness); + + // Check that we never return locally owned or artificial cells + // What is left should only be the ghost cells + Assert(contains_locally_owned_cells(ghost_cell_layer_within_distance) == false, + ExcMessage("Ghost cells within layer_thickness contains locally owned cells.")); + Assert(contains_artificial_cells(ghost_cell_layer_within_distance) == false, + ExcMessage("Ghost cells within layer_thickness contains artificial cells." + "The function compute_ghost_cell_layer_within_distance " + "is probably called while using parallel::distributed::Triangulation. " + "In such case please refer to the description of this function.")); + + return ghost_cell_layer_within_distance; + } + + + + template < class MeshType> + std::pair< Point, Point > + compute_bounding_box + ( const MeshType &mesh, + const std::function &predicate ) + { + std::vector locally_active_vertices_on_subdomain (mesh.get_triangulation().n_vertices(), + false); + + const unsigned int spacedim = MeshType::space_dimension; + + // Two extreme points can define the bounding box + // around the active cells that conform to the given predicate. + Point maxp, minp; + + // initialize minp and maxp with the first predicate cell center + for ( typename MeshType::active_cell_iterator + cell = mesh.begin_active(); + cell != mesh.end(); ++cell) + if ( predicate(cell)) + { + minp = cell->center(); + maxp = cell->center(); + break; + } + + // Run through all the cells to check if it belongs to predicate domain, + // if it belongs to the predicate domain, extend the bounding box. + for ( typename MeshType::active_cell_iterator + cell = mesh.begin_active(); + cell != mesh.end(); ++cell) + if (predicate(cell)) // True predicate --> Part of subdomain + for (unsigned int v=0; v::vertices_per_cell; ++v) + if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] == false) + { + locally_active_vertices_on_subdomain[cell->vertex_index(v)] = true; + for ( unsigned int d=0; dvertex(v)[d]); + maxp[d] = std::max( maxp[d], cell->vertex(v)[d]); + } + } + + return std::make_pair( minp, maxp ); + } + + + + template + std::list > + get_finest_common_cells (const MeshType &mesh_1, + const MeshType &mesh_2) + { + Assert (have_same_coarse_mesh (mesh_1, mesh_2), + ExcMessage ("The two meshes must be represent triangulations that " + "have the same coarse meshes")); + + // the algorithm goes as follows: + // first, we fill a list with pairs + // of iterators common to the two + // meshes on the coarsest + // level. then we traverse the + // list; each time, we find a pair + // of iterators for which both + // correspond to non-active cells, + // we delete this item and push the + // pairs of iterators to their + // children to the back. if these + // again both correspond to + // non-active cells, we will get to + // the later on for further + // consideration + typedef + std::list > + CellList; + + CellList cell_list; + + // first push the coarse level cells + typename MeshType::cell_iterator + cell_1 = mesh_1.begin(0), + cell_2 = mesh_2.begin(0); + for (; cell_1 != mesh_1.end(0); ++cell_1, ++cell_2) + cell_list.emplace_back (cell_1, cell_2); + + // then traverse list as described + // above + typename CellList::iterator cell_pair = cell_list.begin(); + while (cell_pair != cell_list.end()) + { + // if both cells in this pair + // have children, then erase + // this element and push their + // children instead + if (cell_pair->first->has_children() + && + cell_pair->second->has_children()) + { + Assert(cell_pair->first->refinement_case()== + cell_pair->second->refinement_case(), ExcNotImplemented()); + for (unsigned int c=0; cfirst->n_children(); ++c) + cell_list.emplace_back (cell_pair->first->child(c), + cell_pair->second->child(c)); + + // erasing an iterator + // keeps other iterators + // valid, so already + // advance the present + // iterator by one and then + // delete the element we've + // visited before + const typename CellList::iterator previous_cell_pair = cell_pair; + ++cell_pair; + + cell_list.erase (previous_cell_pair); + } + else + // both cells are active, do + // nothing + ++cell_pair; + } + + // just to make sure everything is ok, + // validate that all pairs have at least one + // active iterator or have different + // refinement_cases + for (cell_pair = cell_list.begin(); cell_pair != cell_list.end(); ++cell_pair) + Assert (cell_pair->first->active() + || + cell_pair->second->active() + || + (cell_pair->first->refinement_case() + != cell_pair->second->refinement_case()), + ExcInternalError()); + + return cell_list; + } + + + + template + bool + have_same_coarse_mesh (const Triangulation &mesh_1, + const Triangulation &mesh_2) + { + // make sure the two meshes have + // the same number of coarse cells + if (mesh_1.n_cells (0) != mesh_2.n_cells (0)) + return false; + + // if so, also make sure they have + // the same vertices on the cells + // of the coarse mesh + typename Triangulation::cell_iterator + cell_1 = mesh_1.begin(0), + cell_2 = mesh_2.begin(0), + endc = mesh_1.end(0); + for (; cell_1!=endc; ++cell_1, ++cell_2) + for (unsigned int v=0; v::vertices_per_cell; ++v) + if (cell_1->vertex(v) != cell_2->vertex(v)) + return false; + + // if we've gotten through all + // this, then the meshes really + // seem to have a common coarse + // mesh + return true; + } + + + + template + bool + have_same_coarse_mesh (const MeshType &mesh_1, + const MeshType &mesh_2) + { + return have_same_coarse_mesh (mesh_1.get_triangulation(), + mesh_2.get_triangulation()); + } + + + + template + std::pair::active_cell_iterator, Point > + find_active_cell_around_point (const hp::MappingCollection &mapping, + const hp::DoFHandler &mesh, + const Point &p) + { + Assert ((mapping.size() == 1) || + (mapping.size() == mesh.get_fe_collection().size()), + ExcMessage ("Mapping collection needs to have either size 1 " + "or size equal to the number of elements in " + "the FECollection.")); + + typedef typename hp::DoFHandler::active_cell_iterator cell_iterator; + + std::pair > best_cell; + //If we have only one element in the MappingCollection, + //we use find_active_cell_around_point using only one + //mapping. + if (mapping.size() == 1) + best_cell = find_active_cell_around_point(mapping[0], mesh, p); + else + { + + + // The best distance is set to the + // maximum allowable distance from + // the unit cell; we assume a + // max. deviation of 1e-10 + double best_distance = 1e-10; + int best_level = -1; + + + // Find closest vertex and determine + // all adjacent cells + unsigned int vertex = find_closest_vertex(mesh, p); + + std::vector adjacent_cells_tmp = + find_cells_adjacent_to_vertex(mesh, vertex); + + // Make sure that we have found + // at least one cell adjacent to vertex. + Assert(adjacent_cells_tmp.size()>0, ExcInternalError()); + + // Copy all the cells into a std::set + std::set adjacent_cells(adjacent_cells_tmp.begin(), adjacent_cells_tmp.end()); + std::set searched_cells; + + // Determine the maximal number of cells + // in the grid. + // As long as we have not found + // the cell and have not searched + // every cell in the triangulation, + // we keep on looking. + const unsigned int n_cells = mesh.get_triangulation().n_cells(); + bool found = false; + unsigned int cells_searched = 0; + while (!found && cells_searched < n_cells) + { + typename std::set::const_iterator + cell = adjacent_cells.begin(), + endc = adjacent_cells.end(); + for (; cell != endc; ++cell) + { + try + { + const Point p_cell = mapping[(*cell)->active_fe_index()].transform_real_to_unit_cell(*cell, p); + + + // calculate the infinity norm of + // the distance vector to the unit cell. + const double dist = GeometryInfo::distance_to_unit_cell(p_cell); + + // We compare if the point is inside the + // unit cell (or at least not too far + // outside). If it is, it is also checked + // that the cell has a more refined state + if (dist < best_distance || + (dist == best_distance && (*cell)->level() > best_level)) + { + found = true; + best_distance = dist; + best_level = (*cell)->level(); + best_cell = std::make_pair(*cell, p_cell); + } + } + catch (typename MappingQGeneric::ExcTransformationFailed &) + { + // ok, the transformation + // failed presumably + // because the point we + // are looking for lies + // outside the current + // cell. this means that + // the current cell can't + // be the cell around the + // point, so just ignore + // this cell and move on + // to the next + } + } + //udpate the number of cells searched + cells_searched += adjacent_cells.size(); + // if we have not found the cell in + // question and have not yet searched every + // cell, we expand our search to + // all the not already searched neighbors of + // the cells in adjacent_cells. + if (!found && cells_searched < n_cells) + { + find_active_cell_around_point_internal + (mesh, searched_cells, adjacent_cells); + } + + } + } + + AssertThrow (best_cell.first.state() == IteratorState::valid, + ExcPointNotFound(p)); + + return best_cell; + } + + + template + std::vector + get_patch_around_cell(const typename MeshType::active_cell_iterator &cell) + { + Assert (cell->is_locally_owned(), + ExcMessage ("This function only makes sense if the cell for " + "which you are asking for a patch, is locally " + "owned.")); + + std::vector patch; + patch.push_back (cell); + for (unsigned int face_number=0; face_number::faces_per_cell; ++face_number) + if (cell->face(face_number)->at_boundary()==false) + { + if (cell->neighbor(face_number)->has_children() == false) + patch.push_back (cell->neighbor(face_number)); + else + // the neighbor is refined. in 2d/3d, we can simply ask for the children + // of the neighbor because they can not be further refined and, + // consequently, the children is active + if (MeshType::dimension > 1) + { + for (unsigned int subface=0; subfaceface(face_number)->n_children(); ++subface) + patch.push_back (cell->neighbor_child_on_subface (face_number, subface)); + } + else + { + // in 1d, we need to work a bit harder: iterate until we find + // the child by going from cell to child to child etc + typename MeshType::cell_iterator neighbor + = cell->neighbor (face_number); + while (neighbor->has_children()) + neighbor = neighbor->child(1-face_number); + + Assert (neighbor->neighbor(1-face_number) == cell, ExcInternalError()); + patch.push_back (neighbor); + } + } + return patch; + } + + + + template + std::vector + get_cells_at_coarsest_common_level ( + const std::vector &patch) + { + Assert (patch.size() > 0, ExcMessage("Vector containing patch cells should not be an empty vector!")); + // In order to extract the set of cells with the coarsest common level from the give vector of cells: + // First it finds the number associated with the minimum level of refinmenet, namely "min_level" + int min_level = patch[0]->level(); + + for (unsigned int i=0; ilevel() ); + std::set uniform_cells; + typename std::vector::const_iterator patch_cell; + // it loops through all cells of the input vector + for (patch_cell=patch.begin(); patch_cell!=patch.end () ; ++patch_cell) + { + // If the refinement level of each cell i the loop be equal to the min_level, so that + // that cell inserted into the set of uniform_cells, as the set of cells with the coarsest common refinement level + if ((*patch_cell)->level() == min_level) + uniform_cells.insert (*patch_cell); + else + // If not, it asks for the parent of the cell, until it finds the parent cell + // with the refinement level equal to the min_level and inserts that parent cell into the + // the set of uniform_cells, as the set of cells with the coarsest common refinement level. + { + typename Container::cell_iterator parent = *patch_cell; + + while (parent->level() > min_level) + parent = parent-> parent(); + uniform_cells.insert (parent); + } + } + + return std::vector (uniform_cells.begin(), + uniform_cells.end()); + } + + + + template + void build_triangulation_from_patch(const std::vector &patch, + Triangulation &local_triangulation, + std::map::active_cell_iterator, + typename Container::active_cell_iterator> &patch_to_global_tria_map) + + { + const std::vector uniform_cells = + get_cells_at_coarsest_common_level (patch); + // First it creates triangulation from the vector of "uniform_cells" + local_triangulation.clear(); + std::vector > vertices; + const unsigned int n_uniform_cells=uniform_cells.size(); + std::vector > cells(n_uniform_cells); + unsigned int k=0;// for enumerating cells + unsigned int i=0;// for enumerating vertices + typename std::vector::const_iterator uniform_cell; + for (uniform_cell=uniform_cells.begin(); uniform_cell!=uniform_cells.end(); ++uniform_cell) + { + for (unsigned int v=0; v::vertices_per_cell; ++v) + { + Point position=(*uniform_cell)->vertex (v); + bool repeat_vertex=false; + + for (unsigned int m=0; m::cell_iterator, + typename Container::cell_iterator> patch_to_global_tria_map_tmp; + for (typename Triangulation::cell_iterator coarse_cell = local_triangulation.begin(); + coarse_cell != local_triangulation.end(); ++coarse_cell, ++index) + { + patch_to_global_tria_map_tmp.insert (std::make_pair(coarse_cell, uniform_cells[index])); + // To ensure that the cells with the same coordinates (here, we compare their centers) are mapped into each other. + + Assert(coarse_cell->center().distance( uniform_cells[index]->center())<=1e-15*coarse_cell->diameter(), + ExcInternalError()); + } + bool refinement_necessary; + // In this loop we start to do refinement on the above coarse triangulation to reach + // to the same level of refinement as the patch cells are really on + do + { + refinement_necessary = false; + for (typename Triangulation::active_cell_iterator + active_tria_cell = local_triangulation.begin_active(); + active_tria_cell != local_triangulation.end(); ++active_tria_cell) + { + if (patch_to_global_tria_map_tmp[active_tria_cell]->has_children()) + { + active_tria_cell -> set_refine_flag(); + refinement_necessary = true; + } + else for (unsigned int i=0; i::vertices_per_cell; ++v) + active_tria_cell->vertex(v) = patch[i]->vertex(v); + + Assert(active_tria_cell->center().distance(patch_to_global_tria_map_tmp[active_tria_cell]->center()) + <=1e-15*active_tria_cell->diameter(), ExcInternalError()); + + active_tria_cell->set_user_flag(); + break; + } + } + } + + if (refinement_necessary) + { + local_triangulation.execute_coarsening_and_refinement (); + + for (typename Triangulation::cell_iterator + cell = local_triangulation.begin(); + cell != local_triangulation.end(); ++cell) + { + + if (patch_to_global_tria_map_tmp.find(cell)!=patch_to_global_tria_map_tmp.end()) + { + if (cell-> has_children()) + { + // Note: Since the cell got children, then it should not be in the map anymore + // children may be added into the map, instead + + // these children may not yet be in the map + for (unsigned int c=0; cn_children(); ++c) + { + if (patch_to_global_tria_map_tmp.find(cell->child(c)) == + patch_to_global_tria_map_tmp.end()) + { + patch_to_global_tria_map_tmp.insert (std::make_pair(cell->child(c), + patch_to_global_tria_map_tmp[cell]->child(c))); + + // One might be tempted to assert that the cell + // being added here has the same center as the + // equivalent cell in the global triangulation, + // but it may not be the case. For triangulations + // that have been perturbed or smoothed, the cell + // indices and levels may be the same, but the + // vertex locations may not. We adjust + // the vertices of the cells that have no + // children (ie the active cells) to be + // consistent with the global triangulation + // later on and add assertions at that time + // to guarantee the cells in the + // local_triangulation are physically at the same + // locations of the cells in the patch of the + // global triangulation. + + } + } + // The parent cell whose children were added + // into the map should be deleted from the map + patch_to_global_tria_map_tmp.erase(cell); + } + } + } + } + + } + while (refinement_necessary); + + + // Last assertion check to make sure we have the right cells and centers + // in the map, and hence the correct vertices of the triangulation + for (typename Triangulation::cell_iterator + cell = local_triangulation.begin(); + cell != local_triangulation.end(); ++cell) + { + if (cell->user_flag_set() ) + { + Assert(patch_to_global_tria_map_tmp.find(cell) != patch_to_global_tria_map_tmp.end(), + ExcInternalError() ); + + Assert(cell->center().distance( patch_to_global_tria_map_tmp[cell]->center())<=1e-15*cell->diameter(), + ExcInternalError()); + } + } + + + typename std::map::cell_iterator, + typename Container::cell_iterator>::iterator map_tmp_it = + patch_to_global_tria_map_tmp.begin(),map_tmp_end = patch_to_global_tria_map_tmp.end(); + // Now we just need to take the temporary map of pairs of type cell_iterator "patch_to_global_tria_map_tmp" + // making pair of active_cell_iterators so that filling out the final map "patch_to_global_tria_map" + for (; map_tmp_it!=map_tmp_end; ++map_tmp_it) + patch_to_global_tria_map[map_tmp_it->first] = map_tmp_it->second; + } + + + + + template + std::map< types::global_dof_index,std::vector > + get_dof_to_support_patch_map(DoFHandlerType &dof_handler) + { + + // This is the map from global_dof_index to + // a set of cells on patch. We first map into + // a set because it is very likely that we + // will attempt to add a cell more than once + // to a particular patch and we want to preserve + // uniqueness of cell iterators. std::set does this + // automatically for us. Later after it is all + // constructed, we will copy to a map of vectors + // since that is the prefered output for other + // functions. + std::map< types::global_dof_index,std::set > dof_to_set_of_cells_map; + + std::vector local_dof_indices; + std::vector local_face_dof_indices; + std::vector local_line_dof_indices; + + // a place to save the dof_handler user flags and restore them later + // to maintain const of dof_handler. + std::vector user_flags; + + + // in 3d, we need pointers from active lines to the + // active parent lines, so we construct it as needed. + std::map lines_to_parent_lines_map; + if (DoFHandlerType::dimension == 3) + { + + // save user flags as they will be modified and then later restored + dof_handler.get_triangulation().save_user_flags(user_flags); + const_cast &>(dof_handler.get_triangulation()).clear_user_flags (); + + + typename DoFHandlerType::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + // We only want lines that are locally_relevant + // although it doesn't hurt to have lines that + // are children of ghost cells since there are + // few and we don't have to use them. + if (cell->is_artificial() == false) + { + for (unsigned int l=0; l::lines_per_cell; ++l) + if (cell->line(l)->has_children()) + for (unsigned int c=0; cline(l)->n_children(); ++c) + { + lines_to_parent_lines_map[cell->line(l)->child(c)] = cell->line(l); + // set flags to know that child + // line has an active parent. + cell->line(l)->child(c)->set_user_flag(); + } + } + } + } + + + // We loop through all cells and add cell to the + // map for the dofs that it immediately touches + // and then account for all the other dofs of + // which it is a part, mainly the ones that must + // be added on account of adaptivity hanging node + // constraints. + typename DoFHandlerType::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + // Need to loop through all cells that could + // be in the patch of dofs on locally_owned + // cells including ghost cells + if (cell->is_artificial() == false) + { + const unsigned int n_dofs_per_cell = cell->get_fe().dofs_per_cell; + local_dof_indices.resize(n_dofs_per_cell); + + // Take care of adding cell pointer to each + // dofs that exists on cell. + cell->get_dof_indices(local_dof_indices); + for (unsigned int i=0; i< n_dofs_per_cell; ++i ) + dof_to_set_of_cells_map[local_dof_indices[i]].insert(cell); + + // In the case of the adjacent cell (over + // faces or edges) being more refined, we + // want to add all of the children to the + // patch since the support function at that + // dof could be non-zero along that entire + // face (or line). + + // Take care of dofs on neighbor faces + for (unsigned int f=0; f::faces_per_cell; ++f) + { + if (cell->face(f)->has_children()) + { + for (unsigned int c=0; cface(f)->n_children(); ++c) + { + // Add cell to dofs of all subfaces + // + // *-------------------*----------*---------* + // | | add cell | | + // | |<- to dofs| | + // | |of subface| | + // | cell *----------*---------* + // | | add cell | | + // | |<- to dofs| | + // | |of subface| | + // *-------------------*----------*---------* + // + Assert (cell->face(f)->child(c)->has_children() == false, ExcInternalError()); + + const unsigned int n_dofs_per_face = cell->get_fe().dofs_per_face; + local_face_dof_indices.resize(n_dofs_per_face); + + cell->face(f)->child(c)->get_dof_indices(local_face_dof_indices); + for (unsigned int i=0; i< n_dofs_per_face; ++i ) + dof_to_set_of_cells_map[local_face_dof_indices[i]].insert(cell); + } + } + else if ((cell->face(f)->at_boundary() == false) && (cell->neighbor_is_coarser(f))) + { + + // Add cell to dofs of parent face and all + // child faces of parent face + // + // *-------------------*----------*---------* + // | | | | + // | | cell | | + // | add cell | | | + // | to dofs -> *----------*---------* + // | of parent | add cell | | + // | face |<- to dofs| | + // | |of subface| | + // *-------------------*----------*---------* + // + + // Add cell to all dofs of parent face + std::pair neighbor_face_no_subface_no = cell->neighbor_of_coarser_neighbor(f); + unsigned int face_no = neighbor_face_no_subface_no.first; + unsigned int subface = neighbor_face_no_subface_no.second; + + const unsigned int n_dofs_per_face = cell->get_fe().dofs_per_face; + local_face_dof_indices.resize(n_dofs_per_face); + + cell->neighbor(f)->face(face_no)->get_dof_indices(local_face_dof_indices); + for (unsigned int i=0; i< n_dofs_per_face; ++i ) + dof_to_set_of_cells_map[local_face_dof_indices[i]].insert(cell); + + // Add cell to all dofs of children of + // parent face + for (unsigned int c=0; cneighbor(f)->face(face_no)->n_children(); ++c) + { + if (c != subface) // don't repeat work on dofs of original cell + { + const unsigned int n_dofs_per_face = cell->get_fe().dofs_per_face; + local_face_dof_indices.resize(n_dofs_per_face); + + Assert (cell->neighbor(f)->face(face_no)->child(c)->has_children() == false, ExcInternalError()); + cell->neighbor(f)->face(face_no)->child(c)->get_dof_indices(local_face_dof_indices); + for (unsigned int i=0; i::lines_per_cell; ++l) + { + if (cell->line(l)->has_children()) + { + for (unsigned int c=0; cline(l)->n_children(); ++c) + { + Assert (cell->line(l)->child(c)->has_children() == false, ExcInternalError()); + + // dofs_per_line returns number of dofs + // on line not including the vertices of the line. + const unsigned int n_dofs_per_line = 2*cell->get_fe().dofs_per_vertex + + cell->get_fe().dofs_per_line; + local_line_dof_indices.resize(n_dofs_per_line); + + cell->line(l)->child(c)->get_dof_indices(local_line_dof_indices); + for (unsigned int i=0; iline(l)->user_flag_set() == true) + { + typename DoFHandlerType::line_iterator parent_line = lines_to_parent_lines_map[cell->line(l)]; + Assert (parent_line->has_children(), ExcInternalError() ); + + // dofs_per_line returns number of dofs + // on line not including the vertices of the line. + const unsigned int n_dofs_per_line = 2*cell->get_fe().dofs_per_vertex + + cell->get_fe().dofs_per_line; + local_line_dof_indices.resize(n_dofs_per_line); + + parent_line->get_dof_indices(local_line_dof_indices); + for (unsigned int i=0; in_children(); ++c) + { + Assert (parent_line->child(c)->has_children() == false, ExcInternalError()); + + const unsigned int n_dofs_per_line = 2*cell->get_fe().dofs_per_vertex + + cell->get_fe().dofs_per_line; + local_line_dof_indices.resize(n_dofs_per_line); + + parent_line->child(c)->get_dof_indices(local_line_dof_indices); + for (unsigned int i=0; iis_locally_owned() + }// for cells + + + if (DoFHandlerType::dimension == 3) + { + // finally, restore user flags that were changed above + // to when we constructed the pointers to parent of lines + // Since dof_handler is const, we must leave it unchanged. + const_cast &>(dof_handler.get_triangulation()).load_user_flags (user_flags); + } + + // Finally, we copy map of sets to + // map of vectors using the std::vector::assign() function + std::map< types::global_dof_index, std::vector > dof_to_cell_patches; + + typename std::map >::iterator + it = dof_to_set_of_cells_map.begin(), + it_end = dof_to_set_of_cells_map.end(); + for ( ; it!=it_end; ++it) + dof_to_cell_patches[it->first].assign( it->second.begin(), it->second.end() ); + + return dof_to_cell_patches; + } + + /* + * Internally used in collect_periodic_faces + */ + template + void + match_periodic_face_pairs + (std::set > &pairs1, + std::set::type, unsigned int> > &pairs2, + const int direction, + std::vector > &matched_pairs, + const dealii::Tensor<1,CellIterator::AccessorType::space_dimension> &offset, + const FullMatrix &matrix) + { + static const int space_dim = CellIterator::AccessorType::space_dimension; + (void)space_dim; + Assert (0<=direction && direction orientation; + typedef typename std::set + >::const_iterator PairIterator; + for (PairIterator it1 = pairs1.begin(); it1 != pairs1.end(); ++it1) + { + for (PairIterator it2 = pairs2.begin(); it2 != pairs2.end(); ++it2) + { + const CellIterator cell1 = it1->first; + const CellIterator cell2 = it2->first; + const unsigned int face_idx1 = it1->second; + const unsigned int face_idx2 = it2->second; + if (GridTools::orthogonal_equality(orientation, + cell1->face(face_idx1), + cell2->face(face_idx2), + direction, offset, + matrix)) + { + // We have a match, so insert the matching pairs and + // remove the matched cell in pairs2 to speed up the + // matching: + const PeriodicFacePair matched_face = + { + {cell1, cell2}, + {face_idx1, face_idx2}, + orientation, + matrix + }; + matched_pairs.push_back(matched_face); + pairs2.erase(it2); + ++n_matches; + break; + } + } + } + + //Assure that all faces are matched + AssertThrow (n_matches == pairs1.size() && pairs2.size() == 0, + ExcMessage ("Unmatched faces on periodic boundaries")); + } + + + + template + void + collect_periodic_faces + (const MeshType &mesh, + const types::boundary_id b_id, + const int direction, + std::vector > &matched_pairs, + const Tensor<1,MeshType::space_dimension> &offset, + const FullMatrix &matrix) + { + static const int dim = MeshType::dimension; + static const int space_dim = MeshType::space_dimension; + (void)dim; + (void)space_dim; + Assert (0<=direction && direction > pairs1; + std::set > pairs2; + + for (typename MeshType::cell_iterator cell = mesh.begin(0); + cell != mesh.end(0); ++cell) + { + const typename MeshType::face_iterator face_1 = cell->face(2*direction); + const typename MeshType::face_iterator face_2 = cell->face(2*direction+1); + + if (face_1->at_boundary() && face_1->boundary_id() == b_id) + { + const std::pair pair1 + = std::make_pair(cell, 2*direction); + pairs1.insert(pair1); + } + + if (face_2->at_boundary() && face_2->boundary_id() == b_id) + { + const std::pair pair2 + = std::make_pair(cell, 2*direction+1); + pairs2.insert(pair2); + } + } + + Assert (pairs1.size() == pairs2.size(), + ExcMessage ("Unmatched faces on periodic boundaries")); + + Assert (pairs1.size() > 0, + ExcMessage("No new periodic face pairs have been found. " + "Are you sure that you've selected the correct boundary " + "id's and that the coarsest level mesh is colorized?")); + +#ifdef DEBUG + const unsigned int size_old = matched_pairs.size(); +#endif + + // and call match_periodic_face_pairs that does the actual matching: + match_periodic_face_pairs(pairs1, pairs2, direction, matched_pairs, offset, + matrix); + +#ifdef DEBUG + //check for standard orientation + const unsigned int size_new = matched_pairs.size(); + for (unsigned int i = size_old; i < size_new; ++i) + { + Assert(matched_pairs[i].orientation == 1, + ExcMessage("Found a face match with non standard orientation. " + "This function is only suitable for meshes with cells " + "in default orientation")); + } +#endif + } + + + + template + void + collect_periodic_faces + (const MeshType &mesh, + const types::boundary_id b_id1, + const types::boundary_id b_id2, + const int direction, + std::vector > &matched_pairs, + const Tensor<1,MeshType::space_dimension> &offset, + const FullMatrix &matrix) + { + static const int dim = MeshType::dimension; + static const int space_dim = MeshType::space_dimension; + (void)dim; + (void)space_dim; + Assert (0<=direction && direction > pairs1; + std::set > pairs2; + + for (typename MeshType::cell_iterator cell = mesh.begin(0); + cell != mesh.end(0); ++cell) + { + for (unsigned int i = 0; i < GeometryInfo::faces_per_cell; ++i) + { + const typename MeshType::face_iterator face = cell->face(i); + if (face->at_boundary() && face->boundary_id() == b_id1) + { + const std::pair pair1 + = std::make_pair(cell, i); + pairs1.insert(pair1); + } + + if (face->at_boundary() && face->boundary_id() == b_id2) + { + const std::pair pair2 + = std::make_pair(cell, i); + pairs2.insert(pair2); + } + } + } + + Assert (pairs1.size() == pairs2.size(), + ExcMessage ("Unmatched faces on periodic boundaries")); + + Assert (pairs1.size() > 0, + ExcMessage("No new periodic face pairs have been found. " + "Are you sure that you've selected the correct boundary " + "id's and that the coarsest level mesh is colorized?")); + + // and call match_periodic_face_pairs that does the actual matching: + match_periodic_face_pairs(pairs1, pairs2, direction, matched_pairs, offset, + matrix); + } + + + + /* + * Internally used in orthogonal_equality + * + * An orthogonal equality test for points: + * + * point1 and point2 are considered equal, if + * matrix.point1 + offset - point2 + * is parallel to the unit vector in + */ + template + inline bool orthogonal_equality (const Point &point1, + const Point &point2, + const int direction, + const Tensor<1,spacedim> &offset, + const FullMatrix &matrix) + { + Assert (0<=direction && direction distance; + + if (matrix.m() == spacedim) + for (int i = 0; i < spacedim; ++i) + for (int j = 0; j < spacedim; ++j) + distance(i) += matrix(i,j) * point1(j); + else + distance = point1; + + distance += offset - point2; + + for (int i = 0; i < spacedim; ++i) + { + // Only compare coordinate-components != direction: + if (i == direction) + continue; + + if (fabs(distance(i)) > 1.e-10) + return false; + } + + return true; + } + + + /* + * Internally used in orthogonal_equality + * + * A lookup table to transform vertex matchings to orientation flags of + * the form (face_orientation, face_flip, face_rotation) + * + * See the comment on the next function as well as the detailed + * documentation of make_periodicity_constraints and + * collect_periodic_faces for details + */ + template struct OrientationLookupTable {}; + + template <> struct OrientationLookupTable<1> + { + typedef std::array::vertices_per_face> MATCH_T; + static inline std::bitset<3> lookup (const MATCH_T &) + { + // The 1D case is trivial + return 1; // [true ,false,false] + } + }; + + template <> struct OrientationLookupTable<2> + { + typedef std::array::vertices_per_face> MATCH_T; + static inline std::bitset<3> lookup (const MATCH_T &matching) + { + // In 2D matching faces (=lines) results in two cases: Either + // they are aligned or flipped. We store this "line_flip" + // property somewhat sloppy as "face_flip" + // (always: face_orientation = true, face_rotation = false) + + static const MATCH_T m_tff = {{ 0, 1 }}; + if (matching == m_tff) return 1; // [true ,false,false] + static const MATCH_T m_ttf = {{ 1, 0 }}; + if (matching == m_ttf) return 3; // [true ,true ,false] + Assert(false, ExcInternalError()); + // what follows is dead code, but it avoids warnings about the lack + // of a return value + return 0; + } + }; + + template <> struct OrientationLookupTable<3> + { + typedef std::array::vertices_per_face> MATCH_T; + static inline std::bitset<3> lookup (const MATCH_T &matching) + { + // The full fledged 3D case. *Yay* + // See the documentation in include/deal.II/base/geometry_info.h + // as well as the actual implementation in source/grid/tria.cc + // for more details... + + static const MATCH_T m_tff = {{ 0, 1, 2, 3 }}; + if (matching == m_tff) return 1; // [true ,false,false] + static const MATCH_T m_tft = {{ 1, 3, 0, 2 }}; + if (matching == m_tft) return 5; // [true ,false,true ] + static const MATCH_T m_ttf = {{ 3, 2, 1, 0 }}; + if (matching == m_ttf) return 3; // [true ,true ,false] + static const MATCH_T m_ttt = {{ 2, 0, 3, 1 }}; + if (matching == m_ttt) return 7; // [true ,true ,true ] + static const MATCH_T m_fff = {{ 0, 2, 1, 3 }}; + if (matching == m_fff) return 0; // [false,false,false] + static const MATCH_T m_fft = {{ 2, 3, 0, 1 }}; + if (matching == m_fft) return 4; // [false,false,true ] + static const MATCH_T m_ftf = {{ 3, 1, 2, 0 }}; + if (matching == m_ftf) return 2; // [false,true ,false] + static const MATCH_T m_ftt = {{ 1, 0, 3, 2 }}; + if (matching == m_ftt) return 6; // [false,true ,true ] + Assert(false, ExcInternalError()); + // what follows is dead code, but it avoids warnings about the lack + // of a return value + return 0; + } + }; + + + + template + inline bool + orthogonal_equality (std::bitset<3> &orientation, + const FaceIterator &face1, + const FaceIterator &face2, + const int direction, + const Tensor<1,FaceIterator::AccessorType::space_dimension> &offset, + const FullMatrix &matrix) + { + Assert(matrix.m() == matrix.n(), + ExcMessage("The supplied matrix must be a square matrix")); + + static const int dim = FaceIterator::AccessorType::dimension; + + // Do a full matching of the face vertices: + + std:: + array::vertices_per_face> matching; + + std::set face2_vertices; + for (unsigned int i = 0; i < GeometryInfo::vertices_per_face; ++i) + face2_vertices.insert(i); + + for (unsigned int i = 0; i < GeometryInfo::vertices_per_face; ++i) + { + for (std::set::iterator it = face2_vertices.begin(); + it != face2_vertices.end(); + ++it) + { + if (orthogonal_equality(face1->vertex(i),face2->vertex(*it), + direction, offset, matrix)) + { + matching[i] = *it; + face2_vertices.erase(it); + break; // jump out of the innermost loop + } + } + } + + // And finally, a lookup to determine the ordering bitmask: + if (face2_vertices.empty()) + orientation = OrientationLookupTable::lookup(matching); + + return face2_vertices.empty(); + } + + + + template + inline bool + orthogonal_equality (const FaceIterator &face1, + const FaceIterator &face2, + const int direction, + const Tensor<1,FaceIterator::AccessorType::space_dimension> &offset, + const FullMatrix &matrix) + { + // Call the function above with a dummy orientation array + std::bitset<3> dummy; + return orthogonal_equality (dummy, face1, face2, direction, offset, matrix); + } + + + +} + + +#include "grid_tools_dof_handlers.inst" + + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/grid/grid_tools_dof_handlers.inst.in b/source/grid/grid_tools_dof_handlers.inst.in new file mode 100644 index 0000000000..1444b26d71 --- /dev/null +++ b/source/grid/grid_tools_dof_handlers.inst.in @@ -0,0 +1,242 @@ + +for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS) +{ + +#if deal_II_dimension <= deal_II_space_dimension + namespace GridTools \{ + + template + unsigned int + find_closest_vertex (const X &, + const Point &, + const std::vector &); + + template + unsigned int + find_closest_vertex (const Mapping&, + const X &, + const Point &, + const std::vector &); + + template + std::vector::type> + find_cells_adjacent_to_vertex(const X &, + const unsigned int); + + template + dealii::internal::ActiveCellIterator::type + find_active_cell_around_point (const X &, + const Point &, + const std::vector &); + + template + std::pair::type, Point > + find_active_cell_around_point (const Mapping &, + const X &, + const Point &, + const std::vector &); + + template + std::vector::type> + compute_active_cell_halo_layer (const X &, + const std::function::type&)> &); + + template + std::vector + compute_cell_halo_layer_on_level (const X &, + const std::function &, + const unsigned int); + + template + std::vector::type> + compute_ghost_cell_halo_layer (const X &); + + + template + std::vector::type> + compute_active_cell_layer_within_distance (const X &, + const std::function::type&)> &, + const double); + + + template + std::vector::type> + compute_ghost_cell_layer_within_distance (const X &, const double); + + + template + std::pair< Point, Point > + compute_bounding_box (const X &, + const std::function::type&)> &); + + + template + std::list > + get_finest_common_cells (const X &mesh_1, + const X &mesh_2); + + + template + bool + have_same_coarse_mesh (const X &mesh_1, + const X &mesh_2); + \} + +#endif +} + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) +{ +#if deal_II_dimension <= deal_II_space_dimension + namespace GridTools \{ + + template + std::pair::active_cell_iterator, + Point > + find_active_cell_around_point + (const hp::MappingCollection &, + const hp::DoFHandler &, + const Point &); + + \} +#endif +} + + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; Container : TRIANGULATION_AND_DOFHANDLER_TEMPLATES) +{ +#if deal_II_dimension <= deal_II_space_dimension + namespace GridTools \{ + + template + std::vector::active_cell_iterator> + get_patch_around_cell > + (const Container::active_cell_iterator &cell); + + template + std::vector< Container::cell_iterator> + get_cells_at_coarsest_common_level > ( + const std::vector< Container::active_cell_iterator> & patch_cells); + + template + void build_triangulation_from_patch > ( + const std::vector::active_cell_iterator> &patch, + Triangulation::dimension,Container::space_dimension> &local_triangulation, + std::map::active_cell_iterator, + Container::active_cell_iterator > &patch_to_global_tria_map); + + \} +#endif +} + + + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; Container : DOFHANDLER_TEMPLATES) +{ +#if deal_II_dimension <= deal_II_space_dimension + namespace GridTools \{ + + template + std::map< types::global_dof_index,std::vector::active_cell_iterator> > + get_dof_to_support_patch_map > + (Container &dof_handler); + + \} +#endif +} + + + +// instantiate the following functions only for the "sequential" containers. this +// is a misnomer here, however: the point is simply that we only instantiate +// these functions for certain *iterator* types, and the iterator types are +// the same for sequential and parallel containers; consequently, we get duplicate +// instantiation errors if we instantiate for *all* container types, rather than +// only the sequential ones +for (X : SEQUENTIAL_TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS) +{ +#if deal_II_dimension <= deal_II_space_dimension + namespace GridTools \{ + + template + bool orthogonal_equality (std::bitset<3> &, + const X::active_face_iterator&, + const X::active_face_iterator&, + const int, + const Tensor<1,deal_II_space_dimension> &, + const FullMatrix &); + + template + bool orthogonal_equality (std::bitset<3> &, + const X::face_iterator&, + const X::face_iterator&, + const int, + const Tensor<1,deal_II_space_dimension> &, + const FullMatrix &); + + template + bool orthogonal_equality (const X::active_face_iterator&, + const X::active_face_iterator&, + const int, + const Tensor<1,deal_II_space_dimension> &, + const FullMatrix &); + + template + bool orthogonal_equality (const X::face_iterator&, + const X::face_iterator&, + const int, + const Tensor<1,deal_II_space_dimension> &, + const FullMatrix &); + + template + void collect_periodic_faces (const X &, + const types::boundary_id, + const types::boundary_id, + const int, + std::vector > &, + const Tensor<1,X::space_dimension> &, + const FullMatrix &); + + template + void collect_periodic_faces (const X &, + const types::boundary_id, + const int, + std::vector > &, + const Tensor<1,X::space_dimension> &, + const FullMatrix &); + + \} +#endif +} + +// TODO the text above the last instantiation block implies that this should not be necessary... is it? +for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS) +{ +#if deal_II_dimension <= deal_II_space_dimension +#if deal_II_dimension >= 2 + + namespace GridTools \{ + template + void + collect_periodic_faces > + (const parallel::distributed::Triangulation &, + const types::boundary_id, + const types::boundary_id, + const int, + std::vector::cell_iterator> > &, + const Tensor<1,parallel::distributed::Triangulation::space_dimension> &, + const FullMatrix &); + + template + void + collect_periodic_faces > + (const parallel::distributed::Triangulation &, + const types::boundary_id, + const int, + std::vector::cell_iterator> > &, + const Tensor<1,parallel::distributed::Triangulation::space_dimension> &, + const FullMatrix &); + \} +#endif +#endif +}