Generalize vertex/line/quad_index

[peterrum]

Rely in internal::TriaAccessorImplementation::Implementation::quad_index(), internal::TriaAccessorImplementation::Implementation::line_index(), and internal::TriaAccessorImplementation::Implementation::vertex_index() completely on GeometryInfo. Still to do (in follow-up PRs):

• make GeometryInfo dynamic
• merge orientation/flip/rotation to a single number

internal::TriaAccessorImplementation::Implementation::quad_index():

      template <int dim, int spacedim>
      inline static unsigned int
      quad_index(const TriaAccessor<3, dim, spacedim> &accessor,
                 const unsigned int                    i)
      {
        return accessor.tria->levels[accessor.present_level]
          ->cells.cells[accessor.present_index]
          .face(i);
      }

internal::TriaAccessorImplementation::Implementation::line_index():

      template <int dim, int spacedim>
      inline static unsigned int
      line_index(const TriaAccessor<2, dim, spacedim> &accessor,
                 const unsigned int                    i)
      {
        return accessor.objects().cells[accessor.present_index].face(i);
      }

      template <int dim, int spacedim>
      inline static unsigned int
      line_index(const TriaAccessor<3, dim, spacedim> &accessor,
                 const unsigned int                    i)
      {
        const auto pair = GeometryInfo<3>::standard_line_to_quad_line_index(i);
        const auto quad_index = pair[0];
        const auto line_index = GeometryInfo<3>::standard_to_real_face_line(
          pair[1],
          accessor.face_orientation(quad_index),
          accessor.face_flip(quad_index),
          accessor.face_rotation(quad_index));

        return accessor.quad(quad_index)->line_index(line_index);
      }

internal::TriaAccessorImplementation::Implementation::vertex_index():

      template <int dim, int spacedim>
      inline static unsigned int
      vertex_index(const TriaAccessor<1, dim, spacedim> &accessor,
                   const unsigned int                    corner)
      {
        return accessor.objects().cells[accessor.present_index].face(corner);
      }

      template <int dim, int spacedim>
      inline static unsigned int
      vertex_index(const TriaAccessor<2, dim, spacedim> &accessor,
                   const unsigned int                    corner)
      {
        const auto pair =
          GeometryInfo<2>::standard_corner_to_line_vertex_index(corner);
        const auto line_index   = pair[0];
        const auto vertex_index = GeometryInfo<2>::standard_to_real_line_vertex(
          pair[1], accessor.line_orientation(line_index));

        return accessor.line(line_index)->vertex_index(vertex_index);
      }

      template <int dim, int spacedim>
      inline static unsigned int
      vertex_index(const TriaAccessor<3, dim, spacedim> &accessor,
                   const unsigned int                    corner)
      {
        const auto pair =
          GeometryInfo<3>::standard_corner_to_quad_vertex_index(corner);
        const auto face_index   = pair[0];
        const auto vertex_index = GeometryInfo<3>::standard_to_real_face_vertex(
          pair[1],
          accessor.face_orientation(face_index),
          accessor.face_flip(face_index),
          accessor.face_rotation(face_index));

        return accessor.quad(face_index)->vertex_index(vertex_index);
      }

[peterrum]

Extension to TET/TRI:

• recursive entity definition

std::pair<unsinged int, unsigned int>
GeometryInfo<3>::standard_line_to_quad_line_index(line)
{
  static const std::pair<unsigned int, unsigned int> table[6] =
    {{0, 2}, {0, 1}, {0, 0}, {1, 2}, {1, 1}, {2, 1}};

  return table[line];
}

std::pair<unsinged int, unsigned int>
GeometryInfo<2>::standard_corner_to_line_vertex_index(vertex)
{
  static const std::pair<unsigned int, unsigned int> table[3] =
    {{0, 2}, {0, 1}, {1, 1}};

  return table[vertex];
}

std::pair<unsinged int, unsigned int>
GeometryInfo<3>::standard_corner_to_quad_vertex_index(vertex)
{
  static const std::pair<unsigned int, unsigned int> table[4] =
    {{0, 0}, {0, 2}, {0, 1}, {1, 2}};

  return table[vertex];
}

• orientation

  unsigned int
  GeometryInfo<3>::standard_to_real_face_line(line, orientation)
  {
  static const unsigned int table[6][3] = 
  {{0, 1, 2}, {2, 1, 0}, {0, 2, 1}, {1, 2, 0}, {2, 0, 1}, {1, 0, 2}};
  
  return table[orientation][line];
  }

  unsigned int
  GeometryInfo<2>::standard_to_real_line_vertex(vertex, orientation)
  {
  static const unsigned int table[2][2] = {{1, 0}, {0, 1}};
  
  return table[orientation][vertex];
  }

unsigned int
GeometryInfo<3>::standard_corner_to_quad_vertex_index(vertex, orientation)
{
  static const unsigned int table[6][3] = 
    {{0, 1, 2}, {0, 2, 1}, {1, 0, 2}, {1, 2, 0}, {2, 0, 1}, {2, 1, 0}};

  return table[orientation][vertex];
}

Question: are all 6 orientations needed and/or possible?

[peterrum]

Integration of different entity types:

      template <int dim, int spacedim>
      inline static unsigned int
      vertex_index(const TriaAccessor<3, dim, spacedim> &accessor,
                   const unsigned int                    corner)
      {
        if(accessor->use_entity_table == false) // default case
          {
            const auto pair =
              GeometryInfo<3>::standard_corner_to_quad_vertex_index(corner);
            const auto face_index   = pair[0];
            const auto vertex_index = GeometryInfo<3>::standard_to_real_face_vertex(
              pair[1], accessor.orientation_raw(face_index));

            return accessor.quad(face_index)->vertex_index(vertex_index);
          }
        else
          {
            const auto & geometry_info = 
              entities[entity_type[accessor->index]]; // TODO

            const auto pair =
              geometry_info.standard_corner_to_quad_vertex_index(corner);
            const auto face_index   = pair[0];
            const auto vertex_index = geometry_info.standard_to_real_face_vertex(
              pair[1], accessor.orientation_raw(face_index));

            return accessor.quad(face_index)->vertex_index(vertex_index);
          }
      }

or alternatively:

      template <int dim, int spacedim, typename T>
      inline static unsigned int
      vertex_index(const TriaAccessor<3, dim, spacedim> &accessor,
                   const unsigned int                    corner,
                   const T                              &geometry_info)
      {
        const auto pair =
          geometry_info.standard_corner_to_quad_vertex_index(corner);
        const auto face_index   = pair[0];
        const auto vertex_index = geometry_info.standard_to_real_face_vertex(
          pair[1], accessor.orientation_raw(face_index));

        return accessor.quad(face_index)->vertex_index(vertex_index);
      } 

      template <int dim, int spacedim>
      inline static unsigned int
      vertex_index(const TriaAccessor<3, dim, spacedim> &accessor,
                   const unsigned int                    corner)
      {
        if(accessor->use_entity_table == false) // default case
          return vertex_index(accessor, corner, GeometryInfo<3>);
        else
          return vertex_index(accessor, corner, 
                              entities[entity_type[accessor->index]] /*TODO*/);
      }

Question 1: do we want to support different entity types?

Question 2: do we want a dynamic GeometryInfo?

Question 3: or should we use a jump table: hex vs. tet?