lagrangepoints.impl.hh

Go to the documentation of this file.

#pragma once
 
#include <cassert>
#include <array>
 
#include <dune/common/exceptions.hh>
#include <dune/common/version.hh>
#include <dune/geometry/dimension.hh>
#include <dune/geometry/type.hh>
#include <dune/localfunctions/lagrange/equidistantpoints.hh>
 
namespace Dune::Vtk {
namespace Impl {
 
template <class K>
  template <class Points>
void LagrangePointSetBuilder<K,0>::operator() (GeometryType gt, int /*order*/, Points& points) const
{
  assert(gt.isVertex());
  points.push_back(LP{Vec{},Key{0,0,0}});
}
 
 
template <class K>
  template <class Points>
void LagrangePointSetBuilder<K,1>::operator() (GeometryType gt, int order, Points& points) const
{
  assert(gt.isLine());
 
  if (order == 0) {
    auto refElem = Dune::referenceElement<K,1>(gt);
    points.push_back(LP{refElem.position(0,0), Key{0,0,0}});
    return;
  }
 
  // Vertex nodes
  points.push_back(LP{Vec{0.0}, Key{0,dim,0}});
  points.push_back(LP{Vec{1.0}, Key{1,dim,0}});
 
  if (order > 1) {
    // Inner nodes
    Vec p{0.0};
    for (int i = 0; i < order-1; i++)
    {
      p[0] += 1.0 / order;
      points.push_back(LP{p,Key{0,dim-1,(unsigned int)(i)}});
    }
  }
}
 
 
template <class K>
  template <class Points>
void LagrangePointSetBuilder<K,2>::operator() (GeometryType gt, int order, Points& points) const
{
  if (order == 0) {
    auto refElem = Dune::referenceElement<K,2>(gt);
    points.push_back(LP{refElem.position(0,0), Key{0,0,0}});
    return;
  }
 
  if (gt.isTriangle())
    buildTriangle(order, points);
  else if (gt.isQuadrilateral())
    buildQuad(order, order, points);
  else {
    DUNE_THROW(Dune::NotImplemented,
      "Lagrange points not yet implemented for this GeometryType.");
  }
 
  [[maybe_unused]] std::size_t nPoints = numLagrangePoints(gt, order);
  assert(points.size() == nPoints);
}
 
 
// Construct the point set in a triangle element.
// Loop from the outside to the inside
template <class K>
  template <class Points>
void LagrangePointSetBuilder<K,2>::buildTriangle (int order, Points& points) const
{
  std::size_t nPoints = numLagrangePoints(GeometryTypes::simplex(2), order);
  points.reserve(nPoints);
 
  const int nVertexDOFs = 3;
  const int nEdgeDOFs = 3 * (order-1);
 
  static const unsigned int vertexPerm[3] = {0,1,2};
  static const unsigned int edgePerm[3]   = {0,2,1};
 
  auto calcKey = [&](int p) -> Key
  {
    if (p < nVertexDOFs) {
      return Key{vertexPerm[p], dim, 0};
    }
    else if (p < nVertexDOFs+nEdgeDOFs) {
      unsigned int entityIndex = (p - nVertexDOFs) / (order-1);
      unsigned int index = (p - nVertexDOFs) % (order-1);
      return Key{edgePerm[entityIndex], dim-1, index};
    }
    else {
      unsigned int index = p - (nVertexDOFs + nEdgeDOFs);
      return Key{0, dim-2, index};
    }
  };
 
  std::array<int,3> bindex;
 
  K order_d = K(order);
  for (std::size_t p = 0; p < nPoints; ++p) {
    barycentricIndex(p, bindex, order);
    Vec point{bindex[0] / order_d, bindex[1] / order_d};
    points.push_back(LP{point, calcKey(p)});
  }
}
 
 
// "Barycentric index" is a triplet of integers, each running from 0 to
// <Order>. It is the index of a point on the triangle in barycentric
// coordinates.
template <class K>
void LagrangePointSetBuilder<K,2>::barycentricIndex (int index, std::array<int,3>& bindex, int order)
{
  int max = order;
  int min = 0;
 
  // scope into the correct triangle
  while (index != 0 && index >= 3 * order)
  {
    index -= 3 * order;
    max -= 2;
    min++;
    order -= 3;
  }
 
  // vertex DOFs
  if (index < 3)
  {
    bindex[index] = bindex[(index + 1) % 3] = min;
    bindex[(index + 2) % 3] = max;
  }
  // edge DOFs
  else
  {
    index -= 3;
    int dim = index / (order - 1);
    int offset = (index - dim * (order - 1));
    bindex[(dim + 1) % 3] = min;
    bindex[(dim + 2) % 3] = (max - 1) - offset;
    bindex[dim] = (min + 1) + offset;
  }
}
 
 
// Construct the point set in the quad element
// 1. build equispaced points with index tuple (i,j)
// 2. map index tuple to DOF index and LocalKey
template <class K>
  template <class Points>
void LagrangePointSetBuilder<K,2>::buildQuad(int orderx, int ordery, Points& points) const
{
  if (orderx == 0 && ordery == 0) {
    points.push_back(LP{Vec{0.5, 0.5}});
    return;
  }
  else if (orderx == 0) {
    std::vector< LagrangePoint<K,1> > linePoints;
    LagrangePointSetBuilder<K,1>{}(GeometryTypes::line, ordery, linePoints);
    points.reserve(linePoints.size());
    for (auto& p : linePoints)
      points.push_back(LP{Vec{0.5, p.point()[0]}});
    return;
  }
  else if (ordery == 0) {
    std::vector< LagrangePoint<K,1> > linePoints;
    LagrangePointSetBuilder<K,1>{}(GeometryTypes::line, orderx, linePoints);
    points.reserve(linePoints.size());
    for (auto& p : linePoints)
      points.push_back(LP{Vec{p.point()[0], 0.5}});
    return;
  }
 
  std::size_t nPoints = numLagrangePoints(GeometryTypes::line, orderx);
  nPoints *= numLagrangePoints(GeometryTypes::line, ordery);
  points.resize(nPoints);
 
  assert(orderx > 0 && ordery > 0);
  std::array<int,2> orders{orderx, ordery};
  std::array<Vec,4> nodes{Vec{0., 0.}, Vec{1., 0.}, Vec{1., 1.}, Vec{0., 1.}};
 
  for (int n = 0; n <= orders[1]; ++n) {
    for (int m = 0; m <= orders[0]; ++m) {
      // int idx = pointIndexFromIJ(m,n,orders);
 
      const K r = K(m) / orders[0];
      const K s = K(n) / orders[1];
      Vec p = K(1.0-r) * (nodes[3] * s + nodes[0] * K(1.0 - s)) +
              r *        (nodes[2] * s + nodes[1] * K(1.0 - s));
 
      auto [idx,key] = calcQuadKey(m,n,orders);
      points[idx] = LP{p, key};
      // points[idx] = LP{p, calcQuadKey(n,m,orders)};
    }
  }
}
 
 
// Obtain the VTK DOF index of the node (i,j) in the quad element
// and construct a LocalKey
template <class K>
std::pair<int,typename LagrangePointSetBuilder<K,2>::Key>
LagrangePointSetBuilder<K,2>::calcQuadKey (int i, int j, std::array<int,2> order)
{
  const bool ibdy = (i == 0 || i == order[0]);
  const bool jbdy = (j == 0 || j == order[1]);
  const int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0); // How many boundaries do we lie on at once?
 
