neighbors for DX in 2D

git-svn-id: https://svn.dealii.org/trunk@5334 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/base/include/base/data_out_base.h b/deal.II/base/include/base/data_out_base.h
index 6dd7367714e6542748c1fc27bddacd8c1e4e6b8a..5099334abd355dbea48995681e5d0af1fa735a4d 100644 (file)
--- a/deal.II/base/include/base/data_out_base.h
+++ b/deal.II/base/include/base/data_out_base.h
@@ -479,6 +479,11 @@ class DataOutBase
                                          */
        unsigned int memory_consumption () const;
        
+                                        /**
+                                         * Value for no neighbor.
+                                         */
+       static const unsigned int no_neighbor = static_cast<unsigned int>(-1);
+       
                                         /**
                                          * Exception
                                          */
@@ -497,18 +502,23 @@ class DataOutBase
                                      */
     struct DXFlags 
     {
-      private:
                                         /**
-                                         * Dummy entry to suppress compiler
-                                         * warnings when copying an empty
-                                         * structure. Remove this member
-                                         * when adding the first flag to
-                                         * this structure (and remove the
-                                         * @p{private} as well).
+                                         * Write in multigrid format.
+                                         * This will be necessary to
+                                         * get streamlines right on
+                                         * locally refined grids.
                                          */
-       int dummy;
-
-      public:
+       bool write_multigrid;
+                                        /**
+                                         * Write neighbor information.
+                                         */
+       bool write_neighbors;
+                                        /**
+                                         * Constructor.
+                                         */
+       DXFlags (const bool write_multigrid = false,
+                const bool write_neighbors = false);
+public:
                                         /**
                                          * Declare all flags with name
                                          * and type as offered by this
@@ -1373,6 +1383,7 @@ class DataOutBase
     static void
     write_gmv_reorder_data_vectors (const typename std::vector<Patch<dim,spacedim> > &patches,
                                    std::vector<std::vector<double> >       &data_vectors);
+
 };
 
 

diff --git a/deal.II/base/source/data_out_base.all_dimensions.cc b/deal.II/base/source/data_out_base.all_dimensions.cc
index 35be247473dc957df49fe132b8d55e55b1381253..b2c5f77b93f0d7e54a9323ef15fb86760907782e 100644 (file)
--- a/deal.II/base/source/data_out_base.all_dimensions.cc
+++ b/deal.II/base/source/data_out_base.all_dimensions.cc
@@ -32,14 +32,26 @@ DataOutBase::PovrayFlags::PovrayFlags (const bool smooth,
 {};
 
 
+DataOutBase::DXFlags::DXFlags (const bool write_multigrid,
+                              const bool write_neighbors) :
+               write_multigrid(write_multigrid),
+               write_neighbors(write_neighbors)
+{}
 
-void DataOutBase::DXFlags::declare_parameters (ParameterHandler &/*prm*/)
-{};
 
+void DataOutBase::DXFlags::declare_parameters (ParameterHandler &prm)
+{
+  prm.declare_entry ("Write multigrid", "true", Patterns::Bool());
+  prm.declare_entry ("Write neighbors", "true", Patterns::Bool());
+};
 
 
-void DataOutBase::DXFlags::parse_parameters (ParameterHandler &/*prm*/)
-{};
+
+void DataOutBase::DXFlags::parse_parameters (ParameterHandler &prm)
+{
+  write_multigrid = prm.get_bool ("Write multigrid");
+  write_neighbors = prm.get_bool ("Write neighbors");
+};
 
 
 

diff --git a/deal.II/base/source/data_out_base.cc b/deal.II/base/source/data_out_base.cc
index 38d70727c943b88c1d7a79fbcaa756764dd629b8..6fc1719f3091bc398d0b37c291a740faa3b650de 100644 (file)
--- a/deal.II/base/source/data_out_base.cc
+++ b/deal.II/base/source/data_out_base.cc
@@ -32,10 +32,15 @@ using namespace std;
 
 template <int dim, int spacedim>
 DataOutBase::Patch<dim,spacedim>::Patch () :
