From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Fri, 29 Sep 2000 15:16:11 +0000 (+0000)
Subject: Write a function that deletes unused vertices.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e6a4625394e6fd94d59056d647d29a598466e3b3;p=dealii-svn.git

Write a function that deletes unused vertices.


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

diff --git a/deal.II/deal.II/source/grid/grid_in.cc b/deal.II/deal.II/source/grid/grid_in.cc
index 7b351e38c3..d38f6a1b25 100644
--- a/deal.II/deal.II/source/grid/grid_in.cc
+++ b/deal.II/deal.II/source/grid/grid_in.cc
@@ -14,8 +14,144 @@
 
 #include <grid/grid_in.h>
 #include <grid/tria.h>
+#include <grid/grid_reordering.h>
 
 #include <map>
+#include <algorithm>
+
+
+template <int dim>
+void
+delete_unused_vertices (vector<Point<dim> >    &vertices,
+			vector<CellData<dim> > &cells,
+			SubCellData            &subcelldata)
+{
+				   // first check which vertices are
+				   // actually used
+  vector<bool> vertex_used (vertices.size(), false);
+  for (unsigned int c=0; c<cells.size(); ++c)
+    for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+      vertex_used[cells[c].vertices[v]] = true;
+
+				   // then renumber the vertices that
+				   // are actually used in the same
+				   // order as they were beforehand
+  const unsigned int invalid_vertex = static_cast<unsigned int>(-1);
+  vector<unsigned int> new_vertex_numbers (vertices.size(), invalid_vertex);
+  unsigned int next_free_number = 0;
+  for (unsigned int i=0; i<vertices.size(); ++i)
+    if (vertex_used[i] == true)
+      {
+	new_vertex_numbers[i] = next_free_number;
+	++next_free_number;
+      };
+
+				   // next replace old vertex numbers
+				   // by the new ones
+  for (unsigned int c=0; c<cells.size(); ++c)
+    for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+      cells[c].vertices[v] = new_vertex_numbers[cells[c].vertices[v]];
+
+				   // same for boundary data
+  for (unsigned int c=0; c<subcelldata.boundary_lines.size(); ++c)
+    for (unsigned int v=0; v<GeometryInfo<1>::vertices_per_cell; ++v)
+      subcelldata.boundary_lines[c].vertices[v]
+	= new_vertex_numbers[subcelldata.boundary_lines[c].vertices[v]];
+  for (unsigned int c=0; c<subcelldata.boundary_quads.size(); ++c)
+    for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
+      subcelldata.boundary_quads[c].vertices[v]
+	= new_vertex_numbers[subcelldata.boundary_quads[c].vertices[v]];
+
+				   // finally copy over the vertices
+				   // which we really need to a new
+				   // array and replace the old one by
+				   // the new one
+  vector<Point<dim> > tmp;
+  tmp.reserve (count(vertex_used.begin(), vertex_used.end(), true));
+  for (unsigned int v=0; v<vertices.size(); ++v)
+    if (vertex_used[v] == true)
+      tmp.push_back (vertices[v]);
+  swap (vertices, tmp);
+};
+
+
+
+#include <fstream>
+void
+debug_output (const vector<CellData<3> > &cells,
+	      const vector<Point<3> >    &vertices,
+	      const string               &filename);
+void
+debug_output (const vector<CellData<1> > &cells,
+	      const vector<Point<1> >    &vertices,
+	      const string               &filename);
+void
+debug_output (const vector<CellData<2> > &cells,
+	      const vector<Point<2> >    &vertices,
+	      const string               &filename)
+{
+  double min_x = vertices[cells[0].vertices[0]](0),
+	 max_x = vertices[cells[0].vertices[0]](0),
+	 min_y = vertices[cells[0].vertices[0]](1),
+	 max_y = vertices[cells[0].vertices[0]](1);
+  
+  ofstream x(filename.c_str());
+  for (unsigned int i=0; i<cells.size(); ++i)
+    {
+      for (unsigned int v=0; v<4; ++v)
+	{
+	  const Point<2> &p = vertices[cells[i].vertices[v]];
+	  
+	  if (p(0) < min_x)
+	    min_x = p(0);
+	  if (p(0) > max_x)
+	    max_x = p(0);
+	  if (p(1) < min_y)
+	    min_y = p(1);
+	  if (p(1) > max_y)
+	    max_y = p(1);
+	};
+
+      x << "# cell " << i << endl;
+      Point<2> center;
+      for (unsigned int f=0; f<4; ++f)
+	center += vertices[cells[i].vertices[f]];
+      center /= 4;
+
+      x << "set label \"" << i << "\" at "
+	<< center(0) << ',' << center(1)
+	<< " center"
+	<< endl;
+
+				       // first two line right direction
+      for (unsigned int f=0; f<2; ++f)
+	x << "set arrow from "
+	  << vertices[cells[i].vertices[f]](0) << ',' << vertices[cells[i].vertices[f]](1)
+	  << " to "
+	  << vertices[cells[i].vertices[(f+1)%4]](0) << ',' << vertices[cells[i].vertices[(f+1)%4]](1)
+	  << endl;
+				       // other two lines reverse direction
+      for (unsigned int f=2; f<4; ++f)
+	x << "set arrow from "
+	  << vertices[cells[i].vertices[(f+1)%4]](0) << ',' << vertices[cells[i].vertices[(f+1)%4]](1)
+	  << " to "
+	  << vertices[cells[i].vertices[f]](0) << ',' << vertices[cells[i].vertices[f]](1)
+	  << endl;
+      x << endl;
+    };
+  
+
+  x << endl
+    << "pl [" << min_x << ':' << max_x << "]["
+    << min_y << ':' << max_y <<  "] "
+    << min_y << endl
+    << "pause -1" << endl;
+};
+
+      
+
+
+
 
 
 template <int dim>
