From: hartmann Date: Wed, 6 Jul 2005 16:34:54 +0000 (+0000) Subject: New GridTools::cell_measure function. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=327be9b69f606d37376ba546e43d617ca3019930;p=dealii-svn.git 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 77109e4c44..102a2c91c2 100644 --- 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 + static + double cell_measure(const std::vector > &all_vertices, + const int vertex_indices[GeometryInfo::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 55bb4bf14c..f3813c1689 100644 --- 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 > &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 +double +GridTools::cell_measure(const std::vector > &all_vertices, + const int [GeometryInfo::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