+ // 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,
grid_generator.cc
grid_tools.cc
grid_tools_cache.cc
+ grid_tools_nontemplates.cc
grid_tools_dof_handlers.cc
tria.cc
)
- 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,
return SubCellData();
}
};
-
-
} // 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;
- };
-
// Transformation to rotate around one of the cartesian axes.
class Rotate3d
{
}
- 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,
} // 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
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);
- 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
&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
const Function<deal_II_dimension> * coefficient,
const bool);
-# endif
-
template Triangulation<deal_II_dimension,
deal_II_space_dimension>::DistortedCellList
fix_up_distorted_child_cells(
&distorted_cells,
Triangulation<deal_II_dimension, deal_II_space_dimension>
&triangulation);
+# endif
# endif
\}
#endif
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
projects / dealii.git / commitdiff