]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement a different scheme for generation of neighborship information. This should...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 10 Oct 2002 17:50:52 +0000 (17:50 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 10 Oct 2002 17:50:52 +0000 (17:50 +0000)
git-svn-id: https://svn.dealii.org/trunk@6624 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/data_out.h
deal.II/deal.II/source/numerics/data_out.cc

index f9a8cda4d2fa5b846fc88c2849222d2f5161740f..ce388d7a50b6fe1e1be2e8a18e2cb735a1ca2243 100644 (file)
@@ -406,13 +406,16 @@ class DataOut_DoFData : public DataOutInterface<patch_dim,patch_space_dim>
                                      * 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<double>           *data = 0,
                   const std::vector<std::string> &names = std::vector<std::string>());
@@ -424,9 +427,9 @@ class DataOut_DoFData : public DataOutInterface<patch_dim,patch_space_dim>
                   const std::vector<std::string> &names = std::vector<std::string>());
 
                                         /**
-                                         * 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<patch_dim,patch_space_dim>
                                          * pointed to remains with
                                          * the caller of this class.
                                          */
-      const Vector<double> *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<double>* block_data;
-
-                                      /**
-                                       * True if stored vector is a block vector.
-                                       */
-      bool has_block;
+        const Vector<double> *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<double>* 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<dim,dim>
                    << "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<typename DoFHandler<dim>::cell_iterator,
-      typename std::vector< ::DataOutBase::Patch<dim> >::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<dim,dim>
        unsigned int n_subdivisions;
        std::vector<double>          patch_values;
        std::vector<Vector<double> > patch_values_system;
-       PatchMap* patch_map;
-       Data ()
-         {}
+
+        std::vector<std::vector<unsigned int> > *cell_to_patch_index_map;
     };
                                     /**
                                      * Builds every @p{n_threads}'s
index 91b3f161b64ec87107b8cac38222eb16f7bd4607..95a0ba9169f7da32b31b6807a45351f3d41a83c0 100644 (file)
@@ -369,10 +369,6 @@ void DataOut<dim>::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<GeometryInfo<dim>::vertices_per_cell; ++vertex)
        patch->vertices[vertex] = cell->vertex(vertex);
       
@@ -430,7 +426,78 @@ void DataOut<dim>::build_some_patches (Data data)
                } 
            };
        };
-                                      // next cell (patch) in this thread
+
+                                       // now fill the neighbors fields
+      for (unsigned int f=0; f<GeometryInfo<dim>::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<dim>::cell_iterator
+            neighbor = cell->neighbor(f);
+          if ((*data.cell_to_patch_index_map)[neighbor->level()][neighbor->index()]
+              ==
+              ::DataOutBase::Patch<dim>::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;
           (i<data.n_threads)&&(cell != this->dofs->end()); ++i)
        {
@@ -483,20 +550,55 @@ void DataOut<dim>::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<std::vector<unsigned int> > cell_to_patch_index_map;
+  cell_to_patch_index_map.resize (this->dofs->get_tria().n_levels());
+  for (unsigned int l=0; l<this->dofs->get_tria().n_levels(); ++l)
+    cell_to_patch_index_map[l].resize (this->dofs->get_tria().n_cells(l),
+                                       ::DataOutBase::Patch<dim>::no_neighbor);
+                                  
   unsigned int n_patches = 0;
   for (typename DoFHandler<dim>::cell_iterator cell=first_cell();
        cell != this->dofs->end();
        cell = next_cell(cell))
-    ++n_patches;
+    {
+      Assert (static_cast<unsigned int>(cell->level()) <
+              cell_to_patch_index_map.size(),
+              ExcInternalError());
+      Assert (static_cast<unsigned int>(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<dim>  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; i<this->patches.size(); ++i)
+    this->patches[i].patch_index = i;
   
-  PatchMap patch_map;
-    
-  std::vector<Data> thread_data(n_threads);
 
-                                  // init data for the threads
+                                  // init data for the threads    
+  std::vector<Data> thread_data(n_threads);
   for (unsigned int i=0;i<n_threads;++i)
     {
       thread_data[i].n_threads      = n_threads;
@@ -504,119 +606,21 @@ void DataOut<dim>::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<n_q_points; ++k)
        thread_data[i].patch_values_system[k].reinit(n_components);
     }
 
-                                  // 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
-  ::DataOutBase::Patch<dim>  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<n_threads;++l)
+  for (unsigned int l=0; l<n_threads; ++l)
     Threads::spawn (thread_manager,
                    Threads::encapsulate (&DataOut<dim>::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<dim> >::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<dim>::cell_iterator
-       cell = map_entry->first;
-      const typename std::vector< ::DataOutBase::Patch<dim> >::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<GeometryInfo<dim>::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<dim> >::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());
-           }
-       }
-    }
 };
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.