mesh_part.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/conformal_mesh.hpp>
#include <kernel/geometry/structured_mesh.hpp>
#include <kernel/geometry/index_calculator.hpp>
#include <kernel/geometry/attribute_set.hpp>
#include <kernel/geometry/intern/simple_target_refiner.hpp>
#include <kernel/geometry/intern/standard_target_refiner.hpp>
#include <kernel/geometry/intern/structured_target_refiner.hpp>
#include <kernel/geometry/intern/standard_attrib_refiner.hpp>
#include <kernel/geometry/intern/standard_index_refiner.hpp>
#include <kernel/geometry/intern/standard_vertex_refiner.hpp>
#include <kernel/geometry/intern/index_set_filler.hpp>
#include <kernel/geometry/intern/target_set_computer.hpp>
 
// includes, system
#include <map>
#include <array>
#include <memory>
 
namespace FEAT
{
  namespace Geometry
  {
    template<typename MeshType_>
    class MeshPart
    {
    public:
      typedef MeshType_ MeshType;
      typedef typename MeshType::ShapeType ShapeType;
      typedef typename MeshType::IndexSetHolderType ParentIndexSetHolderType;
      typedef IndexSetHolder<ShapeType> IndexSetHolderType;
      typedef TargetSetHolder<ShapeType> TargetSetHolderType;
      static constexpr int shape_dim = ShapeType::dimension;
 
      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 MeshType::VertexSetType::CoordType AttributeDataType;
      typedef AttributeSet<AttributeDataType> AttributeSetType;
 
      typedef std::map<String, std::unique_ptr<AttributeSetType>> AttributeSetContainer;
      typedef typename AttributeSetContainer::iterator AttributeSetIterator;
      typedef typename AttributeSetContainer::const_iterator AttributeSetConstIterator;
      typedef typename AttributeSetContainer::reverse_iterator AttributeSetReverseIterator;
 
      template<int cell_dim_>
      struct TargetSet
      {
        typedef FEAT::Geometry::TargetSet Type;
      }; // struct TargetSet<...>
 
    protected:
      std::array<Index, shape_dim + 1> _num_entities;
      std::unique_ptr<IndexSetHolderType> _index_set_holder;
      TargetSetHolderType _target_set_holder;
      AttributeSetContainer _mesh_attributes;
 
    public:
      explicit MeshPart(const Index num_entities[], bool create_topology = false) :
        _index_set_holder(),
        _target_set_holder(num_entities),
        _mesh_attributes()
      {
        for(std::size_t i(0); i <= shape_dim; ++i)
          _num_entities[i] = num_entities[i];
 
        if(create_topology)
          _index_set_holder.reset(new IndexSetHolderType(num_entities));
      }
 
      explicit MeshPart(Factory<MeshPart>& factory) :
        _index_set_holder(),
        _target_set_holder(Intern::NumEntitiesWrapper<shape_dim>(factory).num_entities),
        _mesh_attributes()
      {
        // Compute entity counts
        Intern::NumEntitiesWrapper<shape_dim>::apply(factory, _num_entities);
 
        // Fill target sets
        factory.fill_target_sets(_target_set_holder);
 
        // Fill index sets
        factory.fill_index_sets(_index_set_holder);
 
        // Fill attribute sets
        factory.fill_attribute_sets(_mesh_attributes);
      }
 
      MeshPart(MeshPart&& other) = default;
 
      MeshPart& operator=(MeshPart&& other) = default;
 
      MeshPart(const MeshPart& other) = delete;
 
      MeshPart& operator=(const MeshPart&) = delete;
 
      virtual ~MeshPart()
      {
        // unique_ptr takes care of everything
      }
 
      void clone(const MeshPart& other)
      {
        this->_num_entities = other._num_entities;
 
        if(other._index_set_holder)
        {
          if(this->_index_set_holder)
            this->_index_set_holder->clone(*other._index_set_holder);
          else
            this->_index_set_holder.reset(new IndexSetHolderType(other._index_set_holder->clone()));
        }
        else if(this->_index_set_holder)
        {
          this->_index_set_holder.reset();
        }
 
        this->_target_set_holder.clone(other._target_set_holder);
 
        this->_mesh_attributes.clear();
 
        for(auto it = other._mesh_attributes.begin(); it != other._mesh_attributes.end(); ++it)
          this->add_attribute(std::unique_ptr<AttributeSetType>(new AttributeSetType(it->second->clone())), it->first);
      }
 
