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 <typename ITERATOR>
- void
- identify_periodic_vertices_recursively(const GridTools::PeriodicFacePair<ITERATOR> &periodic,
- std::vector<unsigned int> &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<ITERATOR> periodic_child = periodic;
-
- // find appropriate pairs of child elements
- for (unsigned int cf=0; cf<periodic.cell[0]->face(periodic.face_idx[0])->n_children(); ++cf)
- {
- const unsigned int child_index_0 =
- GeometryInfo<dim>::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<dim>::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<GeometryInfo<dim-1>::vertices_per_cell; ++v)
- {
- // take possible non-standard orientation of face on cell[0] into
- // account
- const unsigned int vface0 =
- GeometryInfo<dim>::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 <int dim, int spacedim>
tria.save_refine_flags (flags_before[1]);
std::vector<unsigned int> topological_vertex_numbering(tria.n_vertices());
- std::vector<unsigned int> topological_vertex_numbering_cmp(tria.n_vertices());
for (unsigned int i=0; i<topological_vertex_numbering.size(); ++i)
- {
- topological_vertex_numbering[i] = i;
- topological_vertex_numbering_cmp[i] = i;
- }
- //std::cout << std::endl;
- for (unsigned int i=0; i<periodic_face_pairs_level_0.size(); ++i)
- {
- identify_periodic_vertices_recursively(periodic_face_pairs_level_0[i],
- topological_vertex_numbering_cmp);
- }
+ topological_vertex_numbering[i] = i;
// combine 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
topological_vertex_numbering[cell_1->face(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;*/
}
}
}
ExcInternalError());
}
- for (unsigned int i=0; i<topological_vertex_numbering_cmp.size(); ++i)
- {
- const unsigned int j = topological_vertex_numbering_cmp[i];
- if (j != i)
- Assert(topological_vertex_numbering_cmp[j] == j,
- ExcInternalError());
- }
-
- if (topological_vertex_numbering != topological_vertex_numbering_cmp)
- {
- for (unsigned int i=0; i<topological_vertex_numbering.size(); ++i)
- {
- if (topological_vertex_numbering[i] != topological_vertex_numbering_cmp[i])
- std::cout << i
- << " old: " << topological_vertex_numbering_cmp[i]
- << " new: " << topological_vertex_numbering[i]
- << std::endl;
- }
- Assert( false, ExcInternalError());
- }
-
// this code is replicated from grid/tria.cc but using an indirection
// for periodic boundary conditions
- namespace
- {
- /**
- * ensures 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
- * one is active
- */
- template <typename ITERATOR>
- void
- mark_periodic_vertices_recursively(const GridTools::PeriodicFacePair<ITERATOR> &periodic,
- const int target_level,
- std::vector<bool> &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<ITERATOR> periodic_child = periodic;
-
- // find appropriate pairs of child elements
- for (unsigned int cf=0; cf<periodic.cell[0]->face(periodic.face_idx[0])->n_children(); ++cf)
- {
- const unsigned int child_index_0 =
- GeometryInfo<dim>::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<dim>::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<GeometryInfo<dim-1>::vertices_per_cell; ++v)
- {
- // take possible non-standard orientation of face on cell[0] into
- // account
- const unsigned int vface0 =
- GeometryInfo<dim>::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 <typename ITERATOR>
- void
- set_periodic_ghost_neighbors_recursively(const GridTools::PeriodicFacePair<ITERATOR> &periodic,
- const int target_level,
- std::map<unsigned int, std::set<dealii::types::subdomain_id> > &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<ITERATOR> periodic_child = periodic;
-
- // find appropriate pairs of child elements
- for (unsigned int cf=0; cf<periodic.cell[0]->face(periodic.face_idx[0])->n_children(); ++cf)
- {
- unsigned int c=0;
- for (; c<periodic.cell[0]->n_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 (; c<periodic.cell[1]->n_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<GeometryInfo<ITERATOR::AccessorType::dimension-1>::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
vertices_with_ghost_neighbors[cell->vertex_index(v)]
.insert (cell->level_subdomain_id());
- std::map<unsigned int, std::set<dealii::types::subdomain_id> >
- vertices_with_ghost_neighbors_cmp = vertices_with_ghost_neighbors;
-
-
- for (unsigned int i=0; i<this->periodic_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<cell_iterator, unsigned int>,
std::pair<std::pair<cell_iterator,unsigned int>, std::bitset<3> > >::const_iterator it;
}
}
}
-
- if (vertices_with_ghost_neighbors != vertices_with_ghost_neighbors_cmp)
- {
- std::cout << "Print vertices_with_ghost_neighbors: " << std::endl;
- std::map<unsigned int, std::set<dealii::types::subdomain_id> >::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<dealii::types::subdomain_id>::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<dealii::types::subdomain_id>::const_iterator it2;
- for (it2 = it->second.begin(); it2!=it->second.end(); ++it2)
- std::cout << *it2 << " ";
- std::cout << std::endl;
- }
- Assert(false, ExcInternalError());
- }
}
for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
marked_vertices[cell->vertex_index(v)] = true;
- std::vector<bool> marked_vertices_cmp = marked_vertices;
-
- for (unsigned int i=0; i<this->periodic_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
}
}
- Assert(marked_vertices == marked_vertices_cmp, ExcInternalError());
return marked_vertices;
}
projects / dealii.git / commitdiff