conformal_mesh.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/geometry/factory.hpp>
#include <kernel/geometry/mesh_permutation.hpp>
#include <kernel/geometry/intern/facet_neighbors.hpp>
#include <kernel/geometry/intern/standard_index_refiner.hpp>
#include <kernel/geometry/intern/standard_vertex_refiner.hpp>
#include <kernel/geometry/index_calculator.hpp>
#include <kernel/geometry/facet_flipper.hpp>
#include <kernel/util/math.hpp>
 
// includes, system
#include <memory>
 
namespace FEAT
{
  namespace Geometry
  {
    template<
      typename Shape_,
      int num_coords_ = Shape_::dimension,
      typename Coord_ = Real>
    class ConformalMesh
    {
      static_assert(num_coords_ >= Shape_::dimension, "invalid number of coordinates");
 
    public:
      typedef Shape_ ShapeType;
 
      typedef Coord_ CoordType;
 
      typedef VertexSet<num_coords_, Coord_> VertexSetType;
 
      typedef typename VertexSetType::VertexType VertexType;
 
      typedef IndexSetHolder<ShapeType> IndexSetHolderType;
 
      static constexpr int shape_dim = ShapeType::dimension;
      static constexpr int world_dim = VertexSetType::num_coords;
 
      static constexpr bool is_structured = false;
 
      typedef MeshPermutation<ShapeType> MeshPermutationType;
 
      template<
        int cell_dim_,
        int face_dim_>
      struct IndexSet
      {
        static_assert(cell_dim_ <= shape_dim, "invalid cell dimension");
        static_assert(face_dim_ < cell_dim_, "invalid face/cell dimension");
        static_assert(face_dim_ >= 0, "invalid face dimension");
 
        typedef FEAT::Geometry::IndexSet
        <
          Shape::FaceTraits
          <
            typename Shape::FaceTraits< ShapeType, cell_dim_> ::ShapeType,
            face_dim_
          >
          ::count
        > Type;
      }; // struct IndexSet<...>
 
      typedef typename IndexSet<shape_dim, shape_dim-1>::Type NeighborSetType;
 
    protected:
      Index _num_entities[shape_dim + 1];
 
      VertexSetType _vertex_set;
 
      IndexSetHolderType _index_set_holder;
 
      NeighborSetType _neighbors;
 
      MeshPermutationType _permutation;
 
    public:
      explicit ConformalMesh(const Index num_entities[]) :
        _vertex_set(num_entities[0]),
        _index_set_holder(num_entities),
        _neighbors(num_entities[shape_dim]),
        _permutation()
      {
        for(int i(0); i <= shape_dim; ++i)
        {
          //XASSERTM(num_entities[i] > 0, "Number of entities must not be zero!");
          _num_entities[i] = num_entities[i];
        }
      }
 
      explicit ConformalMesh(Factory<ConformalMesh>& factory) :
        _vertex_set(factory.get_num_entities(0)),
        _index_set_holder(Intern::NumEntitiesWrapper<shape_dim>(factory).num_entities),
        _neighbors(Intern::NumEntitiesWrapper<shape_dim>(factory).num_entities[shape_dim]),
        _permutation()
      {
        // Compute entity counts
        Intern::NumEntitiesWrapper<shape_dim>::apply(factory, _num_entities);
 
        // fill vertex set
        factory.fill_vertex_set(_vertex_set);
 
        // fill index sets
        factory.fill_index_sets(_index_set_holder);
 
        // Recompute entity counts. This is mainly for the MeshStreamerFactory, as complete information about the
        // number of edges/faces might not be known until after fill_index_sets is called
        Intern::NumEntitiesWrapper<shape_dim>::apply(factory, _num_entities);
 
        // Fill neighbor information. This needs facet at cell information, so it needs to be called after
        // fill_index_sets() etc.
        fill_neighbors();
      }
 
      ConformalMesh(ConformalMesh&& other) :
        _vertex_set(std::forward<VertexSetType>(other._vertex_set)),
        _index_set_holder(std::forward<IndexSetHolderType>(other._index_set_holder)),
        _neighbors(std::forward<NeighborSetType>(other._neighbors)),
        _permutation(std::forward<MeshPermutationType>(other._permutation))
      {
        for(int i(0); i <= shape_dim; ++i)
        {
          _num_entities[i] = other.get_num_entities(i);
        }
      }
 