      MeshPart clone() const
      {
        MeshPart mp(this->_num_entities.data(), this->_index_set_holder.get() != nullptr);
        if(this->_index_set_holder)
          mp._index_set_holder->clone(*this->_index_set_holder);
        mp._target_set_holder.clone(this->_target_set_holder);
        for(auto it = this->_mesh_attributes.begin(); it != this->_mesh_attributes.end(); ++it)
          mp.add_attribute(std::unique_ptr<AttributeSetType>(new AttributeSetType(it->second->clone())), it->first);
        return mp;
      }
 
      std::size_t bytes() const
      {
        std::size_t my_bytes(0);
 
        for(const auto& it : _mesh_attributes)
        {
          if(it.second != nullptr)
            my_bytes += it.second->bytes();
        }
 
        my_bytes += _target_set_holder.bytes();
        my_bytes += (_index_set_holder ? _index_set_holder->bytes() : std::size_t(0));
 
        return my_bytes;
      }
 
      bool has_topology() const
      {
        return _index_set_holder.operator bool();
      }
 
      Index get_num_entities(int dim) const
      {
        XASSERT(dim >= 0);
        XASSERT(dim <= shape_dim);
        return _num_entities[std::size_t(dim)];
      }
 
      AttributeSetType* find_attribute(const String& identifier)
      {
        AttributeSetIterator it(_mesh_attributes.find(identifier));
        return (it != _mesh_attributes.end()) ? (*it).second.get() : nullptr;
      }
 
      const AttributeSetType* find_attribute(const String& identifier) const
      {
        AttributeSetConstIterator it(_mesh_attributes.find(identifier));
        return (it != _mesh_attributes.end()) ? (*it).second.get(): nullptr;
      }
 
      const AttributeSetContainer& get_mesh_attributes() const
      {
        return _mesh_attributes;
      }
 
      int get_num_attributes() const
      {
        return int(_mesh_attributes.size());
      }
 
      virtual bool add_attribute(std::unique_ptr<AttributeSetType> attribute, const String& identifier)
      {
        XASSERTM(bool(attribute), "cannot add empty attribute set");
        // note: unique_ptr needs to be moved here, otherwise emplace tries to copy instead
        return _mesh_attributes.emplace(identifier, std::move(attribute)).second;
      }
 
      template<int cell_dim_, int face_dim_>
      typename IndexSet<cell_dim_, face_dim_>::Type& get_index_set()
      {
        XASSERTM(has_topology(), "Requested index_set of MeshPart without topology!");
        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
      {
        XASSERTM(has_topology(), "Requested index_set of MeshPart without topology!");
        return _index_set_holder->template get_index_set_wrapper<cell_dim_>().template get_index_set<face_dim_>();
      }
 
 
      IndexSetHolderType* get_topology()
      {
        return _index_set_holder.get();
      }
 
      const IndexSetHolderType* get_topology() const
      {
        return _index_set_holder.get();
      }
 
      template<int cell_dim_>
      typename TargetSet<cell_dim_>::Type& get_target_set()
      {
        static_assert(cell_dim_ >= 0, "invalid cell dimension");
        static_assert(cell_dim_ <= shape_dim, "invalid cell dimension");
        return _target_set_holder.template get_target_set<cell_dim_>();
      }
 
      template<int cell_dim_>
      const typename TargetSet<cell_dim_>::Type& get_target_set() const
      {
        static_assert(cell_dim_ >= 0, "invalid cell dimension");
        static_assert(cell_dim_ <= shape_dim, "invalid cell dimension");
        return _target_set_holder.template get_target_set<cell_dim_>();
      }
 
      TargetSetHolderType& get_target_set_holder()
      {
        return _target_set_holder;
      }
 
      const TargetSetHolderType& get_target_set_holder() const
      {
        return _target_set_holder;
      }
 
      static String name()
      {
        return "MeshPart<...>";
      }
 
      void permute(const MeshPermutation<ShapeType>& mesh_perm)
      {
        this->_target_set_holder.permute_map(mesh_perm.get_inv_perms());
      }
 
      template<int end_dim_, int current_dim_ = ShapeType::dimension>
      void deduct_target_sets_from_bottom(const ParentIndexSetHolderType& parent_ish)
      {
        Intern::template TargetSetComputer<end_dim_, current_dim_>::bottom_to_top(_target_set_holder, parent_ish);
 
        // Update num_entities information from the possibly modified _target_set_holder
        for(int i(end_dim_); i <= current_dim_; ++i)
          _num_entities[std::size_t(i)] = _target_set_holder.get_num_entities(i);
      }
 
      template<int end_dim_, int current_dim_ = 0>
      void deduct_target_sets_from_top(const ParentIndexSetHolderType& parent_ish)
      {
        Intern::template TargetSetComputer<end_dim_, current_dim_>::top_to_bottom(_target_set_holder, parent_ish);
 
        // Update num_entities information from the possibly modified _target_set_holder
        for(int i(current_dim_); i <= end_dim_; ++i)
          _num_entities[std::size_t(i)] = _target_set_holder.get_num_entities(i);
      }
 
      void deduct_topology(const ParentIndexSetHolderType& parent_ish)
      {
        if(_index_set_holder == nullptr)
          _index_set_holder.reset(new IndexSetHolderType(_num_entities));
 
        Intern::IndexSetFiller<ShapeType::dimension>::fill_ish(
          *_index_set_holder, _target_set_holder, parent_ish);
 
        // build redundant index sets
        RedundantIndexSetBuilder<ShapeType>::compute(*_index_set_holder);
      }
    }; // class MeshPart<ConformalMesh<Shape_>>
 