  int dof = 0;
  unsigned int entityIndex = 0;
  unsigned int index = 0;
 
  if (nbdy == 2) // Vertex DOF
  {
    dof = (i ? (j ? 2 : 1) : (j ? 3 : 0));
    entityIndex = (j ? (i ? 3 : 2) : (i ? 1 : 0));
    return std::make_pair(dof,Key{entityIndex, dim, 0});
  }
 
  int offset = 4;
  if (nbdy == 1) // Edge DOF
  {
    if (!ibdy) {
      dof = (i - 1) + (j ? order[0]-1 + order[1]-1 : 0) + offset;
      entityIndex = j ? 3 : 2;
      index = i-1;
    }
    else if (!jbdy) {
      dof = (j - 1) + (i ? order[0]-1 : 2 * (order[0]-1) + order[1]-1) + offset;
      entityIndex = i ? 1 : 0;
      index = j-1;
    }
    return std::make_pair(dof, Key{entityIndex, dim-1, index});
  }
 
  offset += 2 * (order[0]-1 + order[1]-1);
 
  // nbdy == 0: Face DOF
  dof = offset + (i - 1) + (order[0]-1) * ((j - 1));
  Key innerKey = LagrangePointSetBuilder<K,2>::calcQuadKey(i-1,j-1,{order[0]-2, order[1]-2}).second;
  return std::make_pair(dof, Key{0, dim-2, innerKey.index()});
}
 
 
// Lagrange points on 3d geometries
template <class K>
  template <class Points>
void LagrangePointSetBuilder<K,3>::operator() (GeometryType gt, unsigned int order, Points& points) const
{
  if (order == 0) {
    auto refElem = Dune::referenceElement<K,3>(gt);
    points.push_back(LP{refElem.position(0,0), Key{0,0,0}});
    return;
  }
 
  if (gt.isTetrahedron())
    buildTetra(order, points);
  else if (gt.isHexahedron())
    buildHex(order, points);
  else if (gt.isPrism())
    buildWedge(order, order, points);
  else {
    DUNE_THROW(Dune::NotImplemented,
      "Lagrange points not yet implemented for this GeometryType.");
  }
 
  [[maybe_unused]] std::size_t nPoints = numLagrangePoints(gt, order);
  assert(points.size() == nPoints);
}
 
 
// Construct the point set in the tetrahedron element
// 1. construct barycentric (index) coordinates
// 2. obtains the DOF index, LocalKey and actual coordinate from barycentric index
template <class K>
  template <class Points>
void LagrangePointSetBuilder<K,3>::buildTetra (int order, Points& points) const
{
  std::size_t nPoints = numLagrangePoints(GeometryTypes::simplex(3), order);
  points.reserve(nPoints);
 
  const int nVertexDOFs = 4;
  const int nEdgeDOFs = 6 * (order-1);
  const int nFaceDOFs = 4 * (order-1)*(order-2)/2;
 
  static const unsigned int vertexPerm[4] = {0,1,2,3};
  static const unsigned int edgePerm[6]   = {0,2,1,3,4,5};
  static const unsigned int facePerm[4]   = {1,2,0,3};
 
  auto calcKey = [&](int p) -> Key
  {
    if (p < nVertexDOFs) {
      return Key{vertexPerm[p], dim, 0};
    }
    else if (p < nVertexDOFs+nEdgeDOFs) {
      unsigned int entityIndex = (p - nVertexDOFs) / (order-1);
      unsigned int index = (p - nVertexDOFs) % (order-1);
      return Key{edgePerm[entityIndex], dim-1, index};
    }
    else if (p < nVertexDOFs+nEdgeDOFs+nFaceDOFs) {
      unsigned int index = (p - (nVertexDOFs + nEdgeDOFs)) % ((order-1)*(order-2)/2);
      unsigned int entityIndex = (p - (nVertexDOFs + nEdgeDOFs)) / ((order-1)*(order-2)/2);
      return Key{facePerm[entityIndex], dim-2, index};
    }
    else {
      unsigned int index = p - (nVertexDOFs + nEdgeDOFs + nFaceDOFs);
      return Key{0, dim-3, index};
    }
  };
 
  std::array<int,4> bindex;
 
  K order_d = K(order);
  for (std::size_t p = 0; p < nPoints; ++p) {
    barycentricIndex(p, bindex, order);
    Vec point{bindex[0] / order_d, bindex[1] / order_d, bindex[2] / order_d};
    points.push_back(LP{point, calcKey(p)});
  }
}
 
 
// "Barycentric index" is a set of 4 integers, each running from 0 to
// <Order>. It is the index of a point in the tetrahedron in barycentric
// coordinates.
template <class K>
void LagrangePointSetBuilder<K,3>::barycentricIndex (int p, std::array<int,4>& bindex, int order)
{
  const int nVertexDOFs = 4;
  const int nEdgeDOFs = 6 * (order-1);
 
  static const int edgeVertices[6][2]   = {{0,1}, {1,2}, {2,0}, {0,3}, {1,3}, {2,3}};
  static const int linearVertices[4][4] = {{0,0,0,1}, {1,0,0,0}, {0,1,0,0}, {0,0,1,0}};
  static const int vertexMaxCoords[4]   = {3,0,1,2};
  static const int faceBCoords[4][3]    = {{0,2,3}, {2,0,1}, {2,1,3}, {1,0,3}};
  static const int faceMinCoord[4]      = {1,3,0,2};
 
  int max = order;
  int min = 0;
 
  // scope into the correct tetra
  while (p >= 2 * (order * order + 1) && p != 0 && order > 3)
  {
    p -= 2 * (order * order + 1);
    max -= 3;
    min++;
    order -= 4;
  }
 
  // vertex DOFs
  if (p < nVertexDOFs)
  {
    for (int coord = 0; coord < 4; ++coord)
      bindex[coord] = (coord == vertexMaxCoords[p] ? max : min);
  }
  // edge DOFs
  else if (p < nVertexDOFs+nEdgeDOFs)
  {
    int edgeId = (p - nVertexDOFs) / (order-1);
    int vertexId = (p - nVertexDOFs) % (order-1);
    for (int coord = 0; coord < 4; ++coord)
    {
      bindex[coord] = min +
        (linearVertices[edgeVertices[edgeId][0]][coord] * (max - min - 1 - vertexId) +
          linearVertices[edgeVertices[edgeId][1]][coord] * (1 + vertexId));
    }
  }
  // face DOFs
  else
  {
    int faceId = (p - (nVertexDOFs+nEdgeDOFs)) / ((order-2)*(order-1)/2);
    int vertexId = (p - (nVertexDOFs+nEdgeDOFs)) % ((order-2)*(order-1)/2);
 
    std::array<int,3> projectedBIndex;
    if (order == 3)
      projectedBIndex[0] = projectedBIndex[1] = projectedBIndex[2] = 0;
    else
      LagrangePointSetBuilder<K,2>::barycentricIndex(vertexId, projectedBIndex, order-3);
 
    for (int i = 0; i < 3; i++)
      bindex[faceBCoords[faceId][i]] = (min + 1 + projectedBIndex[i]);
 
    bindex[faceMinCoord[faceId]] = min;
  }
}
 
 
// Construct the point set in the hexahedral element
// 1. build equispaced points with index tuple (i,j,k)
// 2. map index tuple to DOF index and LocalKey
template <class K>
  template <class Points>
void LagrangePointSetBuilder<K,3>::buildHex (int order, Points& points) const
{
  std::size_t nPoints = numLagrangePoints(GeometryTypes::cube(3), order);
  points.resize(nPoints);
 
  std::array<int,3> orders{order, order, order};
  std::array<Vec,8> nodes{Vec{0., 0., 0.}, Vec{1., 0., 0.}, Vec{1., 1., 0.}, Vec{0., 1., 0.},
                          Vec{0., 0., 1.}, Vec{1., 0., 1.}, Vec{1., 1., 1.}, Vec{0., 1., 1.}};
 
  for (int p = 0; p <= orders[2]; ++p) {
    for (int n = 0; n <= orders[1]; ++n) {
      for (int m = 0; m <= orders[0]; ++m) {
        const K r = K(m) / orders[0];
        const K s = K(n) / orders[1];
        const K t = K(p) / orders[2];
        Vec point = K(1.0-r) * ((nodes[3] * K(1.0-t) + nodes[7] * t) * s + (nodes[0] * K(1.0-t) + nodes[4] * t) * K(1.0-s)) +
                    r *        ((nodes[2] * K(1.0-t) + nodes[6] * t) * s + (nodes[1] * K(1.0-t) + nodes[5] * t) * K(1.0-s));
 
        auto [idx,key] = calcHexKey(m,n,p,orders);
        points[idx] = LP{point, key};
      }
    }
  }
}
 
 
// Obtain the VTK DOF index of the node (i,j,k) in the hexahedral element
template <class K>
std::pair<int,typename LagrangePointSetBuilder<K,3>::Key>
LagrangePointSetBuilder<K,3>::calcHexKey (int i, int j, int k, std::array<int,3> order)
{
  const bool ibdy = (i == 0 || i == order[0]);
  const bool jbdy = (j == 0 || j == order[1]);
  const bool kbdy = (k == 0 || k == order[2]);
  const int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0) + (kbdy ? 1 : 0); // How many boundaries do we lie on at once?
 
