From fa65fe3b3e9f703328efa9f1b3c2589321623a36 Mon Sep 17 00:00:00 2001 From: David Wells Date: Sat, 18 Apr 2020 16:08:37 -0400 Subject: [PATCH] Move some non-template functions into a new file. This will let us split up grid_tools.cc with SPLIT_INSTANTIATIONS_COUNT in the next step. --- include/deal.II/grid/grid_tools.h | 13 ++ source/grid/CMakeLists.txt | 1 + source/grid/grid_tools.cc | 296 +----------------------- source/grid/grid_tools.inst.in | 3 +- source/grid/grid_tools_nontemplates.cc | 304 +++++++++++++++++++++++++ 5 files changed, 324 insertions(+), 293 deletions(-) create mode 100644 source/grid/grid_tools_nontemplates.cc diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 1729391102..b8cb641a6e 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -3132,6 +3132,19 @@ namespace GridTools + // This specialization is defined here so that the general template in the + // source file doesn't need to have further 1D overloads for the internal + // functions it calls. + template <> + inline Triangulation<1, 1>::DistortedCellList + fix_up_distorted_child_cells(const Triangulation<1, 1>::DistortedCellList &, + Triangulation<1, 1> &) + { + return {}; + } + + + template void transform(const Predicate & predicate, diff --git a/source/grid/CMakeLists.txt b/source/grid/CMakeLists.txt index f33e978527..f904610a9f 100644 --- a/source/grid/CMakeLists.txt +++ b/source/grid/CMakeLists.txt @@ -37,6 +37,7 @@ SET(_separate_src grid_generator.cc grid_tools.cc grid_tools_cache.cc + grid_tools_nontemplates.cc grid_tools_dof_handlers.cc tria.cc ) diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 3b16e0471a..8857aed1c5 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -187,233 +187,6 @@ namespace GridTools - template <> - double - cell_measure<1>( - const std::vector> &all_vertices, - const unsigned int (&vertex_indices)[GeometryInfo<1>::vertices_per_cell]) - { - return all_vertices[vertex_indices[1]][0] - - all_vertices[vertex_indices[0]][0]; - } - - - - 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>( - const std::vector> &all_vertices, - const unsigned int (&vertex_indices)[GeometryInfo<3>::vertices_per_cell]) - { - // note that this is the - // cell_measure based on the new - // deal.II numbering. When called - // from inside GridReordering make - // sure that you reorder the - // vertex_indices before - 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)}; - - /* - This is the same Maple script as in the barycenter method above - except of that here the shape functions tphi[0]-tphi[7] are ordered - according to the lexicographic numbering. - - x := array(0..7): - y := array(0..7): - z := array(0..7): - tphi[0] := (1-xi)*(1-eta)*(1-zeta): - tphi[1] := xi*(1-eta)*(1-zeta): - tphi[2] := (1-xi)* eta*(1-zeta): - tphi[3] := xi* eta*(1-zeta): - tphi[4] := (1-xi)*(1-eta)*zeta: - tphi[5] := xi*(1-eta)*zeta: - tphi[6] := (1-xi)* eta*zeta: - tphi[7] := xi* eta*zeta: - x_real := sum(x[s]*tphi[s], s=0..7): - y_real := sum(y[s]*tphi[s], s=0..7): - z_real := sum(z[s]*tphi[s], s=0..7): - with (linalg): - J := matrix(3,3, [[diff(x_real, xi), diff(x_real, eta), diff(x_real, - zeta)], [diff(y_real, xi), diff(y_real, eta), diff(y_real, zeta)], - [diff(z_real, xi), diff(z_real, eta), diff(z_real, zeta)]]): - detJ := det (J): - - measure := simplify ( int ( int ( int (detJ, xi=0..1), eta=0..1), - zeta=0..1)): - - readlib(C): - - C(measure, optimized); - - The C code produced by this maple script is further optimized by - hand. In particular, division by 12 is performed only once, not - hundred of times. - */ - - const double t3 = y[3] * x[2]; - const double t5 = z[1] * x[5]; - const double t9 = z[3] * x[2]; - const double t11 = x[1] * y[0]; - const double t14 = x[4] * y[0]; - const double t18 = x[5] * y[7]; - const double t20 = y[1] * x[3]; - const double t22 = y[5] * x[4]; - const double t26 = z[7] * x[6]; - const double t28 = x[0] * y[4]; - const double t34 = - z[3] * x[1] * y[2] + t3 * z[1] - t5 * y[7] + y[7] * x[4] * z[6] + - t9 * y[6] - t11 * z[4] - t5 * y[3] - t14 * z[2] + z[1] * x[4] * y[0] - - t18 * z[3] + t20 * z[0] - t22 * z[0] - y[0] * x[5] * z[4] - t26 * y[3] + - t28 * z[2] - t9 * y[1] - y[1] * x[4] * z[0] - t11 * z[5]; - const double t37 = y[1] * x[0]; - const double t44 = x[1] * y[5]; - const double t46 = z[1] * x[0]; - const double t49 = x[0] * y[2]; - const double t52 = y[5] * x[7]; - const double t54 = x[3] * y[7]; - const double t56 = x[2] * z[0]; - const double t58 = x[3] * y[2]; - const double t64 = -x[6] * y[4] * z[2] - t37 * z[2] + t18 * z[6] - - x[3] * y[6] * z[2] + t11 * z[2] + t5 * y[0] + - t44 * z[4] - t46 * y[4] - t20 * z[7] - t49 * z[6] - - t22 * z[1] + t52 * z[3] - t54 * z[2] - t56 * y[4] - - t58 * z[0] + y[1] * x[2] * z[0] + t9 * y[7] + t37 * z[4]; - const double t66 = x[1] * y[7]; - const double t68 = y[0] * x[6]; - const double t70 = x[7] * y[6]; - const double t73 = z[5] * x[4]; - const double t76 = x[6] * y[7]; - const double t90 = x[4] * z[0]; - const double t92 = x[1] * y[3]; - const double t95 = -t66 * z[3] - t68 * z[2] - t70 * z[2] + t26 * y[5] - - t73 * y[6] - t14 * z[6] + t76 * z[2] - t3 * z[6] + - x[6] * y[2] * z[4] - z[3] * x[6] * y[2] + t26 * y[4] - - t44 * z[3] - x[1] * y[2] * z[0] + x[5] * y[6] * z[4] + - t54 * z[5] + t90 * y[2] - t92 * z[2] + t46 * y[2]; - const double t102 = x[2] * y[0]; - const double t107 = y[3] * x[7]; - const double t114 = x[0] * y[6]; - const double t125 = - y[0] * x[3] * z[2] - z[7] * x[5] * y[6] - x[2] * y[6] * z[4] + - t102 * z[6] - t52 * z[6] + x[2] * y[4] * z[6] - t107 * z[5] - t54 * z[6] + - t58 * z[6] - x[7] * y[4] * z[6] + t37 * z[5] - t114 * z[4] + t102 * z[4] - - z[1] * x[2] * y[0] + t28 * z[6] - y[5] * x[6] * z[4] - - z[5] * x[1] * y[4] - t73 * y[7]; - const double t129 = z[0] * x[6]; - const double t133 = y[1] * x[7]; - const double t145 = y[1] * x[5]; - const double t156 = t90 * y[6] - t129 * y[4] + z[7] * x[2] * y[6] - - t133 * z[5] + x[5] * y[3] * z[7] - t26 * y[2] - - t70 * z[3] + t46 * y[3] + z[5] * x[7] * y[4] + - z[7] * x[3] * y[6] - t49 * z[4] + t145 * z[7] - - x[2] * y[7] * z[6] + t70 * z[5] + t66 * z[5] - - z[7] * x[4] * y[6] + t18 * z[4] + x[1] * y[4] * z[0]; - const double t160 = x[5] * y[4]; - const double t165 = z[1] * x[7]; - const double t178 = z[1] * x[3]; - const double t181 = - t107 * z[6] + t22 * z[7] + t76 * z[3] + t160 * z[1] - x[4] * y[2] * z[6] + - t70 * z[4] + t165 * y[5] + x[7] * y[2] * z[6] - t76 * z[5] - t76 * z[4] + - t133 * z[3] - t58 * z[1] + y[5] * x[0] * z[4] + t114 * z[2] - t3 * z[7] + - t20 * z[2] + t178 * y[7] + t129 * y[2]; - const double t207 = t92 * z[7] + t22 * z[6] + z[3] * x[0] * y[2] - - x[0] * y[3] * z[2] - z[3] * x[7] * y[2] - t165 * y[3] - - t9 * y[0] + t58 * z[7] + y[3] * x[6] * z[2] + - t107 * z[2] + t73 * y[0] - x[3] * y[5] * z[7] + - t3 * z[0] - t56 * y[6] - z[5] * x[0] * y[4] + - t73 * y[1] - t160 * z[6] + t160 * z[0]; - const double t228 = -t44 * z[7] + z[5] * x[6] * y[4] - t52 * z[4] - - t145 * z[4] + t68 * z[4] + t92 * z[5] - t92 * z[0] + - t11 * z[3] + t44 * z[0] + t178 * y[5] - t46 * y[5] - - t178 * y[0] - t145 * z[0] - t20 * z[5] - t37 * z[3] - - t160 * z[7] + t145 * z[3] + x[4] * y[6] * z[2]; - - return (t34 + t64 + t95 + t125 + t156 + t181 + t207 + t228) / 12.; - } - - template Vector compute_aspect_ratio_of_cells(const Triangulation &triangulation, @@ -631,8 +404,6 @@ namespace GridTools return SubCellData(); } }; - - } // namespace @@ -989,26 +760,6 @@ namespace GridTools }; - // the following class is only - // needed in 2d, so avoid trouble - // with compilers warning otherwise - class Rotate2d - { - public: - explicit Rotate2d(const double angle) - : angle(angle) - {} - Point<2> - operator()(const Point<2> &p) const - { - return {std::cos(angle) * p(0) - std::sin(angle) * p(1), - std::sin(angle) * p(0) + std::cos(angle) * p(1)}; - } - - private: - const double angle; - }; - // Transformation to rotate around one of the cartesian axes. class Rotate3d { @@ -1068,25 +819,6 @@ namespace GridTools } - template <> - void - rotate(const double angle, Triangulation<2> &triangulation) - { - transform(Rotate2d(angle), triangulation); - } - - template <> - void - rotate(const double angle, Triangulation<3> &triangulation) - { - (void)angle; - (void)triangulation; - - AssertThrow( - false, ExcMessage("GridTools::rotate() is not available for dim = 3.")); - } - - template void rotate(const double angle, @@ -1141,19 +873,6 @@ namespace GridTools } // namespace - - // Implementation for 1D only - template <> - void - laplace_transform(const std::map> &, - Triangulation<1> &, - const Function<1> *, - const bool) - { - Assert(false, ExcNotImplemented()); - } - - // Implementation for dimensions except 1 template void @@ -1162,6 +881,9 @@ namespace GridTools const Function * coefficient, const bool solve_for_absolute_positions) { + if (dim == 1) + Assert(false, ExcNotImplemented()); + // first provide everything that is needed for solving a Laplace // equation. FE_Q q1(1); @@ -3748,16 +3470,6 @@ namespace GridTools - void - fix_up_faces(const dealii::Triangulation<1, 1>::cell_iterator &, - std::integral_constant, - std::integral_constant) - { - // nothing to do for the faces of cells in 1d - } - - - // possibly fix up the faces of a cell by moving around its mid-points template void @@ -3810,6 +3522,8 @@ namespace GridTools &distorted_cells, Triangulation & /*triangulation*/) { + static_assert(dim != 1 && spacedim != 1, + "This function is only valid when dim != 1 or spacedim != 1."); typename Triangulation::DistortedCellList unfixable_subset; // loop over all cells that we have to fix up diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index ff06ead6ba..824c610b40 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -357,8 +357,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const Function * coefficient, const bool); -# endif - template Triangulation::DistortedCellList fix_up_distorted_child_cells( @@ -367,6 +365,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) &distorted_cells, Triangulation &triangulation); +# endif # endif \} #endif diff --git a/source/grid/grid_tools_nontemplates.cc b/source/grid/grid_tools_nontemplates.cc new file mode 100644 index 0000000000..706dbc4d42 --- /dev/null +++ b/source/grid/grid_tools_nontemplates.cc @@ -0,0 +1,304 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2001 - 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#include + +#include + +#include + +// GridTools functions that are template specializations (i.e., only compiled +// once without expand_instantiations) + +DEAL_II_NAMESPACE_OPEN + + +namespace GridTools +{ + template <> + double + cell_measure<1>( + const std::vector> &all_vertices, + const unsigned int (&vertex_indices)[GeometryInfo<1>::vertices_per_cell]) + { + return all_vertices[vertex_indices[1]][0] - + all_vertices[vertex_indices[0]][0]; + } + + + + 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>( + const std::vector> &all_vertices, + const unsigned int (&vertex_indices)[GeometryInfo<3>::vertices_per_cell]) + { + // note that this is the + // cell_measure based on the new + // deal.II numbering. When called + // from inside GridReordering make + // sure that you reorder the + // vertex_indices before + 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)}; + + /* + This is the same Maple script as in the barycenter method above + except of that here the shape functions tphi[0]-tphi[7] are ordered + according to the lexicographic numbering. + + x := array(0..7): + y := array(0..7): + z := array(0..7): + tphi[0] := (1-xi)*(1-eta)*(1-zeta): + tphi[1] := xi*(1-eta)*(1-zeta): + tphi[2] := (1-xi)* eta*(1-zeta): + tphi[3] := xi* eta*(1-zeta): + tphi[4] := (1-xi)*(1-eta)*zeta: + tphi[5] := xi*(1-eta)*zeta: + tphi[6] := (1-xi)* eta*zeta: + tphi[7] := xi* eta*zeta: + x_real := sum(x[s]*tphi[s], s=0..7): + y_real := sum(y[s]*tphi[s], s=0..7): + z_real := sum(z[s]*tphi[s], s=0..7): + with (linalg): + J := matrix(3,3, [[diff(x_real, xi), diff(x_real, eta), diff(x_real, + zeta)], [diff(y_real, xi), diff(y_real, eta), diff(y_real, zeta)], + [diff(z_real, xi), diff(z_real, eta), diff(z_real, zeta)]]): + detJ := det (J): + + measure := simplify ( int ( int ( int (detJ, xi=0..1), eta=0..1), + zeta=0..1)): + + readlib(C): + + C(measure, optimized); + + The C code produced by this maple script is further optimized by + hand. In particular, division by 12 is performed only once, not + hundred of times. + */ + + const double t3 = y[3] * x[2]; + const double t5 = z[1] * x[5]; + const double t9 = z[3] * x[2]; + const double t11 = x[1] * y[0]; + const double t14 = x[4] * y[0]; + const double t18 = x[5] * y[7]; + const double t20 = y[1] * x[3]; + const double t22 = y[5] * x[4]; + const double t26 = z[7] * x[6]; + const double t28 = x[0] * y[4]; + const double t34 = + z[3] * x[1] * y[2] + t3 * z[1] - t5 * y[7] + y[7] * x[4] * z[6] + + t9 * y[6] - t11 * z[4] - t5 * y[3] - t14 * z[2] + z[1] * x[4] * y[0] - + t18 * z[3] + t20 * z[0] - t22 * z[0] - y[0] * x[5] * z[4] - t26 * y[3] + + t28 * z[2] - t9 * y[1] - y[1] * x[4] * z[0] - t11 * z[5]; + const double t37 = y[1] * x[0]; + const double t44 = x[1] * y[5]; + const double t46 = z[1] * x[0]; + const double t49 = x[0] * y[2]; + const double t52 = y[5] * x[7]; + const double t54 = x[3] * y[7]; + const double t56 = x[2] * z[0]; + const double t58 = x[3] * y[2]; + const double t64 = -x[6] * y[4] * z[2] - t37 * z[2] + t18 * z[6] - + x[3] * y[6] * z[2] + t11 * z[2] + t5 * y[0] + + t44 * z[4] - t46 * y[4] - t20 * z[7] - t49 * z[6] - + t22 * z[1] + t52 * z[3] - t54 * z[2] - t56 * y[4] - + t58 * z[0] + y[1] * x[2] * z[0] + t9 * y[7] + t37 * z[4]; + const double t66 = x[1] * y[7]; + const double t68 = y[0] * x[6]; + const double t70 = x[7] * y[6]; + const double t73 = z[5] * x[4]; + const double t76 = x[6] * y[7]; + const double t90 = x[4] * z[0]; + const double t92 = x[1] * y[3]; + const double t95 = -t66 * z[3] - t68 * z[2] - t70 * z[2] + t26 * y[5] - + t73 * y[6] - t14 * z[6] + t76 * z[2] - t3 * z[6] + + x[6] * y[2] * z[4] - z[3] * x[6] * y[2] + t26 * y[4] - + t44 * z[3] - x[1] * y[2] * z[0] + x[5] * y[6] * z[4] + + t54 * z[5] + t90 * y[2] - t92 * z[2] + t46 * y[2]; + const double t102 = x[2] * y[0]; + const double t107 = y[3] * x[7]; + const double t114 = x[0] * y[6]; + const double t125 = + y[0] * x[3] * z[2] - z[7] * x[5] * y[6] - x[2] * y[6] * z[4] + + t102 * z[6] - t52 * z[6] + x[2] * y[4] * z[6] - t107 * z[5] - t54 * z[6] + + t58 * z[6] - x[7] * y[4] * z[6] + t37 * z[5] - t114 * z[4] + t102 * z[4] - + z[1] * x[2] * y[0] + t28 * z[6] - y[5] * x[6] * z[4] - + z[5] * x[1] * y[4] - t73 * y[7]; + const double t129 = z[0] * x[6]; + const double t133 = y[1] * x[7]; + const double t145 = y[1] * x[5]; + const double t156 = t90 * y[6] - t129 * y[4] + z[7] * x[2] * y[6] - + t133 * z[5] + x[5] * y[3] * z[7] - t26 * y[2] - + t70 * z[3] + t46 * y[3] + z[5] * x[7] * y[4] + + z[7] * x[3] * y[6] - t49 * z[4] + t145 * z[7] - + x[2] * y[7] * z[6] + t70 * z[5] + t66 * z[5] - + z[7] * x[4] * y[6] + t18 * z[4] + x[1] * y[4] * z[0]; + const double t160 = x[5] * y[4]; + const double t165 = z[1] * x[7]; + const double t178 = z[1] * x[3]; + const double t181 = + t107 * z[6] + t22 * z[7] + t76 * z[3] + t160 * z[1] - x[4] * y[2] * z[6] + + t70 * z[4] + t165 * y[5] + x[7] * y[2] * z[6] - t76 * z[5] - t76 * z[4] + + t133 * z[3] - t58 * z[1] + y[5] * x[0] * z[4] + t114 * z[2] - t3 * z[7] + + t20 * z[2] + t178 * y[7] + t129 * y[2]; + const double t207 = t92 * z[7] + t22 * z[6] + z[3] * x[0] * y[2] - + x[0] * y[3] * z[2] - z[3] * x[7] * y[2] - t165 * y[3] - + t9 * y[0] + t58 * z[7] + y[3] * x[6] * z[2] + + t107 * z[2] + t73 * y[0] - x[3] * y[5] * z[7] + + t3 * z[0] - t56 * y[6] - z[5] * x[0] * y[4] + + t73 * y[1] - t160 * z[6] + t160 * z[0]; + const double t228 = -t44 * z[7] + z[5] * x[6] * y[4] - t52 * z[4] - + t145 * z[4] + t68 * z[4] + t92 * z[5] - t92 * z[0] + + t11 * z[3] + t44 * z[0] + t178 * y[5] - t46 * y[5] - + t178 * y[0] - t145 * z[0] - t20 * z[5] - t37 * z[3] - + t160 * z[7] + t145 * z[3] + x[4] * y[6] * z[2]; + + return (t34 + t64 + t95 + t125 + t156 + t181 + t207 + t228) / 12.; + } + + + + namespace + { + // the following class is only + // needed in 2d, so avoid trouble + // with compilers warning otherwise + class Rotate2d + { + public: + explicit Rotate2d(const double angle) + : angle(angle) + {} + Point<2> + operator()(const Point<2> &p) const + { + return {std::cos(angle) * p(0) - std::sin(angle) * p(1), + std::sin(angle) * p(0) + std::cos(angle) * p(1)}; + } + + private: + const double angle; + }; + } // namespace + + + + template <> + void + rotate(const double angle, Triangulation<2> &triangulation) + { + transform(Rotate2d(angle), triangulation); + } + + + + template <> + void + rotate(const double angle, Triangulation<3> &triangulation) + { + (void)angle; + (void)triangulation; + + AssertThrow( + false, ExcMessage("GridTools::rotate() is not available for dim = 3.")); + } +} /* namespace GridTools */ + +DEAL_II_NAMESPACE_CLOSE -- 2.39.5