      ConformalMesh& operator=(ConformalMesh&& other)
      {
        // avoid self-move
        if(this == &other)
          return *this;
 
        _vertex_set = std::forward<VertexSetType>(other._vertex_set);
        _index_set_holder = std::forward<IndexSetHolderType>(other._index_set_holder);
        _neighbors = std::forward<NeighborSetType>(other._neighbors);
        _permutation = std::forward<MeshPermutationType>(other._permutation);
 
        for(int i(0); i <= shape_dim; ++i)
        {
          _num_entities[i] = other.get_num_entities(i);
        }
 
        return *this;
      }
 
      ConformalMesh(const ConformalMesh&) = delete;
 
      ConformalMesh& operator=(const ConformalMesh&) = delete;
 
      virtual ~ConformalMesh()
      {
      }
 
      void clone(const ConformalMesh& other)
      {
        for(int i(0); i <= shape_dim; ++i)
          this->_num_entities[i] = other._num_entities[i];
        this->_vertex_set = other._vertex_set.clone();
        this->_index_set_holder.clone(other._index_set_holder);
        this->_neighbors = other._neighbors.clone();
        this->_permutation.clone(other._permutation);
      }
 
      ConformalMesh clone() const
      {
        ConformalMesh mesh(this->_num_entities);
        mesh._vertex_set = this->_vertex_set.clone();
        mesh._index_set_holder.clone(this->_index_set_holder);
        mesh._neighbors = this->_neighbors.clone();
        mesh._permutation.clone(this->_permutation);
        return mesh;
      }
 
      std::size_t bytes() const
      {
        return _vertex_set.bytes() + _index_set_holder.bytes() + _neighbors.bytes() + _permutation.bytes();
      }
 
      Index get_num_entities(int dim) const
      {
        XASSERT(dim >= 0);
        XASSERT(dim <= shape_dim);
        return _num_entities[dim];
      }
 
      Index get_num_vertices() const
      {
        return _num_entities[0];
      }
 
      Index get_num_elements() const
      {
        return _num_entities[shape_dim];
      }
 
      bool is_permuted() const
      {
        return !this->_permutation.empty();
      }
 
      const MeshPermutationType& get_mesh_permutation() const
      {
        return this->_permutation;
      }
 
      void create_permutation(PermutationStrategy strategy)
      {
        // make sure that we don't already have a permutation
        XASSERTM(this->_permutation.empty(), "mesh is already permuted!");
 
        // create the permutation
        this->_permutation.create(strategy, this->_index_set_holder, this->_vertex_set);
 
        // permute vertex set
        this->_vertex_set.permute(this->_permutation.get_perm(0));
 
        // permute index sets
        this->_index_set_holder.permute(this->_permutation.get_perms(), this->_permutation.get_inv_perms());
 
        // rebuild neighbor index set
        fill_neighbors();
      }
 
      void set_permutation(MeshPermutationType&& mesh_perm)
      {
        // make sure that we don't already have a permutation
        XASSERTM(this->_permutation.empty(), "mesh is already permuted!");
 
        // check the dimensions
        XASSERTM(mesh_perm.validate_sizes(this->_num_entities) == 0, "mesh permutation has invalid size!");
 
        // save the permutation
        this->_permutation = std::forward<MeshPermutationType>(mesh_perm);
 
        // permute vertex set
        this->_vertex_set.permute(this->_permutation.get_perm(0));
 
        // permute index sets
        this->_index_set_holder.permute(this->_permutation.get_perms(), this->_permutation.get_inv_perms());
 
        // rebuild neighbor index set
        fill_neighbors();
      }
 
      bool validate_element_coloring() const
      {
        // no coloring?
        const std::vector<Index>& coloring = this->get_mesh_permutation().get_element_coloring();
        if(coloring.empty())
          return true;
 
        // get vertices-at-element index set
        const auto& verts_at_elem = this->template get_index_set<shape_dim, 0>();
 
        // render transpose
        Adjacency::Graph elems_at_vert(Adjacency::RenderType::transpose, verts_at_elem);
 
        // loop over all color blocks
        for(std::size_t icol(0); icol+1u < coloring.size(); ++icol)
        {
          // get the bounds of our current color block
          const Index iel_beg = coloring[icol];
          const Index iel_end = coloring[icol+1u];
 
          // loop over all elements in the current color block
          for(Index iel(iel_beg); iel < iel_end; ++iel)
          {
            // loop over all vertices adjacent to this element
            for(int ivt(0); ivt < verts_at_elem.num_indices; ++ivt)
            {
              // loop over all elements adjacent to this vertex
              const Index ivtx = verts_at_elem(iel, ivt);
              for(auto it = elems_at_vert.image_begin(ivtx); it != elems_at_vert.image_end(ivtx); ++it)
              {
                // two adjacent element must not be in the same color block
                if((iel_beg <= *it) && (*it < iel_end) && (*it != iel))
                  return false; // invalid coloring
              }
            }
          }
        } // next color block
 
        // ok, coloring is valid
        return true;
      }
 