  int dof = 0;
  unsigned int entityIndex = 0;
  unsigned int index = 0;
 
  if (nbdy == 3) // Vertex DOF
  {
    dof = (i ? (j ? 2 : 1) : (j ? 3 : 0)) + (k ? 4 : 0);
    entityIndex = (i ? 1 : 0) + (j ? 2 : 0) + (k ? 4 : 0);
    return std::make_pair(dof, Key{entityIndex, dim, 0});
  }
 
  int offset = 8;
  if (nbdy == 2) // Edge DOF
  {
    entityIndex = (k ? 8 : 4);
    if (!ibdy)
    { // On i axis
      dof = (i - 1) + (j ? order[0]-1 + order[1]-1 : 0) + (k ? 2 * (order[0]-1 + order[1]-1) : 0) + offset;
      index = i;
      entityIndex += (i ? 1 : 0);
    }
    else if (!jbdy)
    { // On j axis
      dof = (j - 1) + (i ? order[0]-1 : 2 * (order[0]-1) + order[1]-1) + (k ? 2 * (order[0]-1 + order[1]-1) : 0) + offset;
      index = j;
      entityIndex += (j ? 3 : 2);
    }
    else
    { // !kbdy, On k axis
      offset += 4 * (order[0]-1) + 4 * (order[1]-1);
      dof = (k - 1) + (order[2]-1) * (i ? (j ? 3 : 1) : (j ? 2 : 0)) + offset;
      index = k;
      entityIndex = (i ? 1 : 0) + (j ? 2 : 0);
    }
    return std::make_pair(dof, Key{entityIndex, dim-1, index});
  }
 
