From: guido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Thu, 13 Dec 2001 18:58:50 +0000 (+0000)
Subject: neighbors for DX in 2D
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=158f979a78378c4a61b7c835d715c5573e7e9bf3;p=dealii-svn.git

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 6dd7367714..5099334abd 100644
--- 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 35be247473..b2c5f77b93 100644
--- 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 38d70727c9..6fc1719f30 100644
--- 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 d74caa2e80..fca43aa905 100644
--- 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 e073e301ad..84cbcc1a55 100644
--- 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());
+	    }
+	}
+    }
 };