* grid cells are badly ordered this may mean that individual parts of
* the forest stored on a local machine may be split across coarse grid
* cells that are not geometrically close. Consequently, we apply a
- * Cuthill-McKee preordering to ensure that the part of the forest
- * stored by p4est is located on geometrically close coarse grid cells.
+ * hierarchical preordering according to
+ * SparsityTools::reorder_hierarchical() to ensure that the part of the
+ * forest stored by p4est is located on geometrically close coarse grid
+ * cells.
*/
std::vector<types::global_dof_index> coarse_cell_to_p4est_tree_permutation;
std::vector<types::global_dof_index> p4est_tree_to_coarse_cell_permutation;
+ /**
+ * If add_periodicity() is called, this variable stores the given
+ * periodic face pairs on level 0 for later access during the
+ * identification of ghost cells for the multigrid hierarchy.
+ */
+ std::vector<GridTools::PeriodicFacePair<cell_iterator> > periodic_face_pairs_level_0;
+
/**
* Return a pointer to the p4est tree that belongs to the given
* dealii_coarse_cell_index()
void attach_mesh_data();
/**
- * fills a map that, for each vertex, lists all the processors whose
- * subdomains are adjacent to that vertex. Used by
+ * Fills a map that, for each vertex, lists all the processors whose
+ * subdomains are adjacent to that vertex. Used by
* DoFHandler::Policy::ParallelDistributed.
*/
void
(std::map<unsigned int, std::set<dealii::types::subdomain_id> >
&vertices_with_ghost_neighbors);
+ /**
+ * Fills a map that, for each vertex, lists all the processors whose
+ * subdomains are adjacent to that vertex on the given level for the
+ * multigrid hierarchy. Used by DoFHandler::Policy::ParallelDistributed.
+ */
+ void
+ fill_level_vertices_with_ghost_neighbors
+ (const unsigned int level,
+ std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ &vertices_with_ghost_neighbors);
+
+ /**
+ * This method returns a bit vector of length tria.n_vertices()
+ * indicating the locally active vertices on a level, i.e., the vertices
+ * touched by the locally owned level cells for use in geometric
+ * multigrid (possibly including the vertices due to periodic boundary
+ * conditions) are marked by true.
+ *
+ * Used by DoFHandler::Policy::ParallelDistributed.
+ */
+ std::vector<bool>
+ mark_locally_active_vertices_on_level(const unsigned int level) const;
+
template <int, int> friend class dealii::internal::DoFHandler::Policy::ParallelDistributed;
};
(std::map<unsigned int, std::set<dealii::types::subdomain_id> >
&vertices_with_ghost_neighbors);
+ /**
+ * Like above, this method, which is only implemented for dim = 2 or 3,
+ * needs a stub because it is used in dof_handler_policy.cc
+ */
+ void
+ fill_level_vertices_with_ghost_neighbors
+ (const unsigned int level,
+ std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ &vertices_with_ghost_neighbors);
+
+ /**
+ * Like above, this method, which is only implemented for dim = 2 or 3,
+ * needs a stub because it is used in dof_handler_policy.cc
+ */
+ std::vector<bool>
+ mark_locally_active_vertices_on_level(const unsigned int level) const;
+
};
}
}
coarse_cell_to_p4est_tree_permutation.resize (0);
p4est_tree_to_coarse_cell_permutation.resize (0);
+ periodic_face_pairs_level_0.clear();
+
dealii::Triangulation<dim,spacedim>::clear ();
this->update_number_cache ();
// -1. Note that we do not fill other cells we could figure out the
// same way, because we might accidentally set an id for a cell that
// is not a ghost cell.
- for (unsigned int lvl=this->n_levels(); lvl>0;)
+ for (unsigned int lvl=this->n_levels(); lvl>0; )
{
--lvl;
for (typename Triangulation<dim,spacedim>::cell_iterator cell = this->begin(lvl); cell!=this->end(lvl); ++cell)
//step 2: make sure all the neighbors to our level_cells exist. Need
//to look up in p4est...
std::vector<std::vector<bool> > marked_vertices(this->n_levels());
-
- for (unsigned int lvl=0; lvl<this->n_levels(); ++lvl)
- {
- marked_vertices[lvl].resize(this->n_vertices(), false);
-
- for (typename dealii::Triangulation<dim,spacedim>::cell_iterator
- cell = this->begin(lvl);
- cell != this->end(lvl); ++cell)
- if (cell->level_subdomain_id() == this->locally_owned_subdomain())
- for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
- marked_vertices[lvl][cell->vertex_index(v)] = true;
- }
-
-
+ for (unsigned int lvl=0; lvl < this->n_levels(); ++lvl)
+ marked_vertices[lvl] = mark_locally_active_vertices_on_level(lvl);
for (typename Triangulation<dim,spacedim>::cell_iterator cell = this->begin(0); cell!=this->end(0); ++cell)
{
}
+
+ 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;
+
+ // 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
+ mark_periodic_vertices_recursively (periodic_child, target_level,
+ active_vertices_on_level);
+ }
+ 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)
+ if (active_vertices_on_level[periodic.cell[0]->face(periodic.face_idx[0])->vertex_index(v)] ||
+ 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(v)]
+ = 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
}
+
+ /**
+ * Determine the neighboring subdomains that are adjacent to each vertex
+ * on the given multigrid level
+ */
+ template <int dim, int spacedim>
+ void
+ Triangulation<dim,spacedim>::
+ fill_level_vertices_with_ghost_neighbors
+ (const unsigned int level,
+ std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ &vertices_with_ghost_neighbors)
+ {
+ const std::vector<bool> locally_active_vertices =
+ mark_locally_active_vertices_on_level(level);
+ for (cell_iterator cell = this->begin(level); cell != this->end(level); ++cell)
+ if (cell->level_subdomain_id() != dealii::numbers::artificial_subdomain_id
+ && cell->level_subdomain_id() != this->locally_owned_subdomain())
+ for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+ if (locally_active_vertices[cell->vertex_index(v)])
+ vertices_with_ghost_neighbors[cell->vertex_index(v)]
+ .insert (cell->level_subdomain_id());
+
+ for (unsigned int i=0; i<periodic_face_pairs_level_0.size(); ++i)
+ set_periodic_ghost_neighbors_recursively(periodic_face_pairs_level_0[i],
+ level, vertices_with_ghost_neighbors);
+ }
+
+
+
+ template<int dim, int spacedim>
+ std::vector<bool>
+ Triangulation<dim,spacedim>
+ ::mark_locally_active_vertices_on_level (const unsigned int level) const
+ {
+ Assert (dim>1, ExcNotImplemented());
+
+ std::vector<bool> marked_vertices(this->n_vertices(), false);
+ for (cell_iterator cell = this->begin(level); cell != this->end(level); ++cell)
+ if (cell->level_subdomain_id() == this->locally_owned_subdomain())
+ for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+ marked_vertices[cell->vertex_index(v)] = true;
+
+ for (unsigned int i=0; i<periodic_face_pairs_level_0.size(); ++i)
+ mark_periodic_vertices_recursively(periodic_face_pairs_level_0[i],
+ level, marked_vertices);
+
+ return marked_vertices;
+ }
+
+
+
template<int dim, int spacedim>
void
Triangulation<dim,spacedim>::add_periodicity
AssertThrow (false, ExcInternalError());
}
+ periodic_face_pairs_level_0.insert(periodic_face_pairs_level_0.end(),
+ periodicity_vector.begin(),
+ periodicity_vector.end());
+
#else
Assert(false, ExcMessage ("Need p4est version >= 0.3.4.1!"));
#endif
{
Assert (false, ExcNotImplemented());
}
+
+ template <int spacedim>
+ void
+ Triangulation<1,spacedim>::
+ fill_level_vertices_with_ghost_neighbors
+ (const unsigned int /*level*/,
+ std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ &/*vertices_with_ghost_neighbors*/)
+ {
+ Assert (false, ExcNotImplemented());
+ }
+
+ template <int spacedim>
+ std::vector<bool>
+ Triangulation<1,spacedim>::
+ mark_locally_active_vertices_on_level (const unsigned int) const
+ {
+ Assert (false, ExcNotImplemented());
+ return std::vector<bool>();
+ }
}
}
}
- /**
- * Return a vector in which,
- * for every vertex index, we
- * mark whether a locally owned
- * cell is adjacent.
- */
- template <int dim, int spacedim>
- std::vector<bool>
- mark_locally_active_vertices (const parallel::distributed::Triangulation<dim,spacedim> &triangulation)
- {
- std::vector<bool> locally_active_vertices (triangulation.n_vertices(),
- false);
-
- for (typename dealii::Triangulation<dim,spacedim>::active_cell_iterator
- cell = triangulation.begin_active();
- cell != triangulation.end(); ++cell)
- if (cell->subdomain_id() == triangulation.locally_owned_subdomain())
- for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
- locally_active_vertices[cell->vertex_index(v)] = true;
-
- return locally_active_vertices;
- }
-
- /**
- * Return a vector in which,
- * for every vertex index, we
- * mark whether a locally owned
- * cell is adjacent.
- */
- template <int dim, int spacedim>
- std::vector<bool>
- mark_locally_active_vertices_on_level (const parallel::distributed::Triangulation<dim,spacedim> &triangulation, unsigned int level)
- {
- std::vector<bool> locally_active_vertices (triangulation.n_vertices(),
- false);
-
- if (level<triangulation.n_levels())
- for (typename dealii::Triangulation<dim,spacedim>::cell_iterator
- cell = triangulation.begin(level);
- cell != triangulation.end(level); ++cell)
- if (cell->level_subdomain_id() == triangulation.locally_owned_subdomain())
- for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
- locally_active_vertices[cell->vertex_index(v)] = true;
-
- return locally_active_vertices;
- }
template <int spacedim>
}
-
//sending
std::vector<std::vector<char> > sendbuffers (needs_to_get_cells.size());
std::vector<std::vector<char> >::iterator buffer = sendbuffers.begin();
}
-
//sending
std::vector<std::vector<char> > sendbuffers (needs_to_get_cells.size());
std::vector<std::vector<char> >::iterator buffer = sendbuffers.begin();
}
+
template <int dim, int spacedim>
void
ParallelDistributed<dim, spacedim>::
cell->set_user_flag();
}
- //mark the vertices we are interested
- //in, i.e. belonging to own and marked cells
+ //mark the vertices we are interested in, i.e. belonging to own
+ //and marked cells
const std::vector<bool> locally_active_vertices
- = mark_locally_active_vertices_on_level (*tr, level);
+ = tr->mark_locally_active_vertices_on_level (level);
- // add each ghostcells'
- // subdomain to the vertex and
- // keep track of interesting
- // neighbors
+ // add each ghostcells' subdomain to the vertex and keep track of
+ // interesting neighbors
std::map<unsigned int, std::set<dealii::types::subdomain_id> >
vertices_with_ghost_neighbors;
- if (level < tr->n_levels())
- for (typename DoFHandler<dim,spacedim>::level_cell_iterator
- cell = dof_handler.begin(level);
- cell != dof_handler.end(level); ++cell)
- if (cell->level_subdomain_id() != dealii::numbers::artificial_subdomain_id
- && cell->level_subdomain_id() != tr->locally_owned_subdomain())
- for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
- if (locally_active_vertices[cell->vertex_index(v)])
- vertices_with_ghost_neighbors[cell->vertex_index(v)]
- .insert (cell->level_subdomain_id());
-
-
- /* Send and receive cells. After this,
- only the local cells are marked,
- that received new data. This has to
- be communicated in a second
- communication step. */
+ tr->fill_level_vertices_with_ghost_neighbors(level,
+ vertices_with_ghost_neighbors);
+
+ // Send and receive cells. After this, only the local cells are
+ // marked, that received new data. This has to be communicated in
+ // a second communication step.
communicate_mg_dof_indices_on_marked_cells( dof_handler,
vertices_with_ghost_neighbors,
projects / dealii.git / commitdiff