  offset += 4 * (order[0]-1 + order[1]-1 + order[2]-1);
  if (nbdy == 1) // Face DOF
  {
    Key faceKey;
    if (ibdy) // On i-normal face
    {
      dof = (j - 1) + ((order[1]-1) * (k - 1)) + (i ? (order[1]-1) * (order[2]-1) : 0) + offset;
      entityIndex = (i ? 1 : 0);
      faceKey = LagrangePointSetBuilder<K,2>::calcQuadKey(j-1,k-1,{order[1]-2, order[2]-2}).second;
    }
    else {
      offset += 2 * (order[1] - 1) * (order[2] - 1);
      if (jbdy) // On j-normal face
      {
        dof = (i - 1) + ((order[0]-1) * (k - 1)) + (j ? (order[2]-1) * (order[0]-1) : 0) + offset;
        entityIndex = (j ? 3 : 2);
        faceKey = LagrangePointSetBuilder<K,2>::calcQuadKey(i-1,k-1,{order[0]-2, order[2]-2}).second;
      }
      else
      { // kbdy, On k-normal face
        offset += 2 * (order[2]-1) * (order[0]-1);
        dof = (i - 1) + ((order[0]-1) * (j - 1)) + (k ? (order[0]-1) * (order[1]-1) : 0) + offset;
        entityIndex = (k ? 5 : 4);
        faceKey = LagrangePointSetBuilder<K,2>::calcQuadKey(i-1,j-1,{order[0]-2, order[1]-2}).second;
      }
    }
    return std::make_pair(dof, Key{entityIndex, dim-2, faceKey.index()});
  }
 