    template<typename MeshType_>
    class Factory< MeshPart<MeshType_> >
    {
    public:
      typedef typename MeshType_::ShapeType ShapeType;
      typedef MeshPart<MeshType_> MeshPartType;
      typedef typename MeshType_::VertexSetType::CoordType AttributeDataType;
      typedef typename MeshPartType::AttributeSetContainer AttributeSetContainer;
      typedef typename MeshPartType::IndexSetHolderType IndexSetHolderType;
      typedef typename MeshPartType::TargetSetHolderType TargetSetHolderType;
 
    public:
      virtual ~Factory()
      {
      }
 
      virtual Index get_num_entities(int dim) = 0;
 
      virtual void fill_attribute_sets(AttributeSetContainer& attribute_set_holder) = 0;
 
      virtual void fill_index_sets(std::unique_ptr<IndexSetHolderType>& index_set_holder) = 0;
 
      virtual void fill_target_sets(TargetSetHolderType& target_set_holder) = 0;
 
      MeshPartType make()
      {
        return MeshPartType(*this);
      }
 
      std::unique_ptr<MeshPartType> make_unique()
      {
        return std::unique_ptr<MeshPartType>(new MeshPartType(*this));
      }
    }; // class Factory<MeshPart<...>>
 
    namespace Intern
    {
      template<typename ParentMesh_>
      struct TargetSetRefineParentWrapper;
    }
 
