From: Daniel Arndt Date: Tue, 22 Mar 2016 15:35:35 +0000 (+0100) Subject: Use only periodioc face_map for the multigrid structures in distributed/tria.h X-Git-Tag: v8.5.0-rc1~1151^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b4ba1602dfdfe72579542651b86703c00b755fd1;p=dealii.git Use only periodioc face_map for the multigrid structures in distributed/tria.h --- diff --git a/source/distributed/tria.cc b/source/distributed/tria.cc index 2e052d49e2..605000a0bc 100644 --- a/source/distributed/tria.cc +++ b/source/distributed/tria.cc @@ -3124,91 +3124,6 @@ namespace parallel namespace { - // this function combines vertices that have different locations (and - // thus, different vertex_index) but represent the same topological - // entity over periodic boundaries. The vector - // topological_vertex_numbering contains a linear map from 0 to - // n_vertices at input and at output relates periodic vertices with only - // one vertex index. The output is used to always identify the same - // vertex according to the periodicity, e.g. when finding the maximum - // cell level around a vertex. - // - // Example: On a 3D cell with vertices numbered from 0 to 7 and periodic - // boundary conditions in x direction, the vector - // topological_vertex_numbering will contain the numbers - // {0,0,2,2,4,4,6,6} (because the vertex pairs {0,1}, {2,3}, {4,5}, - // {6,7} belong together, respectively). If periodicity is set in x and - // z direction, the output is {0,0,2,2,0,0,2,2}, and if periodicity is - // in all directions, the output is simply {0,0,0,0,0,0,0,0}. - template - void - identify_periodic_vertices_recursively(const GridTools::PeriodicFacePair &periodic, - std::vector &topological_vertex_numbering) - { - const unsigned int dim = ITERATOR::AccessorType::dimension; - - // for hanging nodes we will consider all necessary coupling already - // on the parent level, so we just need to consider neighbors of the - // same level - if (periodic.cell[0]->has_children() && - periodic.cell[1]->has_children()) - { - // copy orientations etc. from parent to child - GridTools::PeriodicFacePair periodic_child = periodic; - - // find appropriate pairs of child elements - for (unsigned int cf=0; cfface(periodic.face_idx[0])->n_children(); ++cf) - { - const unsigned int child_index_0 = - GeometryInfo::child_cell_on_face(periodic.cell[0]->refinement_case(), - periodic.face_idx[0], cf, - periodic.orientation[0], - periodic.orientation[1], - periodic.orientation[2]); - periodic_child.cell[0] = periodic.cell[0]->child(child_index_0); - periodic_child.face_idx[0] = periodic.face_idx[0]; - - // the second face is in standard orientation in terms of the - // periodic face pair - const unsigned int child_index_1 = - GeometryInfo::child_cell_on_face(periodic.cell[1]->refinement_case(), - periodic.face_idx[1], cf); - - periodic_child.cell[1] = periodic.cell[1]->child(child_index_1); - periodic_child.face_idx[1] = periodic.face_idx[1]; - - // recursive call into children - identify_periodic_vertices_recursively (periodic_child, - topological_vertex_numbering); - } - } - - for (unsigned int v=0; v::vertices_per_cell; ++v) - { - // take possible non-standard orientation of face on cell[0] into - // account - const unsigned int vface0 = - GeometryInfo::standard_to_real_face_vertex(v,periodic.orientation[0], - periodic.orientation[1], - periodic.orientation[2]); - const unsigned int vi0 = topological_vertex_numbering[periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(vface0)]; - const unsigned int vi1 = topological_vertex_numbering[periodic.cell[1]->face(periodic.face_idx[1])->vertex_index(v)]; - const unsigned int min_index = std::min(vi0, vi1); - topological_vertex_numbering[periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(vface0)] - = topological_vertex_numbering[periodic.cell[1]->face(periodic.face_idx[1])->vertex_index(v)] - = min_index; - - /* std::cout << "Old identified vertices " << periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(vface0) - << " at " << periodic.cell[0]->face(periodic.face_idx[0])->vertex(vface0) - << " and " << periodic.cell[1]->face(periodic.face_idx[1])->vertex_index(v) - << " at " << periodic.cell[1]->face(periodic.face_idx[1])->vertex(v) - << std::endl;*/ - - } - } - - - // ensures the 2:1 mesh balance for periodic boundary conditions in the // artificial cell layer (the active cells are taken care of by p4est) template @@ -3223,18 +3138,8 @@ namespace parallel tria.save_refine_flags (flags_before[1]); std::vector topological_vertex_numbering(tria.n_vertices()); - std::vector topological_vertex_numbering_cmp(tria.n_vertices()); for (unsigned int i=0; iface(face_no_1)->vertex_index(vface0)] = topological_vertex_numbering[cell_2->face(face_no_2)->vertex_index(v)] = min_index; - /* std::cout << "Newly identified vertices " << cell_1->face(face_no_1)->vertex_index(vface0) - << " at " << cell_1->face(face_no_1)->vertex(vface0) - << " and " << cell_2->face(face_no_2)->vertex_index(v) - << " at " << cell_2->face(face_no_2)->vertex(v) - << std::endl;*/ } } } @@ -3295,27 +3195,6 @@ namespace parallel ExcInternalError()); } - for (unsigned int i=0; i - void - mark_periodic_vertices_recursively(const GridTools::PeriodicFacePair &periodic, - const int target_level, - std::vector &active_vertices_on_level) - { - if (periodic.cell[0]->level() > target_level) - return; - - const unsigned int dim = ITERATOR::AccessorType::dimension; - // for hanging nodes there is nothing to do since we are interested in - // the connections on the same level... - if (periodic.cell[0]->level() < target_level && - periodic.cell[0]->has_children() && - periodic.cell[1]->has_children()) - { - // copy orientations etc. from parent to child - GridTools::PeriodicFacePair periodic_child = periodic; - - // find appropriate pairs of child elements - for (unsigned int cf=0; cfface(periodic.face_idx[0])->n_children(); ++cf) - { - const unsigned int child_index_0 = - GeometryInfo::child_cell_on_face(periodic.cell[0]->refinement_case(), - periodic.face_idx[0], cf, - periodic.orientation[0], - periodic.orientation[1], - periodic.orientation[2]); - periodic_child.cell[0] = periodic.cell[0]->child(child_index_0); - periodic_child.face_idx[0] = periodic.face_idx[0]; - - // the second face is in standard orientation in terms of the - // periodic face pair - const unsigned int child_index_1 = - GeometryInfo::child_cell_on_face(periodic.cell[1]->refinement_case(), - periodic.face_idx[1], cf); - - periodic_child.cell[1] = periodic.cell[1]->child(child_index_1); - periodic_child.face_idx[1] = periodic.face_idx[1]; - - // recursive call into children - mark_periodic_vertices_recursively (periodic_child, target_level, - active_vertices_on_level); - } - return; - } - - if (periodic.cell[0]->level() != target_level) - return; - - for (unsigned int v=0; v::vertices_per_cell; ++v) - { - // take possible non-standard orientation of face on cell[0] into - // account - const unsigned int vface0 = - GeometryInfo::standard_to_real_face_vertex(v,periodic.orientation[0], - periodic.orientation[1], - periodic.orientation[2]); - if (active_vertices_on_level[periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(vface0)] || - active_vertices_on_level[periodic.cell[1]->face(periodic.face_idx[1])->vertex_index(v)]) - active_vertices_on_level[periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(vface0)] - = active_vertices_on_level[periodic.cell[1]->face(periodic.face_idx[1])->vertex_index(v)] - = true; - } - } - - - - /** - * ensures that always both vertices over a periodic face are identified - * together - */ - template - void - set_periodic_ghost_neighbors_recursively(const GridTools::PeriodicFacePair &periodic, - const int target_level, - std::map > &vertices_with_ghost_neighbors) - { - if (periodic.cell[0]->level() > target_level) - return; - - // for hanging nodes there is nothing to do since we are interested in - // the connections on the same level... - if (periodic.cell[0]->level() < target_level && - periodic.cell[0]->has_children() && - periodic.cell[1]->has_children()) - { - // copy orientations etc. from parent to child - GridTools::PeriodicFacePair periodic_child = periodic; - - // find appropriate pairs of child elements - for (unsigned int cf=0; cfface(periodic.face_idx[0])->n_children(); ++cf) - { - unsigned int c=0; - for (; cn_children(); ++c) - { - if (periodic.cell[0]->child(c)->face(periodic.face_idx[0]) == - periodic.cell[0]->face(periodic.face_idx[0])->child(cf)) - break; - } - Assert(c < periodic.cell[0]->n_children(), - ExcMessage("Face child not found")); - periodic_child.cell[0] = periodic.cell[0]->child(c); - periodic_child.face_idx[0] = periodic.face_idx[0]; - - c=0; - for (; cn_children(); ++c) - { - if (periodic.cell[1]->child(c)->face(periodic.face_idx[1]) == - periodic.cell[1]->face(periodic.face_idx[1])->child(cf)) - break; - } - Assert(c < periodic.cell[1]->n_children(), - ExcMessage("Face child not found")); - periodic_child.cell[1] = periodic.cell[1]->child(c); - periodic_child.face_idx[1] = periodic.face_idx[1]; - - // recursive call into children - set_periodic_ghost_neighbors_recursively (periodic_child, target_level, - vertices_with_ghost_neighbors); - } - return; - } - - if (periodic.cell[0]->level() != target_level) - return; - - // TODO: fix non-standard orientation - Assert(periodic.orientation[0] == true && - periodic.orientation[1] == false && - periodic.orientation[2] == false, - ExcNotImplemented()); - - for (unsigned int v=0; v::vertices_per_cell; ++v) - { - const unsigned int - idx0 = periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(v), - idx1 = periodic.cell[1]->face(periodic.face_idx[1])->vertex_index(v); - if (vertices_with_ghost_neighbors.find(idx0) != - vertices_with_ghost_neighbors.end()) - vertices_with_ghost_neighbors[idx1]. - insert(vertices_with_ghost_neighbors[idx0].begin(), - vertices_with_ghost_neighbors[idx0].end()); - if (vertices_with_ghost_neighbors.find(idx1) != - vertices_with_ghost_neighbors.end()) - vertices_with_ghost_neighbors[idx0]. - insert(vertices_with_ghost_neighbors[idx1].begin(), - vertices_with_ghost_neighbors[idx1].end()); - } - } - } - - - /** * Determine the neighboring subdomains that are adjacent to each vertex. * This is achieved via the p4est_iterate/p8est_iterate tool @@ -4785,14 +4503,6 @@ namespace parallel vertices_with_ghost_neighbors[cell->vertex_index(v)] .insert (cell->level_subdomain_id()); - std::map > - vertices_with_ghost_neighbors_cmp = vertices_with_ghost_neighbors; - - - for (unsigned int i=0; iperiodic_face_pairs_level_0.size(); ++i) - set_periodic_ghost_neighbors_recursively(this->periodic_face_pairs_level_0[i], - level, vertices_with_ghost_neighbors_cmp); - //now for the vertices on periodic faces typename std::map, std::pair, std::bitset<3> > >::const_iterator it; @@ -4826,31 +4536,6 @@ namespace parallel } } } - - if (vertices_with_ghost_neighbors != vertices_with_ghost_neighbors_cmp) - { - std::cout << "Print vertices_with_ghost_neighbors: " << std::endl; - std::map >::const_iterator it; - for (it = vertices_with_ghost_neighbors.begin(); it!=vertices_with_ghost_neighbors.end(); ++it) - { - std::cout << "Vertex " << it->first << " is on subdomains "; - std::set::const_iterator it2; - for (it2 = it->second.begin(); it2!=it->second.end(); ++it2) - std::cout << *it2 << " "; - std::cout << std::endl; - } - - std::cout << "Print vertices_with_ghost_neighbors_cmp: " << std::endl; - for (it = vertices_with_ghost_neighbors_cmp.begin(); it!=vertices_with_ghost_neighbors_cmp.end(); ++it) - { - std::cout << "Vertex " << it->first << " is on subdomains "; - std::set::const_iterator it2; - for (it2 = it->second.begin(); it2!=it->second.end(); ++it2) - std::cout << *it2 << " "; - std::cout << std::endl; - } - Assert(false, ExcInternalError()); - } } @@ -4870,12 +4555,6 @@ namespace parallel for (unsigned int v=0; v::vertices_per_cell; ++v) marked_vertices[cell->vertex_index(v)] = true; - std::vector marked_vertices_cmp = marked_vertices; - - for (unsigned int i=0; iperiodic_face_pairs_level_0.size(); ++i) - mark_periodic_vertices_recursively(this->periodic_face_pairs_level_0[i], - level, marked_vertices_cmp); - /** * ensure that if one of the two vertices on a periodic face is marked * as active (i.e., belonging to an owned level cell), also the other @@ -4911,7 +4590,6 @@ namespace parallel } } - Assert(marked_vertices == marked_vertices_cmp, ExcInternalError()); return marked_vertices; }