From: Wolfgang Bangerth Date: Wed, 29 May 2019 15:39:57 +0000 (-0600) Subject: Sort implementations of GridTools::cell_measure by dimension. X-Git-Tag: v9.2.0-rc1~1445^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f093fffc447d17ce5d805fbeb5ea7f7c3927f0de;p=dealii.git Sort implementations of GridTools::cell_measure by dimension. --- diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 6293630325..7e1d8b70c8 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -199,6 +199,65 @@ namespace GridTools + template <> + double + cell_measure<2>( + const std::vector> &all_vertices, + const unsigned int (&vertex_indices)[GeometryInfo<2>::vertices_per_cell]) + { + /* + Get the computation of the measure by this little Maple script. We + use the blinear mapping of the unit quad to the real quad. However, + every transformation mapping the unit faces to straight lines should + do. + + Remember that the area of the quad is given by + \int_K 1 dx dy = \int_{\hat K} |det J| d(xi) d(eta) + + # x and y are arrays holding the x- and y-values of the four vertices + # of this cell in real space. + x := array(0..3); + y := array(0..3); + z := array(0..3); + tphi[0] := (1-xi)*(1-eta): + tphi[1] := xi*(1-eta): + tphi[2] := (1-xi)*eta: + tphi[3] := xi*eta: + x_real := sum(x[s]*tphi[s], s=0..3): + y_real := sum(y[s]*tphi[s], s=0..3): + z_real := sum(z[s]*tphi[s], s=0..3): + + Jxi := ; + Jeta := ; + with(VectorCalculus): + J := CrossProduct(Jxi, Jeta); + detJ := sqrt(J[1]^2 + J[2]^2 +J[3]^2); + + # measure := evalf (Int (Int (detJ, xi=0..1, method = _NCrule ) , + eta=0..1, method = _NCrule ) ): # readlib(C): + + # C(measure, optimized); + + additional optimizaton: divide by 2 only one time + */ + + const double x[4] = {all_vertices[vertex_indices[0]](0), + all_vertices[vertex_indices[1]](0), + all_vertices[vertex_indices[2]](0), + all_vertices[vertex_indices[3]](0)}; + + const double y[4] = {all_vertices[vertex_indices[0]](1), + all_vertices[vertex_indices[1]](1), + all_vertices[vertex_indices[2]](1), + all_vertices[vertex_indices[3]](1)}; + + return (-x[1] * y[0] + x[1] * y[3] + y[0] * x[2] + x[0] * y[1] - + x[0] * y[2] - y[1] * x[3] - x[2] * y[3] + x[3] * y[2]) / + 2; + } + + + template <> double cell_measure<3>( @@ -356,65 +415,6 @@ namespace GridTools - template <> - double - cell_measure<2>( - const std::vector> &all_vertices, - const unsigned int (&vertex_indices)[GeometryInfo<2>::vertices_per_cell]) - { - /* - Get the computation of the measure by this little Maple script. We - use the blinear mapping of the unit quad to the real quad. However, - every transformation mapping the unit faces to straight lines should - do. - - Remember that the area of the quad is given by - \int_K 1 dx dy = \int_{\hat K} |det J| d(xi) d(eta) - - # x and y are arrays holding the x- and y-values of the four vertices - # of this cell in real space. - x := array(0..3); - y := array(0..3); - z := array(0..3); - tphi[0] := (1-xi)*(1-eta): - tphi[1] := xi*(1-eta): - tphi[2] := (1-xi)*eta: - tphi[3] := xi*eta: - x_real := sum(x[s]*tphi[s], s=0..3): - y_real := sum(y[s]*tphi[s], s=0..3): - z_real := sum(z[s]*tphi[s], s=0..3): - - Jxi := ; - Jeta := ; - with(VectorCalculus): - J := CrossProduct(Jxi, Jeta); - detJ := sqrt(J[1]^2 + J[2]^2 +J[3]^2); - - # measure := evalf (Int (Int (detJ, xi=0..1, method = _NCrule ) , - eta=0..1, method = _NCrule ) ): # readlib(C): - - # C(measure, optimized); - - additional optimizaton: divide by 2 only one time - */ - - const double x[4] = {all_vertices[vertex_indices[0]](0), - all_vertices[vertex_indices[1]](0), - all_vertices[vertex_indices[2]](0), - all_vertices[vertex_indices[3]](0)}; - - const double y[4] = {all_vertices[vertex_indices[0]](1), - all_vertices[vertex_indices[1]](1), - all_vertices[vertex_indices[2]](1), - all_vertices[vertex_indices[3]](1)}; - - return (-x[1] * y[0] + x[1] * y[3] + y[0] * x[2] + x[0] * y[1] - - x[0] * y[2] - y[1] * x[3] - x[2] * y[3] + x[3] * y[2]) / - 2; - } - - - template BoundingBox compute_bounding_box(const Triangulation &tria)