  // nbdy == 0: Body DOF
  offset += 2 * ((order[1]-1) * (order[2]-1) + (order[2]-1) * (order[0]-1) + (order[0]-1) * (order[1]-1));
  dof = offset + (i - 1) + (order[0]-1) * ((j - 1) + (order[1]-1) * ((k - 1)));
  Key innerKey = LagrangePointSetBuilder<K,3>::calcHexKey(i-1,j-1,k-1,{order[0]-2, order[1]-2, order[2]-2}).second;
  return std::make_pair(dof, Key{0, dim-3, innerKey.index()});
}
 
template <class K>
  template <class Points>
void LagrangePointSetBuilder<K,3>::buildWedge (int order_tri, int order_quad, Points& points) const
{
  auto refElem = Dune::referenceElement<K,3>(GeometryTypes::prism);
 
  // flip the y-axis
  auto flip = [](auto x) { x[1] = 1-x[1]; return x; };
 
  // transform the coordinates by shifting and scaling
  auto shrink_linexy = [order_tri](auto x) { x *= (1-K(2)/order_tri); x += decltype(x)(K(1)/order_tri); return x; };
  auto shrink_linez = [order_quad](auto x) { x *= (1-K(2)/order_quad); x += decltype(x)(K(1)/order_quad); return x; };
  auto shrink_tri = [order_tri](auto x) { x *= (1-K(3)/order_tri); x += decltype(x)(K(1)/order_tri); return x; };
  auto shrink_quad = [order_tri,order_quad](auto x) {
    x[0] *= (1-K(2)/order_tri);
    x[1] *= (1-K(2)/order_quad);
    x += decltype(x){K(1)/order_tri, K(1)/order_quad};
    return x;
  };
  auto shrink_prism = [order_tri,order_quad](auto x) {
    x[0] *= (1-K(3)/order_tri);
    x[1] *= (1-K(3)/order_tri);
    x[2] *= (1-K(2)/order_quad);
    x += decltype(x){K(1)/order_tri, K(1)/order_tri, K(1)/order_quad};
    return x;
  };
 
  // embed the point set into a sub-area of a subentity of the prism element
  auto embed = [&](auto const& entityPoints, auto const& perm, auto codim, auto shrink) {
    for (int p : perm) {
      unsigned int s = std::abs(p);
      constexpr unsigned int c = codim.value;
      auto geo = refElem.template geometry<c>(s);
 
      unsigned int i = 0;
      for (auto& lp : entityPoints) {
        if (p < 0)
          points.push_back(LP{Vec{flip(geo.global(shrink(lp.point())))}, Key{s,c,i}});
        else
          points.push_back(LP{Vec{geo.global(shrink(lp.point()))}, Key{s,c,i}});
        i++;
      }
    }
  };
 
  std::size_t nPoints = numLagrangePoints(GeometryTypes::simplex(2), order_tri);
  nPoints *= numLagrangePoints(GeometryTypes::line, order_quad);
  points.reserve(nPoints);
 
  if (order_tri >= 0 && order_quad >= 0) {
    // 1. vertex and inner DOFs
    std::vector< LagrangePoint<K,2> > triPoints;
    LagrangePointSetBuilder<K,2>{}(GeometryTypes::simplex(2), std::min(order_tri,1), triPoints);
 
    std::vector< LagrangePoint<K,1> > linePoints;
    LagrangePointSetBuilder<K,1>{}(GeometryTypes::line, std::min(order_quad,1), linePoints);
 
    unsigned int i = 0;
    for (auto& p0 : linePoints)
      for (auto& p1 : triPoints)
        points.push_back(LP{Vec{p1.point()[0], p1.point()[1], p0.point()[0]},Key{i++,3,0}});
  }
 
  if (order_tri > 1) {
    // 2. edge DOFs on horizontal
    std::vector< LagrangePoint<K,1> > linePoints;
    LexicographicPointSetBuilder<K,1>::buildLine(order_tri-2, linePoints);
    embed(linePoints, std::array{3,5,-4,6,8,-7}, Codim<2>{}, shrink_linexy);
  }
 
  if (order_quad > 1) {
    // 2. edge DOFs on vertical edges
    std::vector< LagrangePoint<K,1> > linePoints;
    LexicographicPointSetBuilder<K,1>::buildLine(order_quad-2, linePoints);
    if (order_tri == 0)
      embed(linePoints, std::array{0}, Codim<2>{}, shrink_linez);
    else
      embed(linePoints, std::array{0,1,2}, Codim<2>{}, shrink_linez);
  }
 
  if (order_tri > 2) {
    // 3. triangular face DOFs
    std::vector< LagrangePoint<K,2> > triPoints;
    LagrangePointSetBuilder<K,2>{}(GeometryTypes::simplex(2), order_tri-3, triPoints);
    embed(triPoints, std::array{3,4}, Codim<1>{}, shrink_tri);
  }
 
  if (order_tri > 1 && order_quad > 1) {
    // 4. quad face DOFs
    std::vector< LagrangePoint<K,2> > quadPoints;
    LexicographicPointSetBuilder<K,2>::buildQuad(order_tri-2, order_quad-2, quadPoints);
    embed(quadPoints, std::array{0,2,-1}, Codim<1>{}, shrink_quad);
  }
 
  if (order_tri >= 3 && order_quad >= 2) {
    // 5. volume DOFs
    std::vector< LagrangePoint<K,3> > prismPoints;
    buildWedge(order_tri-3, order_quad-2, prismPoints);
    embed(prismPoints, std::array{0}, Codim<0>{}, shrink_prism);
  }
}
 
 
template <class K, unsigned int dim>
  template <class Points>