+               me(no_neighbor),
                n_subdivisions (1)
+
                                   // all the other data has a
                                   // constructor of its own
 {
+  for (unsigned int i=0;i<GeometryInfo<dim>::faces_per_cell;++i)
+    neighbors[i] = no_neighbor;
+  
   Assert (dim<=spacedim, ExcInvalidCombinationOfDimensions(dim,spacedim));
   Assert (spacedim<=3, ExcNotImplemented());
 };
@@ -469,7 +474,7 @@ void DataOutBase::write_ucd (const typename std::vector<Patch<dim,spacedim> > &p
 template <int dim, int spacedim>
 void DataOutBase::write_dx (const typename std::vector<Patch<dim,spacedim> > &patches,
                            const std::vector<std::string>          &data_names,
-                           const DXFlags                           &/*flags*/,
+                           const DXFlags                           &flags,
                            std::ostream                            &out) 
 {
   AssertThrow (out, ExcIO());
@@ -508,8 +513,10 @@ void DataOutBase::write_dx (const typename std::vector<Patch<dim,spacedim> > &pa
              Assert (false, ExcNotImplemented());
       };
 
-                                  // start with vertices
-  out << "object 1 class array type float rank 1 shape " << spacedim << " items " << n_nodes << " data follows"
+                                  // start with vertices order is
+                                  // lexicographical, x varying
+                                  // fastest
+  out << "object \"vertices\" class array type float rank 1 shape " << spacedim << " items " << n_nodes << " data follows"
       << std::endl;
 
                                   ///////////////////////////////
@@ -613,7 +620,7 @@ void DataOutBase::write_dx (const typename std::vector<Patch<dim,spacedim> > &pa
 
                                   /////////////////////////////////////////
                                   // write cells
-  out << "object 2 class array type int rank 1 shape " << GeometryInfo<dim>::vertices_per_cell
+  out << "object \"cells\" class array type int rank 1 shape " << GeometryInfo<dim>::vertices_per_cell
       << " items " << n_cells << " data follows"
       << std::endl;
   
@@ -644,8 +651,8 @@ void DataOutBase::write_dx (const typename std::vector<Patch<dim,spacedim> > &pa
                  for (unsigned int j=0; j<n_subdivisions; ++j)
                    {
                      out << first_vertex_of_patch+i*(n_subdivisions+1)+j << '\t'
-                         << first_vertex_of_patch+(i+1)*(n_subdivisions+1)+j << '\t'
                          << first_vertex_of_patch+i*(n_subdivisions+1)+j+1 << '\t'
+                         << first_vertex_of_patch+(i+1)*(n_subdivisions+1)+j << '\t'
                          << first_vertex_of_patch+(i+1)*(n_subdivisions+1)+j+1
                          << std::endl;
                    };
@@ -706,12 +713,134 @@ void DataOutBase::write_dx (const typename std::vector<Patch<dim,spacedim> > &pa
   if (dim==3) out << "cubes";
   out << "\"" << std::endl
       << "attribute \"ref\" string \"positions\"" << std::endl;
-  
+
+//TODO:[GK] Patches must be of same size!
+                                  /////////////////////////////
+                                  // write neighbor information
+  if (flags.write_neighbors)
+    {
+      out << "object \"neighbors\" class array type int rank 1 shape "
+         << GeometryInfo<dim>::faces_per_cell
+         << " items " << n_cells
+         << " data follows";
+
+      for (typename std::vector<Patch<dim,spacedim> >::const_iterator
+            patch=patches.begin();
+          patch!=patches.end(); ++patch)
+       {
+         const unsigned int n = patch->n_subdivisions;
+         const unsigned int cells_per_patch
+           = n
+           * ((dim>1) ? n : 1)
+           * ((dim>2) ? n : 1);
+         const unsigned int patch_start = patch->me * cells_per_patch;
+
+         for (unsigned int ix=0;ix<n;++ix)
+           for (unsigned int iy=0;iy<((dim>1) ? n : 1);++iy)
+             for (unsigned int iz=0;iz<((dim>2) ? n : 1);++iz)
+               {
+                 const unsigned int nx = ix*n;
+                 const unsigned int ny = iy;
+                 const unsigned int nz = iz*n*n;
+
+                 out << std::endl;
+                                                  // Direction -x
+                                                  // Last cell in row
+                                                  // of other patch
+                 if (ix==0)
+                   {
+                     const unsigned int nn = patch->neighbors[0];
+                     out << '\t';
+                     if (nn != patch->no_neighbor)
+                       out << (nn*cells_per_patch+ny+nz+n*(n-1));
+                     else
+                       out << "-1";
+                   } else {      
+                     out << '\t'
+                         << patch_start+nx-n+ny+nz;
+                   }
+                                                      // Direction +x
+                                                      // First cell in row
+                                                      // of other patch
+                 if (ix == n-1)
+                   {
+                     const unsigned int nn = patch->neighbors[1];
+                     out << '\t';
+                     if (nn != patch->no_neighbor)
+                       out << (nn*cells_per_patch+ny+nz);
+                     else
+                       out << "-1";
+                   } else {
+                     out << '\t'
+                         << patch_start+nx+n+ny+nz;
+                   }
+                 if (dim<2)
+                   continue;
+                                                  // Direction -y
+                 if (iy==0)
+                   {
+                     const unsigned int nn = patch->neighbors[2];
+                     out << '\t';
+                     if (nn != patch->no_neighbor)
+                       out << (nn*cells_per_patch+nx+nz+1*(n-1));
+                     else
+                       out << "-1";
+                   } else {      
+                     out << '\t'
+                         << patch_start+nx+ny-1+nz;
+                   }
+                                                  // Direction +y
+                 if (iy == n-1)
+                   {
+                     const unsigned int nn = patch->neighbors[3];
+                     out << '\t';
+                     if (nn != patch->no_neighbor)
+                       out << (nn*cells_per_patch+nx+nz);
+                     else
+                       out << "-1";
+                   } else {
+                     out << '\t'
+                         << patch_start+nx+ny+1+nz;
+                   }
+                 if (dim<3)
+                   continue;
+                 
+                                                  // Direction -z
+                 if (iz==0)
+                   {
+                     const unsigned int nn = patch->neighbors[4];
+                     out << '\t';
+                     if (nn != patch->no_neighbor)
+                       out << (nn*cells_per_patch+nx+ny+n*n*(n-1));
+                     else
+                       out << "-1";
+                   } else {      
+                     out << '\t'
+                         << patch_start+nx+ny+nz-n*n;
+                   }
+                                                  // Direction +z
+                 if (iz == n-1)
+                   {
+                     const unsigned int nn = patch->neighbors[4];
+                     out << '\t';
+                     if (nn != patch->no_neighbor)
+                       out << (nn*cells_per_patch+nx+ny);
+                     else
+                       out << "-1";
+                   } else {
+                     out << '\t'
+                         << patch_start+nx+ny+nz+n*n;
+                   }
+               }
+         out << std::endl;       
+       }
+    }
                                   /////////////////////////////
                                   // now write data
   if (n_data_sets != 0)
     {      
-      out << "object 3 class array type float rank 1 shape " << n_data_sets
+      out << "object \"data\" class array type float rank 1 shape "
+         << n_data_sets
          << " items " << n_nodes << " data follows" << std::endl;
 
                                       // loop over all patches
@@ -780,16 +909,19 @@ void DataOutBase::write_dx (const typename std::vector<Patch<dim,spacedim> > &pa
        };
       out << "attribute \"dep\" string \"positions\"" << std::endl;
     } else {
-      out << "object 3 class constantarray type float rank 0 items " << n_nodes << " data follows"
+      out << "object \"data\" class constantarray type float rank 0 items " << n_nodes << " data follows"
          << std::endl << '0' << std::endl;
     }
   
                                   // no model data
   
   out << "object \"deal data\" class field" << std::endl
-      << "component \"positions\" value 1" << std::endl
-      << "component \"connections\" value 2" << std::endl
-      << "component \"data\" value 3" << std::endl;
+      << "component \"positions\" value \"vertices\"" << std::endl
+      << "component \"connections\" value \"cells\"" << std::endl
+      << "component \"data\" value \"data\"" << std::endl;
+
+  if (flags.write_neighbors)
+    out << "component \"neighbors\" value \"neighbors\"" << std::endl;
 
   if (true)
     {

diff --git a/deal.II/deal.II/include/numerics/data_out.h b/deal.II/deal.II/include/numerics/data_out.h
index d74caa2e80f25932feb49710c0cbc474a64d1614..fca43aa905456a1ec3f1279c0ba13ee3d89079eb 100644 (file)
--- a/deal.II/deal.II/include/numerics/data_out.h
+++ b/deal.II/deal.II/include/numerics/data_out.h
@@ -633,6 +633,15 @@ class DataOut : public DataOut_DoFData<dim,dim>
                    << ", is not valid.");
     
   private:
+                                    /**
+                                     * Define a map from cell
+                                     * iterators to patch iterators
+                                     * for finding neighbors.
+                                     */
+    typedef map<typename DoFHandler<dim>::cell_iterator,
+      typename std::vector<typename DataOutBase::Patch<dim> >::iterator>
+      PatchMap;
+    
                                     /**
                                      * All data needed in one thread
                                      * is gathered in the struct
@@ -650,6 +659,7 @@ 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 ()
          {}
     };

diff --git a/deal.II/deal.II/source/numerics/data_out.cc b/deal.II/deal.II/source/numerics/data_out.cc
index e073e301ad6d2710952e4af1baeb81e4d6e96fba..84cbcc1a5525020474b913dbcfae58c335cbdc74 100644 (file)
--- a/deal.II/deal.II/source/numerics/data_out.cc
+++ b/deal.II/deal.II/source/numerics/data_out.cc
@@ -361,6 +361,10 @@ void DataOut<dim>::build_some_patches (Data data)
   for (;cell != dofs->end();)
     {
       Assert (patch != 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);
@@ -467,7 +471,7 @@ void DataOut<dim>::build_patches (const unsigned int n_subdivisions,
                                   // clear the patches array
   if (true)
     {
-      std::vector<DataOutBase::Patch<dim> > dummy;
+      typename std::vector<DataOutBase::Patch<dim> > dummy;
       patches.swap (dummy);
     };
   
@@ -480,6 +484,9 @@ void DataOut<dim>::build_patches (const unsigned int n_subdivisions,
        cell = next_cell(cell))
     ++n_patches;
 
+  
+  PatchMap patch_map;
+    
   std::vector<Data> thread_data(n_threads);
 
                                   // init data for the threads
@@ -490,6 +497,7 @@ 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);
       
@@ -522,6 +530,67 @@ void DataOut<dim>::build_patches (const unsigned int n_subdivisions,
 #else
   build_some_patches(thread_data[0]);
 #endif
+//TODO: New code for neighbors: Parallelize?
+
+                                  // First, number patches
+                                  // consecutively.
+  unsigned int patch_no = 0;
+  typename std::vector<typename DataOutBase::Patch<dim> >::iterator patch;
+  for (patch = patches.begin(); patch != 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.
+  typename PatchMap::iterator map_entry;
+  for (map_entry = patch_map.begin(); map_entry != patch_map.end(); ++map_entry)
+    {
+      typename DoFHandler<dim>::cell_iterator cell = map_entry->first;
+      patch = map_entry->second;
+
+      for (unsigned int i=0;i<GeometryInfo<dim>::faces_per_cell;++i)
+       {
+         if (cell->at_boundary(i))
+           continue;
+         typename DoFHandler<dim>::cell_iterator
+           neighbor = cell->neighbor(i);
+         if (cell->level() != neighbor->level())
+           continue;
+         typename PatchMap::iterator
+           neighbor_entry = patch_map.find(neighbor);
+         if (neighbor_entry == patch_map.end())
+           continue;
+         typename std::vector<typename DataOutBase::Patch<dim> >::iterator
+           neighbor_patch = neighbor_entry->second;
+         
+         switch (dim)
+           {
+             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());
+           }
+       }
+    }
 };