@@ -109,7 +245,7 @@ void GridIn<dim>::read_ucd (istream &in)
 	    else 
 	      {
 						 // no such vertex index
-		Assert (false, ExcInvalidVertexIndex(cell, cells.back().vertices[i]));
+		AssertThrow (false, ExcInvalidVertexIndex(cell, cells.back().vertices[i]));
 		cells.back().vertices[i] = -1;
 	      };
 	}
@@ -133,15 +269,15 @@ void GridIn<dim>::read_ucd (istream &in)
 	    else 
 	      {
 						 // no such vertex index
-		Assert (false,
-			ExcInvalidVertexIndex(cell,
-					      subcelldata.boundary_lines.back().vertices[i]));
+		AssertThrow (false,
+			     ExcInvalidVertexIndex(cell,
+						   subcelldata.boundary_lines.back().vertices[i]));
 		subcelldata.boundary_lines.back().vertices[i] = -1;
 	      };
 	  }
 	else
 					   // cannot read this
-	  Assert (false, ExcUnknownIdentifier(cell_type));
+	  AssertThrow (false, ExcUnknownIdentifier(cell_type));
     };
 
 				   // check that no forbidden arrays are used
@@ -203,7 +339,7 @@ void GridIn<dim>::read_dbmesh (istream &in)
   AssertThrow (line=="Vertices", ExcInvalidDBMESHInput(line));
   
   unsigned int n_vertices;
-  int dummy;
+  double dummy;
   
   in >> n_vertices;
   vector<Point<dim> >     vertices (n_vertices);
@@ -215,6 +351,7 @@ void GridIn<dim>::read_dbmesh (istream &in)
 				       // read Ref phi_i, whatever that may be
       in >> dummy;
     };
+  AssertThrow (in, ExcInvalidDBMeshFormat());
 
   skip_empty_lines(in);
 
@@ -233,6 +370,7 @@ void GridIn<dim>::read_dbmesh (istream &in)
 				       // read Ref phi_i, whatever that may be
       in >> dummy;
     };
+  AssertThrow (in, ExcInvalidDBMeshFormat());
 
   skip_empty_lines(in);
 
@@ -253,6 +391,7 @@ void GridIn<dim>::read_dbmesh (istream &in)
 				       // read Ref phi_i, whatever that may be
       in >> dummy;
     };
+  AssertThrow (in, ExcInvalidDBMeshFormat());
 
   skip_empty_lines(in);
 
@@ -274,12 +413,19 @@ void GridIn<dim>::read_dbmesh (istream &in)
       for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
 	{
 	  in >> cells.back().vertices[i];
+	  
+	  AssertThrow ((cells.back().vertices[i] >= 1)
+		       &&
+		       (static_cast<unsigned int>(cells.back().vertices[i]) <= vertices.size()),
+		       ExcInvalidVertexIndex(cell, cells.back().vertices[i]));
+		  
 	  --cells.back().vertices[i];
 	};
 
 				       // read and discard Ref phi_i
       in >> dummy;
     };
+  AssertThrow (in, ExcInvalidDBMeshFormat());
 
   skip_empty_lines(in);
   
@@ -289,7 +435,7 @@ void GridIn<dim>::read_dbmesh (istream &in)
 				   // clue what they mean. skip them
 				   // all and leave the interpretation
 				   // to other implementors...
-  while (getline(in,line), line.find("End")==string::npos);
+  while (getline(in,line), ((line.find("End")==string::npos) && (in)));
 				   // ok, so we are not at the end of
 				   // the file, that's it, mostly
 
@@ -299,6 +445,12 @@ void GridIn<dim>::read_dbmesh (istream &in)
 
   AssertThrow (in, ExcIO());
 
+				   // do some clean-up on vertices
+  delete_unused_vertices (vertices, cells, subcelldata);
+  
+  debug_output (cells, vertices, "x1");
+  GridReordering<dim>::reorder_cells (cells);
+  debug_output (cells, vertices, "x2");
   tria->create_triangulation (vertices, cells, subcelldata);
 };