doxygen/deal.II/grid__tools__nontemplates_8cc_source.html

// ---------------------------------------------------------------------

//

// Copyright (C) 2001 - 2022 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 <deal.II/physics/transformations.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 ArrayView<const unsigned int> &vertex_indices)

  {

    AssertDimension(vertex_indices.size(), 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 ArrayView<const unsigned int> &vertex_indices)

  {

    if (vertex_indices.size() == 3) // triangle

      {

        const double x[3] = {all_vertices[vertex_indices[0]](0),

                             all_vertices[vertex_indices[1]](0),

                             all_vertices[vertex_indices[2]](0)};


        const double y[3] = {all_vertices[vertex_indices[0]](1),

                             all_vertices[vertex_indices[1]](1),

                             all_vertices[vertex_indices[2]](1)};


        return 0.5 *

               ((x[0] - x[2]) * (y[1] - y[0]) - (x[1] - x[0]) * (y[0] - y[2]));

      }


    AssertDimension(vertex_indices.size(), 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 ArrayView<const unsigned int> &vertex_indices)

  {

    if (vertex_indices.size() == 4) // tetrahedron

      {

        const auto &a = all_vertices[vertex_indices[0]];

        const auto &b = all_vertices[vertex_indices[1]];

        const auto &c = all_vertices[vertex_indices[2]];

        const auto &d = all_vertices[vertex_indices[3]];


        return (1.0 / 6.0) * (d - a) * cross_product_3d(b - a, c - a);

      }

    else if (vertex_indices.size() == 5) // pyramid

      {

        // This remarkably simple formula comes from Equation 4 of

        // "Calculation of the volume of a general hexahedron for flow

        // predictions", Davies and Salmond, AIAA Journal vol. 23 no. 6.

        const auto &x0 = all_vertices[vertex_indices[0]];

        const auto &x1 = all_vertices[vertex_indices[1]];

        const auto &x2 = all_vertices[vertex_indices[2]];

        const auto &x3 = all_vertices[vertex_indices[3]];

        const auto &x4 = all_vertices[vertex_indices[4]];


        const auto v01 = x1 - x0;

        const auto v02 = x2 - x0;

        const auto v03 = x3 - x0;

        const auto v04 = x4 - x0;

        const auto v21 = x2 - x1;


        // doing high - low consistently puts us off by -1 from the original

        // paper in the first term

        return -v04 * cross_product_3d(v21, v03) / 6.0 +

               v03 * cross_product_3d(v01, v02) / 12.0;

      }

    else if (vertex_indices.size() == 6) // wedge

      {

        /* Script used to generate volume code:


           #!/usr/bin/env python

           # coding: utf-8

           import sympy as sp

           from sympy.simplify.cse_main import cse

           n_vertices = 6

           xs = list(sp.symbols(" ".join(["x{}".format(i)

                                          for i in range(n_vertices)])))

           ys = list(sp.symbols(" ".join(["y{}".format(i)

                                          for i in range(n_vertices)])))

           zs = list(sp.symbols(" ".join(["z{}".format(i)

                                          for i in range(n_vertices)])))

           xi, eta, zeta = sp.symbols("xi eta zeta")

           tphi = [(1 - xi - eta)*(1 - zeta),

                   (xi)*(1 - zeta),

                   (eta)*(1 - zeta),

                   (1 - xi - eta)*(zeta),

                   (xi)*(zeta),

                   (eta)*(zeta)]

           x_real = sum(xs[i]*tphi[i] for i in range(n_vertices))

           y_real = sum(ys[i]*tphi[i] for i in range(n_vertices))

           z_real = sum(zs[i]*tphi[i] for i in range(n_vertices))

           J = sp.Matrix([[var.diff(v) for v in [xi, eta, zeta]]

                           for var in [x_real, y_real, z_real]])

           detJ = J.det()

           detJ2 = detJ.expand().collect(zeta).collect(eta).collect(xi)

           for x in xs:

               detJ2 = detJ2.collect(x)

           for y in ys:

               detJ2 = detJ2.collect(y)

           for z in zs:

               detJ2 = detJ2.collect(z)

           measure = sp.integrate(sp.integrate(

                                  sp.integrate(detJ2, (eta, 0, 1 - xi)),

                                               (xi, 0, 1)), (zeta, 0, 1))

           measure2 = measure

           for vs in [xs, ys, zs]:

               for v in vs:

                   measure2 = measure2.collect(v)


           pairs, expression = cse(measure2)

           for vertex_no in range(n_vertices):

               for (coordinate, index) in [('x', 0), ('y', 1), ('z', 2)]:

                   print(

                   "const double {}{} = all_vertices[vertex_indices[{}]][{}];"

                   .format(coordinate, vertex_no, vertex_no, index))


           for pair in pairs:

               print("const double " + sp.ccode(pair[0]) + " = "

               + sp.ccode(pair[1]) + ";")

           print("const double result = " + sp.ccode(expression[0]) + ";")

           print("return result;")

         */

        const double x0  = all_vertices[vertex_indices[0]][0];

        const double y0  = all_vertices[vertex_indices[0]][1];

        const double z0  = all_vertices[vertex_indices[0]][2];

        const double x1  = all_vertices[vertex_indices[1]][0];

        const double y1  = all_vertices[vertex_indices[1]][1];

        const double z1  = all_vertices[vertex_indices[1]][2];

        const double x2  = all_vertices[vertex_indices[2]][0];

        const double y2  = all_vertices[vertex_indices[2]][1];

        const double z2  = all_vertices[vertex_indices[2]][2];

        const double x3  = all_vertices[vertex_indices[3]][0];

        const double y3  = all_vertices[vertex_indices[3]][1];

        const double z3  = all_vertices[vertex_indices[3]][2];

        const double x4  = all_vertices[vertex_indices[4]][0];

        const double y4  = all_vertices[vertex_indices[4]][1];

        const double z4  = all_vertices[vertex_indices[4]][2];

        const double x5  = all_vertices[vertex_indices[5]][0];

        const double y5  = all_vertices[vertex_indices[5]][1];

        const double z5  = all_vertices[vertex_indices[5]][2];

        const double x6  = (1.0 / 12.0) * z1;

        const double x7  = -x6;

        const double x8  = (1.0 / 12.0) * z3;

        const double x9  = x7 + x8;

        const double x10 = (1.0 / 12.0) * z2;

        const double x11 = -x8;

        const double x12 = x10 + x11;

        const double x13 = (1.0 / 6.0) * z2;

        const double x14 = (1.0 / 12.0) * z4;

        const double x15 = (1.0 / 6.0) * z1;

        const double x16 = (1.0 / 12.0) * z5;

        const double x17 = -x16;

        const double x18 = x16 + x7;

        const double x19 = -x14;

        const double x20 = x10 + x19;

        const double x21 = (1.0 / 12.0) * z0;

        const double x22 = x19 + x21;

        const double x23 = -x10;

        const double x24 = x14 + x23;

        const double x25 = (1.0 / 6.0) * z0;

        const double x26 = x17 + x21;

        const double x27 = x23 + x8;

        const double x28 = -x21;

        const double x29 = x16 + x28;

        const double x30 = x17 + x6;

        const double x31 = x14 + x28;

        const double x32 = x11 + x6;

        const double x33 = (1.0 / 6.0) * z5;

        const double x34 = (1.0 / 6.0) * z4;

        const double x35 = (1.0 / 6.0) * z3;

        const double result =

          x0 * (x12 * y5 + x9 * y4 + y1 * (-x13 + x14 + x8) +

                y2 * (x11 + x15 + x17) + y3 * (x18 + x20)) +

          x1 * (x22 * y3 + x24 * y5 + y0 * (x11 + x13 + x19) +

                y2 * (x14 + x16 - x25) + y4 * (x26 + x27)) +

          x2 * (x29 * y3 + x30 * y4 + y0 * (-x15 + x16 + x8) +

                y1 * (x17 + x19 + x25) + y5 * (x31 + x32)) +

          x3 * (x26 * y2 + x31 * y1 + y0 * (x24 + x30) + y4 * (x28 + x33 + x7) +

                y5 * (x10 + x21 - x34)) +

          x4 * (x18 * y2 + x32 * y0 + y1 * (x12 + x29) + y3 * (x21 - x33 + x6) +

                y5 * (x23 + x35 + x7)) +

          x5 * (x20 * y1 + x27 * y0 + y2 * (x22 + x9) + y3 * (x23 + x28 + x34) +

                y4 * (x10 - x35 + x6));

        return result;

      }


    AssertDimension(vertex_indices.size(), 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)};


    /*

      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)

        : rotation_matrix(

            Physics::Transformations::Rotations::rotation_matrix_2d(angle))

      {}

      Point<2>

      operator()(const Point<2> &p) const

      {

        return static_cast<Point<2>>(rotation_matrix * p);

      }


    private:

      const Tensor<2, 2, double> rotation_matrix;

    };

  } // 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

point.h

ArrayView
Definition: array_view.h:82

Point
Definition: point.h:111

Tensor
Definition: tensor.h:503

Triangulation
Definition: tria.h:1122

DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442

DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443

vertex_indices
unsigned int vertex_indices[2]
Definition: grid_tools.cc:1293

grid_tools.h

rotation_matrix
const Tensor< 2, 2, double > rotation_matrix
Definition: grid_tools_nontemplates.cc:440

AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667

StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)

AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583

GridTools
Definition: grid_tools.h:132

GridTools::transform
void transform(const Transformation &transformation, Triangulation< dim, spacedim > &triangulation)

GridTools::cell_measure< 2 >
double cell_measure< 2 >(const std::vector< Point< 2 > > &all_vertices, const ArrayView< const unsigned int > &vertex_indices)
Definition: grid_tools_nontemplates.cc:47

GridTools::cell_measure< 1 >
double cell_measure< 1 >(const std::vector< Point< 1 > > &all_vertices, const ArrayView< const unsigned int > &vertex_indices)
Definition: grid_tools_nontemplates.cc:34

GridTools::rotate
void rotate(const double angle, Triangulation< dim > &triangulation)

GridTools::cell_measure< 3 >
double cell_measure< 3 >(const std::vector< Point< 3 > > &all_vertices, const ArrayView< const unsigned int > &vertex_indices)
Definition: grid_tools_nontemplates.cc:121

Physics
Definition: physics.h:29

triangulation
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69

GeometryInfo
Definition: geometry_info.h:1969

cross_product_3d
constexpr Tensor< 1, dim, typename ProductType< Number1, Number2 >::type > cross_product_3d(const Tensor< 1, dim, Number1 > &src1, const Tensor< 1, dim, Number2 > &src2)
Definition: tensor.h:2660

transformations.h