      bool validate_element_layering() const
      {
        // no layering?
        const std::vector<Index>& layering = this->get_mesh_permutation().get_element_layering();
        if(layering.empty())
          return true;
 
        // get vertices-at-element index set
        const auto& verts_at_elem = this->template get_index_set<shape_dim, 0>();
 
        // render transpose
        Adjacency::Graph elems_at_vert(Adjacency::RenderType::transpose, verts_at_elem);
 
        // loop over all layers
        for(std::size_t ilay(0); ilay+1u < layering.size(); ++ilay)
        {
          // get the bounds of our current layer
          const Index iel_beg = layering[ilay];
          const Index iel_end = layering[ilay+1u];
 
          // get the lower bound for valid neighbors of our current layer = beginning of previous layer
          const Index iel_lower = layering[Math::max(ilay, std::size_t(1)) - 1u];
 
          // get the upper bound for valid neighbors of our current layer = end of next layer
          const Index iel_upper = layering[Math::min(ilay+2u, layering.size()-1u)];
 
          // loop over all elements in the current layer
          for(Index iel(iel_beg); iel < iel_end; ++iel)
          {
            // loop over all vertices adjacent to this element
            for(int ivt(0); ivt < verts_at_elem.num_indices; ++ivt)
            {
              // loop over all elements adjacent to this vertex
              const Index ivtx = verts_at_elem(iel, ivt);
              for(auto it = elems_at_vert.image_begin(ivtx); it != elems_at_vert.image_end(ivtx); ++it)
              {
                // adjacent element outside of adjacent layers?
                if(!((iel_lower <= *it) && (*it < iel_upper)))
                  return false; // invalid layer
              }
            }
          }
        } // next color block
 
        // ok, layering is valid
        return true;
      }
 
      void fill_neighbors()
      {
        // Facet at cell index set
        auto& facet_idx = get_index_set<shape_dim, shape_dim -1>();
 
        XASSERTM(get_num_entities(shape_dim-1) == facet_idx.get_index_bound(), "mesh num_entities / index_set num_entities mismatch");
 
        if(_neighbors.get_num_indices() == Index(0))
          _neighbors = std::move(typename IndexSet<shape_dim, shape_dim-1>::Type(get_num_entities(shape_dim)));
 
        Intern::FacetNeighbors::compute(_neighbors, facet_idx);
      }
 
      typename IndexSet<shape_dim, shape_dim-1>::Type& get_neighbors()
      {
        return _neighbors;
      }
 
      const typename IndexSet<shape_dim, shape_dim-1>::Type& get_neighbors() const
      {
        return _neighbors;
      }
 
      VertexSetType& get_vertex_set()
      {
        return _vertex_set;
      }
 
      const VertexSetType& get_vertex_set() const
      {
        return _vertex_set;
      }
 
      template<
        int cell_dim_,
        int face_dim_>
      typename IndexSet<cell_dim_, face_dim_>::Type& get_index_set()
      {
        return _index_set_holder.template get_index_set_wrapper<cell_dim_>().template get_index_set<face_dim_>();
      }
 
      template<
        int cell_dim_,
        int face_dim_>
      const typename IndexSet<cell_dim_, face_dim_>::Type& get_index_set() const
      {
        return _index_set_holder.template get_index_set_wrapper<cell_dim_>().template get_index_set<face_dim_>();
      }
 
      IndexSetHolderType& get_index_set_holder()
      {
        return _index_set_holder;
      }
 
      const IndexSetHolderType& get_index_set_holder() const
      {
        return _index_set_holder;
      }
 
      IndexSetHolderType* get_topology()
      {
        return &_index_set_holder;
      }
 
      const IndexSetHolderType* get_topology() const
      {
        return &_index_set_holder;
      }
 
      void deduct_topology_from_top()
      {
        RedundantIndexSetBuilder<ShapeType>::compute(_index_set_holder);
        NumEntitiesExtractor<shape_dim>::set_num_entities(_index_set_holder, _num_entities);
        this->fill_neighbors();
        this->reorient_boundary_facets();
      }
 
      void reorient_boundary_facets()
      {
        Geometry::FacetFlipper<ShapeType>::reorient(this->_index_set_holder);
      }
 