    template<typename ParentMesh_>
    class StandardRefinery< MeshPart<ParentMesh_> > :
      public Factory< MeshPart<ParentMesh_> >
    {
    public:
      typedef MeshPart<ParentMesh_> MeshType;
      typedef typename MeshType::ShapeType ShapeType;
      typedef typename MeshType::AttributeSetType AttributeType;
      typedef typename MeshType::AttributeSetContainer AttributeSetContainer;
      typedef typename MeshType::IndexSetHolderType IndexSetHolderType;
      typedef typename MeshType::TargetSetHolderType TargetSetHolderType;
 
      static constexpr int shape_dim = ShapeType::dimension;
 
    protected:
      const MeshType& _coarse_meshpart;
      const ParentMesh_* _parent_mesh;
      const MeshType* _parent_meshpart;
      Index _num_entities_coarse[shape_dim + 1];
      Index _num_entities_fine[shape_dim + 1];
 
    public:
      explicit StandardRefinery(const MeshType& coarse_meshpart, const ParentMesh_& parent_mesh) :
        _coarse_meshpart(coarse_meshpart),
        _parent_mesh(&parent_mesh),
        _parent_meshpart(nullptr)
      {
        // get number of entities in coarse mesh
        for(int i(0); i <= shape_dim; ++i)
        {
          _num_entities_fine[i] = _num_entities_coarse[i] = coarse_meshpart.get_num_entities(i);
        }
 
        // calculate number of entities in fine mesh
        Intern::EntityCountWrapper<Intern::StandardRefinementTraits, ShapeType>::query(_num_entities_fine);
      }
 
      explicit StandardRefinery(const MeshType& coarse_meshpart, const MeshPart<ParentMesh_>& parent_meshpart) :
        _coarse_meshpart(coarse_meshpart),
        _parent_mesh(nullptr),
        _parent_meshpart(&parent_meshpart)
      {
        // get number of entities in coarse mesh
        for(int i(0); i <= shape_dim; ++i)
        {
          _num_entities_fine[i] = _num_entities_coarse[i] = coarse_meshpart.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_attribute_sets(AttributeSetContainer& attribute_container) override
      {
        // Attributes of shape dimension 0 only make sense if we have a mesh topology
        if(!_coarse_meshpart.has_topology())
          return;
 
        // Iterate over the attributes in the coarse mesh
        typename MeshType::AttributeSetConstIterator it(_coarse_meshpart.get_mesh_attributes().begin());
        typename MeshType::AttributeSetConstIterator jt(_coarse_meshpart.get_mesh_attributes().end());
 
        for(; it != jt; ++it)
        {
          // Create a new empty attribute of the desired size
          std::unique_ptr<AttributeType> refined_attribute(new AttributeType(
            get_num_entities(0), it->second->get_dimension()));
 
          // Refine the attribute in the coarse mesh and write the result to the new attribute
          Intern::StandardAttribRefineWrapper<ShapeType, AttributeType>
            ::refine(*refined_attribute, *(it->second), *_coarse_meshpart.get_topology());
 
          // Add the attribute to the corresponding set
          XASSERTM(attribute_container.emplace(it->first, std::move(refined_attribute)).second, "Error refining attribute " + it->first);
        }
      }
 
      virtual void fill_index_sets(std::unique_ptr<IndexSetHolderType>& index_set_holder) override
      {
        XASSERT(index_set_holder.get() == nullptr);
 
        // Only create the topology for the refined mesh if the coarse mesh has a topology
        if(!_coarse_meshpart.has_topology())
          return;
 
        index_set_holder.reset(new IndexSetHolderType(_num_entities_fine));
 
        // refine indices
        Intern::IndexRefineWrapper<ShapeType>::refine(
          *index_set_holder, _num_entities_coarse, *_coarse_meshpart.get_topology());
      }
 
      virtual void fill_target_sets(TargetSetHolderType& target_set_holder) override
      {
        // We have to use a helper class here, because the target set refinement
        // depends heavily on the underlying mesh type. For this, we first need
        // to check whether our parent is the root mesh or another mesh-part.
        if(_parent_mesh != nullptr)
        {
          // the parent is the root mesh
          Intern::TargetSetRefineParentWrapper<ParentMesh_>::fill_target_sets(target_set_holder,
            _coarse_meshpart.get_target_set_holder(), _coarse_meshpart.get_topology(), *_parent_mesh);
        }
        else if(_parent_meshpart != nullptr)
        {
          // the parent is a mesh-part
          Intern::TargetSetRefineParentWrapper<ParentMesh_>::fill_target_sets(target_set_holder,
            _coarse_meshpart.get_target_set_holder(), _coarse_meshpart.get_topology(), *_parent_meshpart);
        }
        else
          XABORTM("no parent present");
      }
 
    }; // class StandardRefinery<MeshPart<...>,...>
 
    namespace Intern
    {
      template<typename Shape_, int num_coords_, typename Coord_>
      struct TargetSetRefineParentWrapper<ConformalMesh<Shape_, num_coords_, Coord_>>
      {
        static constexpr int shape_dim = Shape_::dimension;
 
        // Note: This functions works for both mesh parents and mesh-part parents.
        template<typename TargetSetHolderType_, typename IndexSetHolder_, typename Parent_>
        static void fill_target_sets(
          TargetSetHolderType_& target_set_holder,
          const TargetSetHolderType_& coarse_target_set_holder,
          const IndexSetHolder_* coarse_ish,
          const Parent_& parent)
        {
          // get number of parent entities
          Index num_entities_parent[shape_dim+1];
          for(int i(0); i <= shape_dim; ++i)
            num_entities_parent[i] = parent.get_num_entities(i);
 
          if(coarse_ish != nullptr)
          {
            // The coarse mesh-part has a topology, so refine the target set to match.
            // The parent must have a topology, too, which is non-trivial if the
            // parent is a mesh-part itself.
            const auto* parent_topo = parent.get_topology();
            XASSERTM(parent_topo != nullptr, "mesh-part has topology, but parent doesn't");
            Intern::TargetRefineWrapper<Shape_>::refine(target_set_holder,
              num_entities_parent, coarse_target_set_holder, *coarse_ish, *parent_topo);
          }
          else
          {
            // The coarse mesh-part does not have a topology, stick to simple refinement
            Intern::SimpleTargetRefineWrapper<Shape_>::refine(
              target_set_holder, num_entities_parent, coarse_target_set_holder);
          }
        }
      }; // // struct TargetSetRefineParentWrapper<ConformalMesh<...>>
 
      template<int shape_dim_, int num_coords_, typename Coord_>
      struct TargetSetRefineParentWrapper<StructuredMesh<shape_dim_, num_coords_, Coord_>>
      {
        typedef Shape::Hypercube<shape_dim_> ShapeType;
 
        template<typename TargetSetHolderType_, typename IndexSetHolder_>
        static void fill_target_sets(
          TargetSetHolderType_& target_set_holder,
          const TargetSetHolderType_& coarse_target_set_holder,
          const IndexSetHolder_* coarse_ish,
          const StructuredMesh<shape_dim_, num_coords_, Coord_>& parent)
        {
          // get number of coarse/fine parent slices
          Index num_slices_c[shape_dim_], num_slices_f[shape_dim_];
          for(int i(0); i < shape_dim_; ++i)
          {
            num_slices_c[i] = parent.get_num_slices(i);
            num_slices_f[i] = Index(2) * num_slices_c[i];
          }
 
          if(coarse_ish != nullptr)
          {
            XABORTM("TargetSetRefineParentWrapper not implemented");
            //Intern::TargetRefineWrapper<ShapeType>::refine(
              //target_set_holder, num_entities_parent, coarse_target_set_holder, *coarse_ish, parent_ish);
          }
          else
          {
            Intern::StructuredTargetRefineWrapper<shape_dim_>::refine_simple(
              target_set_holder, coarse_target_set_holder, num_slices_c, num_slices_f);
          }
        }
 
        // Specialization for MeshPart<StructuredMesh<...>>
        // Note: The topology of a MeshPart is always "conforming" (ie "unstructured"),
        // therefore the implementation is identical to the class template specialization
        // for ConformalMesh above.
        template<typename TargetSetHolderType_, typename IndexSetHolder_>
        static void fill_target_sets(
          TargetSetHolderType_& target_set_holder,
          const TargetSetHolderType_& coarse_target_set_holder,
          const IndexSetHolder_* coarse_ish,
          const MeshPart<StructuredMesh<shape_dim_, num_coords_, Coord_>>& parent)
        {
          // get number of parent entities
          Index num_entities_parent[shape_dim_+1];
          for(int i(0); i <= shape_dim_; ++i)
            num_entities_parent[i] = parent.get_num_entities(i);
 
          if(coarse_ish != nullptr)
          {
            // The coarse mesh-part has a topology, so refine the target set to match.
            // The parent must have a topology, too, which is non-trivial if the
            // parent is a mesh-part itself.
            const auto* parent_topo = parent.get_topology();
            XASSERTM(parent_topo != nullptr, "mesh-part has topology, but parent doesn't");
            Intern::TargetRefineWrapper<ShapeType>::refine(target_set_holder,
              num_entities_parent, coarse_target_set_holder, *coarse_ish, *parent_topo);
          }
          else
          {
            // The coarse mesh-part does not have a topology, stick to simple refinement
            Intern::SimpleTargetRefineWrapper<ShapeType>::refine(
              target_set_holder, num_entities_parent, coarse_target_set_holder);
          }
        }
      }; // struct TargetSetRefineParentWrapper<StructuredMesh<...>>
    } // namespace Intern
 
#ifdef FEAT_EICKT
    extern template class MeshPart<ConformalMesh<Shape::Simplex<2>, 2, Real>>;
    extern template class MeshPart<ConformalMesh<Shape::Simplex<3>, 3, Real>>;
    extern template class MeshPart<ConformalMesh<Shape::Hypercube<2>, 2, Real>>;
    extern template class MeshPart<ConformalMesh<Shape::Hypercube<3>, 3, Real>>;
 
    extern template class StandardRefinery<MeshPart<ConformalMesh<Shape::Simplex<2>, 2, Real>>>;
    extern template class StandardRefinery<MeshPart<ConformalMesh<Shape::Simplex<3>, 3, Real>>>;
    extern template class StandardRefinery<MeshPart<ConformalMesh<Shape::Hypercube<2>, 2, Real>>>;
    extern template class StandardRefinery<MeshPart<ConformalMesh<Shape::Hypercube<3>, 3, Real>>>;
#endif // FEAT_EICKT
  } // namespace Geometry
} // namespace FEAT