From 79ea0c58d3831dcc0d130e19722490bf5281f15c Mon Sep 17 00:00:00 2001 From: wolf Date: Thu, 10 Oct 2002 17:50:52 +0000 Subject: [PATCH] Implement a different scheme for generation of neighborship information. This should remove the previously existing bottleneck in generation of patches -- setting neighborship information could take up two thirds or more of the total time, and it was an algorithm at least O(N log N) in complexity, compared to the O(N) of the rest. The new algorithm should bot take a noticable time compared to all the other things that are happening. git-svn-id: https://svn.dealii.org/trunk@6624 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/numerics/data_out.h | 84 +++---- deal.II/deal.II/source/numerics/data_out.cc | 230 ++++++++++---------- 2 files changed, 161 insertions(+), 153 deletions(-) diff --git a/deal.II/deal.II/include/numerics/data_out.h b/deal.II/deal.II/include/numerics/data_out.h index f9a8cda4d2..ce388d7a50 100644 --- a/deal.II/deal.II/include/numerics/data_out.h +++ b/deal.II/deal.II/include/numerics/data_out.h @@ -406,13 +406,16 @@ class DataOut_DoFData : public DataOutInterface * Declare an entry in the list of * data elements. */ - struct DataEntry { + struct DataEntry + { /** - * Constructor. If no arguments are - * given, an invalid object is - * constructed (we need a constructor - * with no explicit arguments for - * STL classes). + * Constructor. If no + * arguments are given, an + * invalid object is + * constructed (we need a + * constructor with no + * explicit arguments for STL + * classes). */ DataEntry (const Vector *data = 0, const std::vector &names = std::vector()); @@ -424,9 +427,9 @@ class DataOut_DoFData : public DataOutInterface const std::vector &names = std::vector()); /** - * Determine an estimate for the - * memory consumption (in bytes) - * of this object. + * Determine an estimate for + * the memory consumption (in + * bytes) of this object. */ unsigned int memory_consumption () const; @@ -437,21 +440,22 @@ class DataOut_DoFData : public DataOutInterface * pointed to remains with * the caller of this class. */ - const Vector *single_data; - - /** - * Pointer to a block vector of - * data. Either this pointer - * or the previous is used, - * depending on the stored - * vector type. - */ - const BlockVector* block_data; - - /** - * True if stored vector is a block vector. - */ - bool has_block; + const Vector *single_data; + + /** + * Pointer to a block vector of + * data. Either this pointer + * or the previous is used, + * depending on the stored + * vector type. + */ + const BlockVector* block_data; + + /** + * True if stored vector is a + * block vector. + */ + bool has_block; /** * Names of the components of this @@ -638,23 +642,24 @@ class DataOut : public DataOut_DoFData << "The number of subdivisions per patch, " << arg1 << ", is not valid."); - /** - * Define a map from cell - * iterators to patch iterators - * for finding neighbors. - */ - typedef std::map::cell_iterator, - typename std::vector< ::DataOutBase::Patch >::iterator> - PatchMap; - private: /** * All data needed in one thread * is gathered in the struct - * Data. - * The data is handled globally - * to avoid allocation of memory - * in the threads. + * Data. The data is handled + * globally to avoid allocation + * of memory in the threads. + * + * The @p{index_to_patch_map} is + * an array that stores for index + * @p{[i][j]} the number of the + * patch that associated with the + * cell with index @p{j} on level + * @p{i}. This information is set + * up prior to generation of the + * patches, and is needed to + * generate neighborship + * information. */ struct Data { @@ -665,9 +670,8 @@ class DataOut : public DataOut_DoFData unsigned int n_subdivisions; std::vector patch_values; std::vector > patch_values_system; - PatchMap* patch_map; - Data () - {} + + std::vector > *cell_to_patch_index_map; }; /** * Builds every @p{n_threads}'s diff --git a/deal.II/deal.II/source/numerics/data_out.cc b/deal.II/deal.II/source/numerics/data_out.cc index 91b3f161b6..95a0ba9169 100644 --- a/deal.II/deal.II/source/numerics/data_out.cc +++ b/deal.II/deal.II/source/numerics/data_out.cc @@ -369,10 +369,6 @@ void DataOut::build_some_patches (Data data) { Assert (patch != this->patches.end(), ExcInternalError()); - // Insert cell-patch pair into - // map for neighbors - data.patch_map->insert(typename PatchMap::value_type(cell, patch)); - for (unsigned int vertex=0; vertex::vertices_per_cell; ++vertex) patch->vertices[vertex] = cell->vertex(vertex); @@ -430,7 +426,78 @@ void DataOut::build_some_patches (Data data) } }; }; - // next cell (patch) in this thread + + // now fill the neighbors fields + for (unsigned int f=0; f::faces_per_cell; ++f) + { + // let's look up whether + // the neighbor behind that + // face is noted in the + // table of cells which we + // treat. this can only + // happen if the neighbor + // exists, and is on the + // same level as this cell, + // but it may also happen + // that the neighbor is not + // a member of the range of + // cells over which we loop + if (cell->at_boundary(f) + || + (cell->neighbor(f)->level() != cell->level())) + continue; + + const typename DoFHandler::cell_iterator + neighbor = cell->neighbor(f); + if ((*data.cell_to_patch_index_map)[neighbor->level()][neighbor->index()] + == + ::DataOutBase::Patch::no_neighbor) + continue; + + // now, there is a + // neighbor, so get its + // patch number and set it + // for the neighbor index + const unsigned int neighbor_patch_index + = this->patches[(*data.cell_to_patch_index_map) + [neighbor->level()][neighbor->index()]].patch_index; + +//TODO:[GK] Shouldn't we use the deal.II (i.e. the unnatural) numbering of the neighbors here as well rather than some new numbering scheme? + + switch (dim) + { + case 1: + patch->neighbors[f] = neighbor_patch_index; + break; + + case 2: + switch (f) + { + case 0: patch->neighbors[2] = neighbor_patch_index; break; + case 1: patch->neighbors[1] = neighbor_patch_index; break; + case 2: patch->neighbors[3] = neighbor_patch_index; break; + case 3: patch->neighbors[0] = neighbor_patch_index; break; + } + break; + case 3: + switch (f) + { + case 0: patch->neighbors[2] = neighbor_patch_index; break; + case 1: patch->neighbors[3] = neighbor_patch_index; break; + case 2: patch->neighbors[4] = neighbor_patch_index; break; + case 3: patch->neighbors[1] = neighbor_patch_index; break; + case 4: patch->neighbors[5] = neighbor_patch_index; break; + case 5: patch->neighbors[0] = neighbor_patch_index; break; + } + break; + + default: + Assert(false, ExcNotImplemented()); + } + }; + + // next cell (patch) in this + // thread for (unsigned int i=0; (idofs->end()); ++i) { @@ -483,20 +550,55 @@ void DataOut::build_patches (const unsigned int n_subdivisions, }; // first count the cells we want to - // create patches of and make sure - // there is enough memory for that + // create patches of. also fill the + // // object that maps the cell + // indices to the patch numbers, as + // this will be needed for + // generation of neighborship + // information + std::vector > cell_to_patch_index_map; + cell_to_patch_index_map.resize (this->dofs->get_tria().n_levels()); + for (unsigned int l=0; ldofs->get_tria().n_levels(); ++l) + cell_to_patch_index_map[l].resize (this->dofs->get_tria().n_cells(l), + ::DataOutBase::Patch::no_neighbor); + unsigned int n_patches = 0; for (typename DoFHandler::cell_iterator cell=first_cell(); cell != this->dofs->end(); cell = next_cell(cell)) - ++n_patches; + { + Assert (static_cast(cell->level()) < + cell_to_patch_index_map.size(), + ExcInternalError()); + Assert (static_cast(cell->index()) < + cell_to_patch_index_map[cell->level()].size(), + ExcInternalError()); + + cell_to_patch_index_map[cell->level()][cell->index()] = n_patches; + ++n_patches; + }; + + // create the patches with default + // values. allocate as many patches + // as are needed, as this reduces + // expensive copying when push_back + // or similar operations are used + // which would regularly overflow + // the allocated amount of memory + // + // then number the patches + // consecutively + ::DataOutBase::Patch default_patch; + default_patch.n_subdivisions = n_subdivisions; + default_patch.data.reinit (n_datasets, n_q_points); + this->patches.insert (this->patches.end(), n_patches, default_patch); + for (unsigned int i=0; ipatches.size(); ++i) + this->patches[i].patch_index = i; - PatchMap patch_map; - - std::vector thread_data(n_threads); - // init data for the threads + // init data for the threads + std::vector thread_data(n_threads); for (unsigned int i=0;i::build_patches (const unsigned int n_subdivisions, thread_data[i].n_components = n_components; thread_data[i].n_datasets = n_datasets; thread_data[i].n_subdivisions = n_subdivisions; - thread_data[i].patch_map = &patch_map; thread_data[i].patch_values.resize (n_q_points); thread_data[i].patch_values_system.resize (n_q_points); + + thread_data[i].cell_to_patch_index_map = &cell_to_patch_index_map; for (unsigned int k=0; k default_patch; - default_patch.n_subdivisions = n_subdivisions; - default_patch.data.reinit (n_datasets, n_q_points); - this->patches.insert (this->patches.end(), n_patches, default_patch); - -#ifdef DEAL_II_USE_MT - Threads::ThreadManager thread_manager; - for (unsigned int l=0;l::build_some_patches) .collect_args (this, thread_data[l])); thread_manager.wait(); - - // just one thread -#else - build_some_patches(thread_data[0]); -#endif -//TODO: New code for neighbors: Parallelize? - - // First, number patches - // consecutively. - unsigned int patch_no = 0; - for (typename std::vector< ::DataOutBase::Patch >::iterator - patch = this->patches.begin(); patch != this->patches.end(); ++patch) - patch->me = patch_no++; - - // Traverse the map of cells - // created in build_some_patches - // and enter numbers of neighboring - // patches on the same level. - for (typename PatchMap::iterator map_entry = patch_map.begin(); - map_entry != patch_map.end(); - ++map_entry) - { - const typename DoFHandler::cell_iterator - cell = map_entry->first; - const typename std::vector< ::DataOutBase::Patch >::iterator - patch = map_entry->second; - - // loop over all faces and see - // whether there's a cell - // behind that for which we - // need to set neighbor - // information - for (unsigned int i=0;i::faces_per_cell;++i) - { - // we can't do anything if - // - there is no neighbor - // - the neighbor is coarser - // - or the cell behind the - // face is refined, i.e. is - // not inserted into the map - if (cell->at_boundary(i) || - (cell->level() != cell->neighbor(i)->level()) || - (patch_map.find(cell->neighbor(i)) == patch_map.end())) - continue; - - // if there is a neighbor, - // get its patch... - const typename - std::vector< ::DataOutBase::Patch >::iterator - neighbor_patch = patch_map.find(cell->neighbor(i))->second; - // ...and set its neighbor - // pointer set -//TODO:[GK] Shouldn't we use the deal.II (i.e. the unnatural) numbering of the neighbors here as well rather than some new numbering scheme? - switch (dim) - { - case 1: - patch->neighbors[i] = neighbor_patch->me; - break; - - case 2: - switch (i) - { - case 0: patch->neighbors[2] = neighbor_patch->me; break; - case 1: patch->neighbors[1] = neighbor_patch->me; break; - case 2: patch->neighbors[3] = neighbor_patch->me; break; - case 3: patch->neighbors[0] = neighbor_patch->me; break; - } - break; - case 3: - switch (i) - { - case 0: patch->neighbors[2] = neighbor_patch->me; break; - case 1: patch->neighbors[3] = neighbor_patch->me; break; - case 2: patch->neighbors[4] = neighbor_patch->me; break; - case 3: patch->neighbors[1] = neighbor_patch->me; break; - case 4: patch->neighbors[5] = neighbor_patch->me; break; - case 5: patch->neighbors[0] = neighbor_patch->me; break; - } - break; - - default: - Assert(false, ExcNotImplemented()); - } - } - } }; -- 2.39.5