void LexicographicPointSetBuilder<K,dim>::buildLine (int a, Points& points)
{
  static_assert(dim == 1);
  using Point = typename Points::value_type;
  using Vec = typename Point::Vector;
  using Key = Dune::LocalKey;
 
  auto coord = [a](int i) { return a == 0 ? K(0.5) : K(i)/a; };
 
  for (int i = 0; i < a+1; ++i) {
    auto key = (i == 0 || i == a)
      ? Key(0u + (i > 0), 1, 0)   // vertex
      : Key(0, 0, i-1);           // edge
    points.push_back(Point{Vec{coord(i)}, key});
  }
}
 
 
template <class K, unsigned int dim>
  template <class Points>
void LexicographicPointSetBuilder<K,dim>::buildQuad (int a, int b, Points& points)
{
  static_assert(dim == 2);
  using Point = typename Points::value_type;
  using Vec = typename Point::Vector;
  using Key = Dune::LocalKey;
 
  auto coord_a = [a](int i) { return a == 0 ? K(0.5) : K(i)/a; };
  auto coord_b = [b](int i) { return b == 0 ? K(0.5) : K(i)/b; };
 
  for (int j = 0; j < b+1; ++j) {
    for (int i = 0; i < a+1; ++i) {
      auto key = ((i == 0 || i == a) && (j == 0 || j == b))
        ? Key((i > 0) + 2u * (j > 0), 2, 0)  // vertex
        : (i == 0 || i == a) ? Key(0u + (i > 0), 1, j-1) // left/right edge
        : (j == 0 || j == b) ? Key(2u + (j > 0), 1, i-1) // bottom/top edge
        : Key(0, 0, i-1);   // face
 
      points.push_back(Point{Vec{coord_a(i), coord_b(j)}, key});
    }
  }
}
 
 
template <class K, unsigned int dim>
  template <class Points>
