]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Move some non-template functions into a new file.
authorDavid Wells <drwells@email.unc.edu>
Sat, 18 Apr 2020 20:08:37 +0000 (16:08 -0400)
committerDavid Wells <drwells@email.unc.edu>
Sun, 19 Apr 2020 14:17:19 +0000 (10:17 -0400)
This will let us split up grid_tools.cc with SPLIT_INSTANTIATIONS_COUNT
in the next step.

include/deal.II/grid/grid_tools.h
source/grid/CMakeLists.txt
source/grid/grid_tools.cc
source/grid/grid_tools.inst.in
source/grid/grid_tools_nontemplates.cc [new file with mode: 0644]

index 1729391102b894199fe3d7e40d6f436e318e0e8f..b8cb641a6e05416eaeddca17ca0b439a370e3541 100644 (file)
@@ -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 <int dim, typename Predicate, int spacedim>
   void
   transform(const Predicate &             predicate,
index f33e978527f38c3477edbeec3bca0fa9261dcb3a..f904610a9f3af0ccbabc85ba792a96918d6e4e5a 100644 (file)
@@ -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
   )
index 3b16e0471afb00821bcbf49591914d98fa2321e7..8857aed1c510ebeebb24517c056fe0261aefb744 100644 (file)
@@ -187,233 +187,6 @@ namespace GridTools
 
 
 
-  template <>
-  double
-  cell_measure<1>(
-    const std::vector<Point<1>> &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<Point<2>> &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 := <diff(x_real,xi)  | diff(y_real,xi) | diff(z_real,xi)>;
-      Jeta := <diff(x_real,eta)| diff(y_real,eta)| diff(z_real,eta)>;
-      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<Point<3>> &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 <int dim>
   Vector<double>
   compute_aspect_ratio_of_cells(const Triangulation<dim> &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 <int dim>
   void
   rotate(const double           angle,
@@ -1141,19 +873,6 @@ namespace GridTools
   } // namespace
 
 
-
-  // Implementation for 1D only
-  template <>
-  void
-  laplace_transform(const std::map<unsigned int, Point<1>> &,
-                    Triangulation<1> &,
-                    const Function<1> *,
-                    const bool)
-  {
-    Assert(false, ExcNotImplemented());
-  }
-
-
   // Implementation for dimensions except 1
   template <int dim>
   void
@@ -1162,6 +881,9 @@ namespace GridTools
                     const Function<dim> *                     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<dim> q1(1);
@@ -3748,16 +3470,6 @@ namespace GridTools
 
 
 
-      void
-      fix_up_faces(const dealii::Triangulation<1, 1>::cell_iterator &,
-                   std::integral_constant<int, 1>,
-                   std::integral_constant<int, 1>)
-      {
-        // 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 <int dim, int spacedim>
       void
@@ -3810,6 +3522,8 @@ namespace GridTools
       &distorted_cells,
     Triangulation<dim, spacedim> & /*triangulation*/)
   {
+    static_assert(dim != 1 && spacedim != 1,
+                  "This function is only valid when dim != 1 or spacedim != 1.");
     typename Triangulation<dim, spacedim>::DistortedCellList unfixable_subset;
 
     // loop over all cells that we have to fix up
index ff06ead6babdb8422067ea21528d45c79fa96bc8..824c610b40783c8a751e22f86b27b19785994811 100644 (file)
@@ -357,8 +357,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
         const Function<deal_II_dimension> *                     coefficient,
         const bool);
 
-#    endif
-
       template Triangulation<deal_II_dimension,
                              deal_II_space_dimension>::DistortedCellList
       fix_up_distorted_child_cells(
@@ -367,6 +365,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
           &distorted_cells,
         Triangulation<deal_II_dimension, deal_II_space_dimension>
           &triangulation);
+#    endif
 #  endif
     \}
 #endif
diff --git a/source/grid/grid_tools_nontemplates.cc b/source/grid/grid_tools_nontemplates.cc
new file mode 100644 (file)
index 0000000..706dbc4
--- /dev/null
@@ -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 <deal.II/base/point.h>
+
+#include <deal.II/grid/grid_tools.h>
+
+#include <vector>
+
+// 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<Point<1>> &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<Point<2>> &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 := <diff(x_real,xi)  | diff(y_real,xi) | diff(z_real,xi)>;
+      Jeta := <diff(x_real,eta)| diff(y_real,eta)| diff(z_real,eta)>;
+      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<Point<3>> &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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.