New GridTools::cell_measure function.

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

diff --git a/deal.II/deal.II/include/grid/grid_tools.h b/deal.II/deal.II/include/grid/grid_tools.h
index 77109e4c4427bc2e66f33219e713646beea471ca..102a2c91c2855078b1247cd0ed057d9eee3652a8 100644 (file)
--- a/deal.II/deal.II/include/grid/grid_tools.h
+++ b/deal.II/deal.II/include/grid/grid_tools.h
@@ -55,6 +55,18 @@ class GridTools
     static
     double diameter (const Triangulation<1> &tria);
 
+                                    /**
+                                     * Return the measure of a cell
+                                     * represented by a subset of
+                                     * vertices in @p all_vertices
+                                     * which is specified by @p
+                                     * vertex_indices.
+                                     */
+    template <int dim>
+    static
+    double cell_measure(const std::vector<Point<dim> > &all_vertices,
+                       const int vertex_indices[GeometryInfo<dim>::vertices_per_cell]);
+
                                     /**
                                      * Transform the vertices of the
                                      * given triangulation by

diff --git a/deal.II/deal.II/source/grid/grid_tools.cc b/deal.II/deal.II/source/grid/grid_tools.cc
index 55bb4bf14c4a537efac66fe8bc58a554a8c7ec6c..f3813c16892635f95c49a2e5c5bbaf1cdd7f4fd8 100644 (file)
--- a/deal.II/deal.II/source/grid/grid_tools.cc
+++ b/deal.II/deal.II/source/grid/grid_tools.cc
@@ -102,6 +102,110 @@ GridTools::diameter (const Triangulation<1> &tria)
 
 
 
+#if deal_II_dimension == 3
+
+template <>
+double
+GridTools::cell_measure(const std::vector<Point<3> > &all_vertices,
+                       const int vertex_indices[GeometryInfo<3>::vertices_per_cell])
+{
+  const double x[8] = { all_vertices[vertex_indices[0]](0),
+                       all_vertices[vertex_indices[1]](0),
+                       all_vertices[vertex_indices[2]](0),
+                       all_vertices[vertex_indices[3]](0),
+                       all_vertices[vertex_indices[4]](0),
+                       all_vertices[vertex_indices[5]](0),
+                       all_vertices[vertex_indices[6]](0),
+                       all_vertices[vertex_indices[7]](0)   };
+  const double y[8] = { all_vertices[vertex_indices[0]](1),
+                       all_vertices[vertex_indices[1]](1),
+                       all_vertices[vertex_indices[2]](1),
+                       all_vertices[vertex_indices[3]](1),
+                       all_vertices[vertex_indices[4]](1),
+                       all_vertices[vertex_indices[5]](1),
+                       all_vertices[vertex_indices[6]](1),
+                       all_vertices[vertex_indices[7]](1)  };
+  const double z[8] = { all_vertices[vertex_indices[0]](2),
+                       all_vertices[vertex_indices[1]](2),
+                       all_vertices[vertex_indices[2]](2),
+                       all_vertices[vertex_indices[3]](2),
+                       all_vertices[vertex_indices[4]](2),
+                       all_vertices[vertex_indices[5]](2),
+                       all_vertices[vertex_indices[6]](2),
+                       all_vertices[vertex_indices[7]](2)  };
+
+/*
+  Get the computation of the measure by the same little Maple script
+  as in TriaObjectAccessor<3, 3>::barycenter(), see tria_accessor.cc.
+*/
+  double s1, s2, s3;
+  
+  s3 = -x[2]*y[7]*z[3]/12.0+x[1]*y[4]*z[0]/12.0-x[0]*y[3]*z[7]/12.0-z[1]*x
+       [2]*y[0]/12.0+x[5]*y[4]*z[0]/12.0+x[1]*y[5]*z[0]/12.0-x[1]*y[2]*z[0]/12.0-y[2]*
+       x[3]*z[7]/12.0-x[1]*y[0]*z[4]/12.0-y[1]*x[0]*z[2]/12.0-y[0]*x[7]*z[3]/12.0+y[5]
+       *x[4]*z[7]/12.0-x[0]*y[3]*z[4]/12.0-y[1]*x[0]*z[3]/12.0-z[0]*x[7]*y[4]/12.0+x
+       [1]*y[0]*z[3]/12.0-x[2]*y[3]*z[0]/12.0-y[1]*x[2]*z[5]/12.0;
+  s2 = s3+y[0]*x[3]*z[4]/12.0+y[1]*x[2]*z[0]/12.0+y[2]*x[7]*z[3]/12.0+z[1]*
+       x[4]*y[0]/12.0-y[5]*x[4]*z[0]/12.0-y[1]*x[5]*z[0]/12.0+x[5]*y[7]*z[4]/12.0+y[1]
+       *x[3]*z[0]/12.0+x[0]*y[4]*z[7]/12.0+z[0]*x[7]*y[3]/12.0-z[2]*x[3]*y[0]/12.0+y
+       [1]*x[0]*z[5]/12.0-x[1]*y[5]*z[2]/12.0+y[1]*x[0]*z[4]/12.0+z[1]*x[2]*y[5]/12.0-
+       x[5]*y[0]*z[4]/12.0+y[2]*x[3]*z[1]/12.0+x[1]*y[0]*z[2]/12.0;
+  s3 = -z[1]*x[0]*y[4]/12.0+z[1]*x[0]*y[2]/12.0-x[0]*y[7]*z[4]/12.0+z[0]*x
+       [4]*y[7]/12.0+x[2]*y[3]*z[7]/12.0+y[2]*x[6]*z[3]/12.0+y[1]*x[6]*z[2]/12.0-y[2]*
+       x[1]*z[3]/12.0+x[1]*y[2]*z[5]/12.0+y[0]*x[3]*z[7]/12.0+y[2]*x[3]*z[0]/12.0-y[0]
+       *x[4]*z[3]/12.0-y[2]*x[0]*z[3]/12.0-y[1]*x[4]*z[0]/12.0+z[1]*x[5]*y[0]/12.0+x
+       [1]*y[2]*z[6]/12.0-x[1]*y[0]*z[5]/12.0-x[2]*y[3]*z[1]/12.0;
+  s1 = s3-z[0]*x[3]*y[4]/12.0+z[1]*x[0]*y[3]/12.0-z[1]*x[5]*y[2]/12.0+z[1]*
+       x[2]*y[6]/12.0-z[2]*x[3]*y[1]/12.0+z[2]*x[3]*y[6]/12.0+x[2]*y[1]*z[3]/12.0-z[1]
+       *x[3]*y[0]/12.0+x[2]*y[3]*z[6]/12.0-z[0]*x[3]*y[7]/12.0+y[5]*x[0]*z[4]/12.0-z
+       [1]*x[6]*y[2]/12.0-z[2]*x[6]*y[3]/12.0+z[2]*x[0]*y[3]/12.0+z[2]*x[3]*y[7]/12.0+
+       z[2]*x[1]*y[3]/12.0+y[1]*x[5]*z[2]/12.0-z[2]*x[7]*y[3]/12.0+s2;
+  s3 = x[0]*y[7]*z[3]/12.0-z[1]*x[0]*y[5]/12.0-x[1]*y[3]*z[0]/12.0-x[1]*y
+       [6]*z[2]/12.0-x[2]*y[6]*z[3]/12.0-x[5]*y[4]*z[7]/12.0+z[0]*x[4]*y[3]/12.0+x[0]*
+       y[4]*z[3]/12.0-x[6]*y[5]*z[7]/12.0-x[6]*y[7]*z[2]/12.0+z[1]*x[6]*y[5]/12.0+y[6]
+       *x[4]*z[7]/12.0+y[1]*x[5]*z[6]/12.0+x[1]*y[6]*z[5]/12.0-y[5]*x[7]*z[4]/12.0+y
+       [0]*x[7]*z[4]/12.0-y[0]*x[4]*z[7]/12.0-z[5]*x[0]*y[4]/12.0;
+  s2 = s3+z[6]*x[3]*y[7]/12.0+z[6]*x[7]*y[4]/12.0-z[6]*x[4]*y[7]/12.0-z[6]*
+       x[7]*y[3]/12.0+y[6]*x[5]*z[7]/12.0-y[6]*x[7]*z[5]/12.0-x[2]*y[5]*z[6]/12.0+z[6]
+       *x[2]*y[7]/12.0+z[6]*x[7]*y[5]/12.0-z[6]*x[5]*y[7]/12.0-z[6]*x[7]*y[2]/12.0+x
+       [6]*y[7]*z[5]/12.0+z[5]*x[7]*y[4]/12.0-z[5]*x[4]*y[7]/12.0+y[5]*x[1]*z[4]/12.0-
+       y[5]*x[6]*z[4]/12.0-y[5]*x[4]*z[1]/12.0+y[5]*x[4]*z[6]/12.0;
+  s3 = y[6]*x[7]*z[3]/12.0+x[6]*y[5]*z[2]/12.0-x[5]*y[6]*z[2]/12.0-x[6]*y
+       [2]*z[5]/12.0+x[5]*y[2]*z[6]/12.0+x[6]*y[2]*z[7]/12.0+z[5]*x[4]*y[0]/12.0-y[6]*
+       x[2]*z[7]/12.0+x[2]*y[6]*z[5]/12.0+y[6]*x[7]*z[2]/12.0-y[6]*x[3]*z[7]/12.0-x[7]
+       *y[4]*z[3]/12.0-x[6]*y[4]*z[7]/12.0-x[5]*y[1]*z[4]/12.0+x[6]*y[7]*z[4]/12.0+x
+       [7]*y[3]*z[4]/12.0-z[1]*x[5]*y[6]/12.0+x[3]*y[4]*z[7]/12.0+x[5]*y[6]*z[4]/12.0;
+
+                                  // the vertices entered into the
+                                  // above field and the vertices
+                                  // used by the maple script use
+                                  // different orderings (i.e. one is
+                                  // the inside-out orientation of
+                                  // the other one). the measure is
+                                  // thus negative. Thus take the
+                                  // negative value of the result
+  return -(s3+x[5]*y[4]*z[1]/12.0-x[5]*y[4]*z[6]/12.0-x[3]*y[7]*z[4]/12.0+x[6]*
+    y[3]*z[7]/12.0-y[1]*x[6]*z[5]/12.0-z[5]*x[1]*y[4]/12.0+z[5]*x[6]*y[4]/12.0+z[5]
+    *x[4]*y[1]/12.0-z[5]*x[4]*y[6]/12.0-x[6]*y[7]*z[3]/12.0-x[4]*y[3]*z[7]/12.0+x
+    [4]*y[7]*z[3]/12.0-x[1]*y[5]*z[6]/12.0-y[6]*x[7]*z[4]/12.0-y[1]*x[2]*z[6]/12.0-
+    y[2]*x[3]*z[6]/12.0+s2+s1+x[2]*y[0]*z[3]/12.0);
+}
+
+#else
+
+template <int dim>
+double
+GridTools::cell_measure(const std::vector<Point<dim> > &all_vertices,
+                       const int [GeometryInfo<dim>::vertices_per_cell])
+{
+  Assert(false, ExcNotImplemented());
+  return 0.;
+}
+
+#endif
+
+
+
 // define some transformations in an anonymous namespace
 #ifdef DEAL_II_ANON_NAMESPACE_BOGUS_WARNING
   namespace TRANS