From: Wolfgang Bangerth Date: Fri, 23 Jun 2017 13:09:53 +0000 (-0600) Subject: Let the DoFHandler policy classes' distribute_dofs() return the NumberCache. X-Git-Tag: v9.0.0-rc1~1475^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a537c62ab8fd9c61d6022bcb0a6dc9e7467f6cc6;p=dealii.git Let the DoFHandler policy classes' distribute_dofs() return the NumberCache. In times of move operators, there is really no reason why these functions should return the NumberCache object by reference, rather than by value. We make this slightly more efficient by providing (defaulted) move operations to NumberCache. However, this change uncovered a subtle problem that sooner or later would certainly have tripped someone up: If any of the policy classes call any function that relies on the correctness of the number cache of the DoFHandler object being worked on, then this is likely going to fail. Indeed, the implementation of distribute_dofs() for shared triangulations has this problem by calling some of the DoFTools functions that query properties of the DoFHandler that aren't all completely consistent yet internally. This happened to work -- more or less by design -- because the NumberCache object that the DoFHandler passes to distribute_dofs() by reference is the one that belongs to the DoFHandler, and distribute_dofs() makes sure that it sets the one member that is necessary to the correct value before calling the DoFTools function. This can only properly be described as 'brittle': small changes to the code will make this fail, and it is not obvious why this works or doesn't. The only proper solution, obviously, is to not call external functions from inside distribute_dofs() until all data structures have been set up completely and are internally consistent. I work around this by duplicating a small amount of code (specifically, the DoFTools::get_subdomain_association() and DoFRenumbering::subdomain_wise() functions), but this effort is made more palatable because these functions are in fact more general than we need here, and we can throw away a good amount of code that deals with either sequential or completely distributed DoFHandlers; we also know that a couple of assertions in these functions must be satisfied here, so we can eliminate them as well. The consequence of this code is slightly more code that is more robust and contains fewer surprises. --- diff --git a/include/deal.II/dofs/dof_handler_policy.h b/include/deal.II/dofs/dof_handler_policy.h index f9c31bbb6c..4e725dfb66 100644 --- a/include/deal.II/dofs/dof_handler_policy.h +++ b/include/deal.II/dofs/dof_handler_policy.h @@ -70,8 +70,8 @@ namespace internal * argument. */ virtual - void - distribute_dofs (NumberCache &number_cache) const = 0; + NumberCache + distribute_dofs () const = 0; /** * Distribute the multigrid dofs on each level of the DoFHandler @@ -113,8 +113,8 @@ namespace internal // documentation is inherited virtual - void - distribute_dofs (NumberCache &number_cache) const; + NumberCache + distribute_dofs () const; // documentation is inherited virtual @@ -158,8 +158,8 @@ namespace internal * number_cache.locally_owned_dofs are updated consistently. */ virtual - void - distribute_dofs (NumberCache &number_cache) const; + NumberCache + distribute_dofs () const; /** * This function is not yet implemented. @@ -201,8 +201,8 @@ namespace internal // documentation is inherited virtual - void - distribute_dofs (NumberCache &number_cache) const; + NumberCache + distribute_dofs () const; // documentation is inherited virtual diff --git a/include/deal.II/dofs/number_cache.h b/include/deal.II/dofs/number_cache.h index c101dfe949..07615788e0 100644 --- a/include/deal.II/dofs/number_cache.h +++ b/include/deal.II/dofs/number_cache.h @@ -39,6 +39,18 @@ namespace internal */ NumberCache (); + /** + * Move constructor. Simply move all members of the referenced + * object to the current object. + */ + NumberCache (NumberCache &&) = default; + + /** + * Move assignment operator. Simply move all members of the referenced + * object to the current object. + */ + NumberCache &operator= (NumberCache &&) = default; + /** * Determine an estimate for the memory consumption (in bytes) of this * object. diff --git a/source/dofs/dof_handler.cc b/source/dofs/dof_handler.cc index d90bd71e20..08f9011c18 100644 --- a/source/dofs/dof_handler.cc +++ b/source/dofs/dof_handler.cc @@ -1098,7 +1098,7 @@ void DoFHandler::distribute_dofs (const FiniteElementdistribute_dofs (number_cache); + number_cache = policy->distribute_dofs (); // initialize the block info object // only if this is a sequential diff --git a/source/dofs/dof_handler_policy.cc b/source/dofs/dof_handler_policy.cc index bf23ce3dd7..681056c601 100644 --- a/source/dofs/dof_handler_policy.cc +++ b/source/dofs/dof_handler_policy.cc @@ -922,9 +922,9 @@ namespace internal template - void + NumberCache Sequential:: - distribute_dofs (NumberCache &number_cache_current) const + distribute_dofs () const { const types::global_dof_index n_dofs = Implementation::distribute_dofs (0, @@ -950,7 +950,8 @@ namespace internal number_cache.locally_owned_dofs_per_processor = std::vector (1, number_cache.locally_owned_dofs); - number_cache_current = number_cache; + + return number_cache; } @@ -1011,7 +1012,7 @@ namespace internal number_cache.locally_owned_dofs_per_processor = std::vector (1, number_cache.locally_owned_dofs); - number_cache_current = number_cache; + number_cache_current = std::move(number_cache); } @@ -1027,55 +1028,243 @@ namespace internal + namespace + { + /** + * This function is a variation of DoFTools::get_dof_subdomain_association() + * with the exception that it is (i) specialized for + * parallel::shared::Triangulation objects, and (ii) does not assume that the + * internal number cache of the DoFHandler has already been set up. In can, + * consequently, be called at a point when we are still distributing degrees + * of freedom. + */ + template + std::vector + get_dof_subdomain_association (const DoFHandlerType &dof_handler, + const types::global_dof_index n_dofs, + const unsigned int n_procs) + { + std::vector subdomain_association (n_dofs, + numbers::invalid_subdomain_id); + std::vector local_dof_indices; + local_dof_indices.reserve (DoFTools::max_dofs_per_cell(dof_handler)); + + // pseudo-randomly assign variables which lie on the interface between + // subdomains to each of the two or more + bool coin_flip = true; + + // loop over all cells and record which subdomain a DoF belongs to. + // toss a coin in case it is on an interface + typename DoFHandlerType::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + // get the owner of the cell; note that we have made sure above that + // all cells are either locally owned or ghosts (not artificial), so + // this call will always yield the true owner + const types::subdomain_id subdomain_id = cell->subdomain_id(); + const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; + local_dof_indices.resize (dofs_per_cell); + cell->get_dof_indices (local_dof_indices); + + // set subdomain ids. if dofs already have their values set then + // they must be on partition interfaces. in that case randomly + // assign them to either the previous association or the current + // one, where we take "random" to be "once this way once that way" + for (unsigned int i=0; i - void + NumberCache ParallelShared:: - distribute_dofs (NumberCache &number_cache) const + distribute_dofs () const { - // If the underlying shared::Tria allows artificial cells, we need to do - // some tricks here to make Sequential algorithms play nicely. - // Namely, we first restore the original partition (without artificial cells) - // and then turn artificial cells on at the end of this function. const parallel::shared::Triangulation *tr = (dynamic_cast*> (&this->dof_handler->get_triangulation())); Assert(tr != nullptr, ExcInternalError()); - typename parallel::shared::Triangulation::active_cell_iterator - cell = this->dof_handler->get_triangulation().begin_active(), - endc = this->dof_handler->get_triangulation().end(); - std::vector current_subdomain_ids(tr->n_active_cells()); - const std::vector &true_subdomain_ids = tr->get_true_subdomain_ids_of_cells(); + + const unsigned int n_procs = Utilities::MPI::n_mpi_processes(tr->get_communicator()); + + // If the underlying shared::Tria allows artificial cells, + // then save the current set of subdomain ids, and set + // subdomain ids to the "true" owner of each cell. we later + // restore these flags + std::vector saved_subdomain_ids; if (tr->with_artificial_cells()) - for (unsigned int index=0; cell != endc; cell++, index++) + { + saved_subdomain_ids.resize (tr->n_active_cells()); + + typename parallel::shared::Triangulation::active_cell_iterator + cell = this->dof_handler->get_triangulation().begin_active(), + endc = this->dof_handler->get_triangulation().end(); + + const std::vector &true_subdomain_ids + = tr->get_true_subdomain_ids_of_cells(); + + for (unsigned int index=0; cell != endc; ++cell, ++index) + { + saved_subdomain_ids[index] = cell->subdomain_id(); + cell->set_subdomain_id(true_subdomain_ids[index]); + } + } + + // first let the sequential algorithm do its magic. it is going to + // enumerate DoFs on all cells, regardless of owner + NumberCache number_cache + = this->Sequential::distribute_dofs (); + + // then re-enumerate them based on their subdomain association. + // for this, we first have to identify for each current DoF + // index which subdomain they belong to. ideally, we would + // like to call DoFRenumbering::subdomain_wise(), but + // because the NumberCache of the current DoFHandler is not + // fully set up yet, we can't quite do that. also, that + // function has to deal with other kinds of triangulations as + // well, whereas we here know what kind of triangulation + // we have and can simplify the code accordingly + const unsigned int n_dofs = number_cache.n_global_dofs; + std::vector new_dof_indices (n_dofs, + numbers::invalid_dof_index); + { + // first get the association of each dof with a subdomain and + // determine the total number of subdomain ids used + const std::vector subdomain_association + = get_dof_subdomain_association (*this->dof_handler, n_dofs, n_procs); + + // then renumber the subdomains by first looking at those belonging + // to subdomain 0, then those of subdomain 1, etc. note that the + // algorithm is stable, i.e. if two dofs i,j have iSequential::renumber_dofs (new_dof_indices, number_cache); + + + // update the number cache. for this, we first have to find the subdomain + // association for each DoF again following renumbering, from which we + // can then compute the IndexSets of locally owned DoFs for all processors. + // all other fields then follow from this + // + // given the way we enumerate degrees of freedom, the locally owned + // ranges must all be contiguous and consecutive. this makes filling + // the IndexSets cheap. an assertion at the top verifies that this + // assumption is true + const std::vector subdomain_association + = get_dof_subdomain_association (*this->dof_handler, n_dofs, n_procs); + + for (unsigned int i=1; i= subdomain_association[i-1], + ExcInternalError()); + + { + number_cache.locally_owned_dofs_per_processor.clear(); + number_cache.locally_owned_dofs_per_processor.resize (n_procs, + IndexSet(n_dofs)); + + // we know that the set of subdomain indices is contiguous from + // the assertion above; find the start and end index for each + // processor, taking into account that sometimes a processor + // may not in fact have any DoFs at all. we do the latter + // by just identifying contiguous ranges of subdomain_ids + // and filling IndexSets for those subdomains; subdomains + // that don't appear will lead to IndexSets that are simply + // never touched and remain empty as initialized above. + unsigned int start_index = 0; + unsigned int end_index = 0; + while (start_index < n_dofs) { - current_subdomain_ids[index] = cell->subdomain_id(); - cell->set_subdomain_id(true_subdomain_ids[index]); + while ((end_index) < n_dofs && + (subdomain_association[end_index] == subdomain_association[start_index])) + ++end_index; + + // we've now identified a range of same indices. set that + // range in the corresponding IndexSet + if (end_index > start_index) + { + const unsigned int subdomain_owner = subdomain_association[start_index]; + number_cache.locally_owned_dofs_per_processor[subdomain_owner] + .add_range (start_index, end_index); + } + + // then move on to thinking about the next range + start_index = end_index; } + } - // let the sequential algorithm do its magic, then sort DoF indices - // by subdomain - this->Sequential::distribute_dofs (number_cache); - DoFRenumbering::subdomain_wise (*this->dof_handler); + // from this, calculate all of the other number cache fields + number_cache.locally_owned_dofs + = number_cache.locally_owned_dofs_per_processor[this->dof_handler->get_triangulation().locally_owned_subdomain()]; + number_cache.n_locally_owned_dofs_per_processor + .resize (number_cache.locally_owned_dofs_per_processor.size()); + for (unsigned int i=0; idof_handler->get_triangulation().locally_owned_subdomain()]; - // dofrenumbering will reset subdomains, this is ugly but we need to do it again: - cell = tr->begin_active(); - if (tr->with_artificial_cells()) - for (unsigned int index=0; cell != endc; cell++, index++) - cell->set_subdomain_id(true_subdomain_ids[index]); + number_cache.n_global_dofs = n_dofs; - number_cache.locally_owned_dofs_per_processor = DoFTools::locally_owned_dofs_per_subdomain (*this->dof_handler); - number_cache.locally_owned_dofs = number_cache.locally_owned_dofs_per_processor[this->dof_handler->get_triangulation().locally_owned_subdomain()]; - number_cache.n_locally_owned_dofs_per_processor.resize (number_cache.locally_owned_dofs_per_processor.size()); - for (unsigned int i = 0; i < number_cache.n_locally_owned_dofs_per_processor.size(); i++) - number_cache.n_locally_owned_dofs_per_processor[i] = number_cache.locally_owned_dofs_per_processor[i].n_elements(); - number_cache.n_locally_owned_dofs = number_cache.n_locally_owned_dofs_per_processor[this->dof_handler->get_triangulation().locally_owned_subdomain()]; - // restore current subdomain ids - cell = tr->begin_active(); + // finally, restore current subdomain ids if (tr->with_artificial_cells()) - for (unsigned int index=0; cell != endc; cell++, index++) - cell->set_subdomain_id(current_subdomain_ids[index]); + { + typename parallel::shared::Triangulation::active_cell_iterator + cell = this->dof_handler->get_triangulation().begin_active(), + endc = this->dof_handler->get_triangulation().end(); + + for (unsigned int index=0; cell != endc; ++cell, ++index) + cell->set_subdomain_id(saved_subdomain_ids[index]); + } + + return number_cache; } @@ -2088,9 +2277,9 @@ namespace internal template - void + NumberCache ParallelDistributed:: - distribute_dofs (NumberCache &number_cache_current) const + distribute_dofs () const { NumberCache number_cache; @@ -2289,7 +2478,7 @@ namespace internal #endif // DEBUG #endif // DEAL_II_WITH_P4EST - number_cache_current = number_cache; + return number_cache; } @@ -2777,7 +2966,7 @@ namespace internal } #endif - number_cache_current = number_cache; + number_cache_current = std::move(number_cache); } } } diff --git a/source/dofs/dof_tools.cc b/source/dofs/dof_tools.cc index b376be5f63..6c03ad4e01 100644 --- a/source/dofs/dof_tools.cc +++ b/source/dofs/dof_tools.cc @@ -1422,9 +1422,11 @@ namespace DoFTools ExcDimensionMismatch(subdomain_association.size(), dof_handler.n_dofs())); - Assert(dof_handler.n_dofs() > 0, - ExcMessage("Number of DoF is not positive. " - "This could happen when the function is called before NumberCache is written.")); + // catch an error that happened in some versions of the shared tria + // distribute_dofs() function where we were trying to call this + // function at a point in time when not all internal DoFHandler + // structures were quite set up yet. + Assert(dof_handler.n_dofs() > 0, ExcInternalError()); // In case this function is executed with parallel::shared::Triangulation // with possibly artifical cells, we need to take "true" subdomain IDs (i.e. without