      void transform(const VertexType& origin, const VertexType& angles, const VertexType& offset)
      {
        _vertex_set.transform(origin, angles, offset);
      }
 
      static String name()
      {
        return "ConformalMesh<"+Shape_::name()+","+stringify(num_coords_)+">";
      }
    }; // class ConformalMesh<...>
 
    /* ************************************************************************************************************* */
 
    template<
      typename Shape_,
      int num_coords_,
      typename CoordType_>
    class Factory< ConformalMesh<Shape_, num_coords_, CoordType_> >
    {
    public:
      typedef ConformalMesh<Shape_, num_coords_, CoordType_> MeshType;
 
      typedef typename MeshType::VertexSetType VertexSetType;
      typedef typename MeshType::IndexSetHolderType IndexSetHolderType;
 
    public:
      virtual ~Factory()
      {
      }
 
      virtual Index get_num_entities(int dim) = 0;
 
      virtual void fill_vertex_set(VertexSetType& vertex_set) = 0;
 
      virtual void fill_index_sets(IndexSetHolderType& index_set_holder) = 0;
 
      MeshType make()
      {
        return MeshType(*this);
      }
 
      std::unique_ptr<MeshType> make_unique()
      {
        return std::unique_ptr<MeshType>(new MeshType(*this));
      }
    }; // class Factory<ConformalMesh<...>>
 
    /* ************************************************************************************************************* */
 
    template<
      typename Shape_,
      int num_coords_,
      typename CoordType_>
    class StandardRefinery<ConformalMesh<Shape_, num_coords_, CoordType_> > :
      public Factory< ConformalMesh<Shape_, num_coords_, CoordType_> >
    {
    public:
      typedef ConformalMesh<Shape_, num_coords_, CoordType_> MeshType;
      typedef typename MeshType::ShapeType ShapeType;
      typedef typename MeshType::VertexSetType VertexSetType;
      typedef typename MeshType::IndexSetHolderType IndexSetHolderType;
 
      static constexpr int shape_dim = ShapeType::dimension;
 
    protected:
      const MeshType& _coarse_mesh;
      Index _num_entities_coarse[shape_dim + 1];
      Index _num_entities_fine[shape_dim + 1];
 
    public:
      explicit StandardRefinery(const MeshType& coarse_mesh) :
        _coarse_mesh(coarse_mesh)
      {
        // get number of entities in coarse mesh
        for(int i(0); i <= shape_dim; ++i)
        {
          _num_entities_fine[i] = coarse_mesh.get_num_entities(i);
          _num_entities_coarse[i] = coarse_mesh.get_num_entities(i);
        }
 
        // calculate number of entities in fine mesh
        Intern::EntityCountWrapper<Intern::StandardRefinementTraits, ShapeType>::query(_num_entities_fine);
      }
 
      virtual ~StandardRefinery()
      {
      }
 
      virtual Index get_num_entities(int dim) override
      {
        return _num_entities_fine[dim];
      }
 
      virtual void fill_vertex_set(VertexSetType& vertex_set) override
      {
        // refine vertices
        Intern::StandardVertexRefineWrapper<ShapeType, VertexSetType>
          ::refine(vertex_set, _coarse_mesh.get_vertex_set(), _coarse_mesh.get_index_set_holder());
      }
 
      virtual void fill_index_sets(IndexSetHolderType& index_set_holder) override
      {
        // refine indices
        Intern::IndexRefineWrapper<ShapeType>
          ::refine(index_set_holder, _num_entities_coarse, _coarse_mesh.get_index_set_holder());
      }
    }; // class StandardRefinery<ConformalMesh<...>,Nil>
 
#ifdef FEAT_EICKT
    extern template class ConformalMesh<Shape::Simplex<2>, 2, Real>;
    extern template class ConformalMesh<Shape::Simplex<3>, 3, Real>;
    extern template class ConformalMesh<Shape::Hypercube<2>, 2, Real>;
    extern template class ConformalMesh<Shape::Hypercube<3>, 3, Real>;
 
    extern template class StandardRefinery<ConformalMesh<Shape::Simplex<2>, 2, Real>>;
    extern template class StandardRefinery<ConformalMesh<Shape::Simplex<3>, 3, Real>>;
    extern template class StandardRefinery<ConformalMesh<Shape::Hypercube<2>, 2, Real>>;
    extern template class StandardRefinery<ConformalMesh<Shape::Hypercube<3>, 3, Real>>;
#endif // FEAT_EICKT
  } // namespace Geometry
} // namespace FEAT