void LexicographicPointSetBuilder<K,dim>::buildHex (int a, int b, int c, Points& points)
{
  static_assert(dim == 3);
  using Point = typename Points::value_type;
  using Vec = typename Point::Vector;
  using Key = Dune::LocalKey;
 
  auto coord_a = [a](int i) { return a == 0 ? K(0.5) : K(i)/a; };
  auto coord_b = [b](int i) { return b == 0 ? K(0.5) : K(i)/b; };
  auto coord_c = [c](int i) { return c == 0 ? K(0.5) : K(i)/c; };
 
  for (int k = 0; k < c+1; ++k) {
    for (int j = 0; j < b+1; ++j) {
      for (int i = 0; i < a+1; ++i) {
        auto key = ((i == 0 || i == a) && (j == 0 || j == b) && (k == 0 | k == c))
          ? Key((i > 0) + 2u * (j > 0) + 4u * (k > 0), 3, 0)  // vertex
          : (i == 0 || i == a) ? Key(0u + (i > 0), 2, (k-1)*(b-1) + (j-1)) // left/right face
          : (j == 0 || j == b) ? Key(2u + (j > 0), 2, (k-1)*(a-1) + (i-1)) // front/rear face
          : (k == 0 || k == c) ? Key(4u + (k > 0), 2, (j-1)*(a-1) + (i-1)) // bottom/top face
          : (i == 0 && j == 0) ? Key(0,1,k-1) // edge 0
          : (i == a && j == 0) ? Key(1,1,k-1) // edge 1
          : (i == 0 && j == b) ? Key(2,1,k-1)
          : (i == a && j == b) ? Key(3,1,k-1)
          : (i == 0 && k == 0) ? Key(4,1,j-1)
          : (i == a && k == 0) ? Key(5,1,j-1)
          : (j == 0 && k == 0) ? Key(6,1,i-1)
          : (j == b && k == 0) ? Key(7,1,i-1)
          : (i == 0 && k == c) ? Key(8,1,j-1)
          : (i == a && k == c) ? Key(9,1,j-1)
          : (j == 0 && k == c) ? Key(10,1,i-1)
          : (j == b && k == c) ? Key(11,1,i-1) // edge 11
          : Key(0, 0, (k-1)*(b-1)*(a-1) + (j-1)*(a-1) + (i-1));   // cell
 
        points.push_back(Point{Vec{coord_a(i), coord_b(j), coord_c(k)}, key});
      }
    }
  }
}
 
}} // end namespace Dune::Vtk::Impl