shape.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
// includes, FEAT
#include <kernel/base_header.hpp>
#include <kernel/util/meta_math.hpp>
#include <kernel/util/string.hpp>
 
namespace FEAT
{
  namespace Shape
  {
    struct Vertex
    {
      static constexpr int dimension = 0;
 
      static String name()
      {
        return "Vertex";
      }
    }; // struct Vertex
 
    template<int dimension_>
    struct Simplex
    {
      static_assert(dimension_ > 0, "parameter dimension_ must be greater than 0");
 
      static constexpr int dimension = dimension_;
 
      static String name()
      {
        return "Simplex<" + stringify(dimension_) + ">";
      }
    }; // class Simplex
 
    template<int dimension_>
    struct Hypercube
    {
      static_assert(dimension_ > 0, "parameter dimension_ must be greater than 0");
 
      static constexpr int dimension = dimension_;
 
      static String name()
      {
        return "Hypercube<" + stringify(dimension_) + ">";
      }
    }; // struct Hypercube
 
    typedef Simplex<2> Triangle;
    typedef Simplex<3> Tetrahedron;
 
    typedef Hypercube<2> Quadrilateral;
    typedef Hypercube<3> Hexahedron;
 
    template<
      typename Shape_,
      int face_dim_>
#ifndef DOXYGEN
    struct FaceTraits;
#else
    struct FaceTraits
    {
      typedef ... ShapeType;
 
      static constexpr int count = ...
    };
#endif // DOXYGEN
 
 
    template<>
    struct FaceTraits<Vertex, 0>
    {
      typedef Vertex ShapeType;
 
      static constexpr int count = 1;
    }; // struct FaceTraits<Vertex,0>
 
    template<
      int cell_dim_,
      int face_dim_>
    struct FaceTraits<Simplex<cell_dim_>, face_dim_>
    {
      static_assert(face_dim_ > 0, "parameter face_dim_ must be greater than 0");
      static_assert(cell_dim_ >= face_dim_, "face_dim_ must not be greater than cell_dim_");
 
      typedef Simplex<face_dim_> ShapeType;
 
      static constexpr int count = MetaMath::Binomial<cell_dim_ + 1, face_dim_ + 1>::value;
    }; // struct FaceTraits<Simplex<...>, ...>
 
    template<int cell_dim_>
    struct FaceTraits<Simplex<cell_dim_>, 0>
    {
      static_assert(cell_dim_ > 0, "parameter cell_dim_ must be greater than 0");
 
      typedef Vertex ShapeType;
 
      static constexpr int count = cell_dim_ + 1;
    }; // struct FaceTraits<Simplex<...>, 0>
 
    template<
      int cell_dim_,
      int face_dim_>
    struct FaceTraits<Hypercube<cell_dim_>, face_dim_>
    {
      static_assert(face_dim_ > 0, "parameter face_dim_ must be greater than 0");
      static_assert(cell_dim_ >= face_dim_, "face_dim_ must not be greater than cell_dim_");
 
      typedef Hypercube<face_dim_> ShapeType;
 
      static constexpr int count = (1 << (cell_dim_ - face_dim_)) * MetaMath::Binomial<cell_dim_, face_dim_>::value;
    }; // struct FaceTraits<Hypercube<...>, ...>
 
    template<int cell_dim_>
    struct FaceTraits<Hypercube<cell_dim_>, 0>
    {
      static_assert(cell_dim_ > 0, "parameter cell_dim_ must be greater than 0");
 
      typedef Vertex ShapeType;
 
      static constexpr int count = (1 << cell_dim_);
    }; // struct FaceTraits<HyperCube<...>, 0>
 
    template<typename Shape_>
#ifndef DOXYGEN
    struct ReferenceCell;
#else
    struct ReferenceCell
    {
      template<typename T_>
      static T_ vertex(int vertex_idx, int coord_idx);
 
      template<typename T_>
      static T_ centre(int coord_idx);
 
      template<typename T_>
      static T_ volume();
 
      static int facet_orientation(int facet_index);
    };
#endif // DOXYGEN
 
 
    template<>
    struct ReferenceCell<Shape::Vertex>
    {
      template<typename T_>
      static T_ vertex(int, int)
      {
        return T_(0);
      }
 
      template<typename T_>
      static T_ centre(int)
      {
        return T_(0);
      }
 
      template<typename T_>
      static T_ volume()
      {
        return T_(0);
      }
 
      // orientation is not implemented due to being nonsensical
 
    };
 
    template<int dim_>
    struct ReferenceCell< Shape::Simplex<dim_> >
    {
      template<typename T_>
      static T_ vertex(int vertex_idx, int coord_idx)
      {
        return (coord_idx + 1) == vertex_idx ? T_(1) : T_(0);
      }
 
      template<typename T_>
      static T_ centre(int)
      {
        return T_(1) / T_(dim_+1);
      }
 
      template<typename T_>
      static T_ volume()
      {
        return T_(1) / T_(MetaMath::Factorial<dim_>::value);
      }
 
      static int facet_orientation(int facet_index)
      {
        ASSERTM((facet_index >= 0), "facet index out of range!");
        ASSERTM((facet_index < FaceTraits<Simplex<dim_>, dim_-1>::count), "facet index out of range!");
        return 1 - (dim_ == 2 ? 0 : dim_ == 3 ? (facet_index & 1) << 1 : ((facet_index & 1) ^ 1) << 1);
      }
    };
 
    template<int dim_>
    struct ReferenceCell< Shape::Hypercube<dim_> >
    {
      template<typename T_>
      static T_ vertex(int vertex_idx, int coord_idx)
      {
        return T_((((vertex_idx >> coord_idx) & 1) << 1) - 1);
      }
 
      template<typename T_>
      static T_ centre(int)
      {
        return T_(0);
      }
 
      template<typename T_>
      static T_ volume()
      {
        return T_(1 << dim_); // = 2^dim
      }
 
      static int facet_orientation(int facet_index)
      {
        ASSERTM((facet_index >= 0), "facet index out of range!");
        ASSERTM((facet_index < FaceTraits<Hypercube<dim_>,dim_-1>::count), "facet index out of range!");
 
        return 1 - (( ((facet_index >> 1) ^ facet_index ^ dim_) & 1 ) << 1);
      }
    };
 
  } // namespace Shape
} // namespace FEAT