From 31d65a4acd6260c6483a264865cc02d088c09672 Mon Sep 17 00:00:00 2001 From: Mathias Anselmann Date: Tue, 3 Nov 2020 12:51:59 +0100 Subject: [PATCH] Adjusted mg_constrained_dofs for usage of own IndexSet, so that the user can add additional constrained DoFs to the MGConstrainedDoFs class. --- .../deal.II/multigrid/mg_constrained_dofs.h | 100 ++++++++++++++++++ 1 file changed, 100 insertions(+) diff --git a/include/deal.II/multigrid/mg_constrained_dofs.h b/include/deal.II/multigrid/mg_constrained_dofs.h index 53e9d88dc7..9fc24d05b5 100644 --- a/include/deal.II/multigrid/mg_constrained_dofs.h +++ b/include/deal.II/multigrid/mg_constrained_dofs.h @@ -18,6 +18,7 @@ #include +#include #include #include @@ -69,6 +70,17 @@ public: void initialize(const DoFHandler &dof); + /** + * Same as the function above, but here the internal AffineConstraints objects + * on each level are initialized with the argument level_relevant_dofs. + * This should define a superset of locally relevant DoFs and allows the user + * to add additional indices to the set of constrained DoFs. + */ + template + void + initialize(const DoFHandler &dof, + const MGLevelObject & level_relevant_dofs); + /** * Fill the internal data structures with information * about Dirichlet boundary dofs. @@ -309,6 +321,94 @@ MGConstrainedDoFs::initialize(const DoFHandler &dof) } +template +inline void +MGConstrainedDoFs::initialize( + const DoFHandler &dof, + const MGLevelObject & level_relevant_dofs) +{ + boundary_indices.clear(); + refinement_edge_indices.clear(); + level_constraints.clear(); + user_constraints.clear(); + + const unsigned int nlevels = dof.get_triangulation().n_global_levels(); + + // At this point level_constraint and refinement_edge_indices are empty. + refinement_edge_indices.resize(nlevels); + level_constraints.resize(nlevels); + user_constraints.resize(nlevels); + for (unsigned int l = level_relevant_dofs.min_level(); + l <= level_relevant_dofs.max_level(); + ++l) + { + level_constraints[l].reinit(level_relevant_dofs[l]); + + // Loop through relevant cells and faces finding those which are periodic + // neighbors. + typename DoFHandler::cell_iterator cell = dof.begin(l), + endc = dof.end(l); + for (; cell != endc; ++cell) + if (cell->level_subdomain_id() != numbers::artificial_subdomain_id) + { + for (auto f : cell->face_indices()) + if (cell->has_periodic_neighbor(f) && + cell->periodic_neighbor(f)->level() == cell->level()) + { + if (cell->is_locally_owned_on_level()) + { + Assert( + cell->periodic_neighbor(f)->level_subdomain_id() != + numbers::artificial_subdomain_id, + ExcMessage( + "Periodic neighbor of a locally owned cell must either be owned or ghost.")); + } + // Cell is a level-ghost and its neighbor is a + // level-artificial cell nothing to do here + else if (cell->periodic_neighbor(f)->level_subdomain_id() == + numbers::artificial_subdomain_id) + { + Assert(cell->is_locally_owned_on_level() == false, + ExcInternalError()); + continue; + } + + const unsigned int dofs_per_face = + dof.get_fe(0).n_dofs_per_face(f); + std::vector dofs_1(dofs_per_face); + std::vector dofs_2(dofs_per_face); + + cell->periodic_neighbor(f) + ->face(cell->periodic_neighbor_face_no(f)) + ->get_mg_dof_indices(l, dofs_1, 0); + cell->face(f)->get_mg_dof_indices(l, dofs_2, 0); + // Store periodicity information in the level + // AffineConstraints object. Skip DoFs for which we've + // previously entered periodicity constraints already; this + // can happen, for example, for a vertex dof at a periodic + // boundary that we visit from more than one cell + for (unsigned int i = 0; i < dofs_per_face; ++i) + if (level_constraints[l].can_store_line(dofs_2[i]) && + level_constraints[l].can_store_line(dofs_1[i]) && + !level_constraints[l].is_constrained(dofs_2[i]) && + !level_constraints[l].is_constrained(dofs_1[i])) + { + level_constraints[l].add_line(dofs_2[i]); + level_constraints[l].add_entry(dofs_2[i], + dofs_1[i], + 1.); + } + } + } + level_constraints[l].close(); + + // Initialize with empty IndexSet of correct size + refinement_edge_indices[l] = IndexSet(dof.n_dofs(l)); + } + + MGTools::extract_inner_interface_dofs(dof, refinement_edge_indices); +} + template inline void -- 2.39.5