patch_meshpart_splitter.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/mesh_part.hpp>
 
// includes, system
#include <map>
 
namespace FEAT
{
  namespace Geometry
  {
    // Forward declarations
    class PatchPartMap;
 
    template<typename, int>
    class PatchPartMapHolder;
 
    template<typename Mesh_>
    class PatchMeshPartSplitter;
 
    template<typename Shape_, int num_coords_, typename Coord_>
    class PatchMeshPartSplitter<Geometry::ConformalMesh<Shape_, num_coords_, Coord_>> :
      public Factory<MeshPart<Geometry::ConformalMesh<Shape_, num_coords_, Coord_>>>
    {
    public:
      typedef Shape_ ShapeType;
      static constexpr int shape_dim = ShapeType::dimension;
      typedef Geometry::ConformalMesh<Shape_, num_coords_, Coord_> MeshType;
      typedef Geometry::MeshPart<MeshType> MeshPartType;
      typedef typename MeshType::VertexSetType::CoordType AttributeDataType;
      typedef typename MeshPartType::AttributeSetContainer AttributeSetContainer;
      typedef typename MeshPartType::IndexSetHolderType IndexSetHolderType;
      typedef typename MeshPartType::TargetSetHolderType TargetSetHolderType;
 
    protected:
      const MeshType& _base_mesh;
      const MeshPartType& _patch_mesh_part;
 
      const AttributeSetContainer* _cur_part_attribute_container;
      const IndexSetHolder<ShapeType>* _cur_part_topology;
      PatchPartMapHolder<ShapeType, ShapeType::dimension> _part_holder;
 
    public:
      explicit PatchMeshPartSplitter(const MeshType& base_mesh, const MeshPartType& patch_mesh_part) :
        _base_mesh(base_mesh),
        _patch_mesh_part(patch_mesh_part),
        _cur_part_attribute_container(nullptr),
        _cur_part_topology(nullptr),
        _part_holder(patch_mesh_part.get_target_set_holder())
        {
        }
 
      virtual ~PatchMeshPartSplitter()
      {
      }
 
      bool build(const MeshPartType& mesh_part)
      {
        _cur_part_topology = mesh_part.get_topology();
        _cur_part_attribute_container = &(mesh_part.get_mesh_attributes());
        return _part_holder.build(mesh_part.get_target_set_holder(), mesh_part.get_topology());
      }
 
      /* *************************************************************************************** */
      /* F A C T O R Y   I N T E R F A C E   I M P L E M E N T A T I O N                         */
      /* *************************************************************************************** */
 
      virtual Index get_num_entities(int dim) override
      {
        return _part_holder.get_num_entities(dim);
      }
 
      virtual void fill_attribute_sets(AttributeSetContainer& attribute_container) override
      {
        _part_holder.fill_attribute_sets(attribute_container, *_cur_part_attribute_container);
      }
 
      virtual void fill_index_sets(std::unique_ptr<IndexSetHolderType>& index_set_holder) override
      {
        XASSERT(index_set_holder == nullptr);
 
        // If the base MeshPart has a topology, create one for the patch MeshPart, too
        if(_part_holder.has_topology())
        {
          // Get num_entities for all shape dimensions
          Index num_entities[shape_dim+1];
          for(int i(0); i < shape_dim+1; ++i)
            num_entities[i] = get_num_entities(i);
 
          // Create this MeshPart's topology
          index_set_holder.reset(new IndexSetHolderType(num_entities));
          // This creates vertex-at-shape information for all shapes
          _part_holder.fill_index_sets(*index_set_holder, *_cur_part_topology);
          // Build the redundant index sets
          RedundantIndexSetBuilder<ShapeType>::compute(*index_set_holder);
        }
      }
 
      virtual void fill_target_sets(TargetSetHolderType& target_set_holder) override
      {
        _part_holder.fill_target_sets(target_set_holder);
      }
    }; // class PatchMeshPartSplitter<ConformalMesh<...>>
 
    class PatchPartMap
    {
    private:
      std::map<Index,Index> _idx_map;
      std::map<Index,Index> _io_map;
      std::vector<Index> _indices;
 
    public:
      explicit PatchPartMap(const TargetSet& target_set) :
        _idx_map(),
        _io_map(),
        _indices()
      {
        // loop over all indices
        for(Index i(0); i < target_set.get_num_entities(); ++i)
        {
          // This maps the BaseMesh entity j = target_set[i] to the patch MeshPart entity i
          XASSERTM(_idx_map.emplace(target_set[i], i).second, "internal error");
        }
      }
 
      bool build(const TargetSet& target_in)
      {
        // Very important: Clear old information
        _indices.clear();
        _io_map.clear();
 
        // We need information for every entity in the BaseMesh MeshPart
        for(Index i(0); i < target_in.get_num_entities(); ++i)
        {
          // Check if the BaseMesh entity referenced by the BaseMesh MeshPart is in the PatchMeshPart
          auto it = _idx_map.find(target_in[i]);
          if(it != _idx_map.end())
          {
            // We create new entries in the PatchMeshPart MeshPart here, so the entity i in the BaseMesh MeshPart
            // refers to the newly created entity in the PatchMeshPart MeshPart
            _io_map.emplace(i, Index(_indices.size()));
            // Add the new entity
            _indices.push_back(it->second);
          }
        }
        return !_indices.empty();
      }
 
      Index size() const
      {
        return Index(_indices.size());
      }
 
      template<typename Attribute_>
      void fill_attribute_set(
        std::map<String, std::unique_ptr<Attribute_>>& attribute_container_out,
        const std::map<String, std::unique_ptr<Attribute_>>& attribute_container_in) const
      {
        for(auto it(attribute_container_in.begin()); it != attribute_container_in.end(); ++it)
        {
          // Create new attribute for the PatchMeshPart MeshPart
          std::unique_ptr<Attribute_> new_attribute(new Attribute_(
            Index(_indices.size()), it->second->get_dimension()));
 
          const int dim = it->second->get_dimension();
 
          // Emplace all elements referred to by the _io_map into new_attribute
          for(auto jt: _io_map)
          {
            for(int k(0); k < dim; ++k)
              new_attribute->operator()(jt.second, k) = it->second->operator()(jt.first, k);
          }
 
          // Push the new attribute to set_out
          XASSERTM(attribute_container_out.insert(std::make_pair(it->first, std::move(new_attribute))).second, "Error inserting new AttributeSet");
        }
      }
 
      void fill_target_set(TargetSet& target_out) const
      {
        XASSERT(target_out.get_num_entities() == size());
        for(Index i(0); i < target_out.get_num_entities(); ++i)
        {
          target_out[i] = _indices[i];
        }
      }
 
      template<typename IndexSetType_>
      void fill_index_set(IndexSetType_& index_set_out, const IndexSetType_& index_set_in,
        const std::map<Index, Index>& vertex_map) const
      {
        // We need to check all shape entities in the MeshPart referring to the BaseMesh
        for(Index k(0); k < index_set_in.get_num_entities(); ++k)
        {
          // Check if the BaseMesh entity referred to by the BaseMesh MeshPart is also present in the PatchMeshPart
          auto kt = _io_map.find(k);
          if(kt != _io_map.end())
          {
            // For all vertices at that entity, we need to find them in the PatchMeshPart's MeshPart
            for(int j(0); j < index_set_in.get_num_indices(); ++j)
            {
              // Fetch the BaseMesh MeshPart's vertex number
              Index i_in(index_set_in(k,j));
              // That vertex is definitely in the PatchMeshPart's MeshPart, but we need to catch any errors here
              auto jt = vertex_map.find(i_in);
              if(jt != vertex_map.end())
                // This is the correct number
                index_set_out(kt->second, j) = jt->second;
              else
                XABORTM("Vertex "+stringify(i_in)+" missing in MeshPart topology!");
            }
          }
        }
      }
 
      const std::map<Index, Index>& get_io_map() const
      {
        return _io_map;
      }
    }; // class PatchPartMap
 
    template<typename Shape_, int dim_ = Shape_::dimension>
    class PatchPartMapHolder :
      public PatchPartMapHolder<Shape_, dim_ - 1>
    {
    public:
      typedef PatchPartMapHolder<Shape_, dim_ - 1> BaseClass;
 
    private:
      PatchPartMap _patch_map;
      bool _has_topology;
 
    public:
      explicit PatchPartMapHolder(const TargetSetHolder<Shape_>& tsh) :
        BaseClass(tsh),
        _patch_map(tsh.template get_target_set<dim_>()),
        _has_topology(false)
      {
      }
 
      bool build(const TargetSetHolder<Shape_>& tsh, const IndexSetHolder<Shape_>* ish)
      {
        ish == nullptr ? _has_topology = false : _has_topology = true;
 
        // Call the BaseClass version (meaning of lower shape dimension) first
        bool b1 = BaseClass::build(tsh, ish);
        bool b2 =  _patch_map.build(tsh.template get_target_set<dim_>());
 
        return  b1 || b2;
      }
 
      Index get_num_entities(int dim) const
      {
        XASSERT(dim <= dim_);
        if(dim == dim_)
          return _patch_map.size();
        else
          return BaseClass::get_num_entities(dim);
      }
 
      bool has_topology() const
      {
        return _has_topology;
      }
 
      template<typename AttributeSetContainer_>
      void fill_attribute_sets(
        AttributeSetContainer_& attribute_container_out,
        const AttributeSetContainer_& attribute_container_in) const
      {
        // Fill attribute sets of one shape dimension lower first
        BaseClass::fill_attribute_sets(attribute_container_out, attribute_container_in);
      }
 
      void fill_target_sets(TargetSetHolder<Shape_>& tsh) const
      {
        BaseClass::fill_target_sets(tsh);
        _patch_map.fill_target_set(tsh.template get_target_set<dim_>());
      }
 
      template<typename IndexSetHolder_>
      void fill_index_sets(IndexSetHolder_& ish, const IndexSetHolder_& ish_in) const
      {
        BaseClass::fill_index_sets(ish, ish_in);
        _patch_map.fill_index_set(ish.template get_index_set<dim_,0>(), ish_in.template get_index_set<dim_,0>(),
        get_vertex_map());
      }
 
      const std::map<Index, Index>& get_vertex_map() const
      {
        return BaseClass::get_vertex_map();
      }
    }; // class PatchPartMapHolder<...>
 
 
    template<typename Shape_>
    class PatchPartMapHolder<Shape_, 0>
    {
    private:
      PatchPartMap _patch_map;
 
    public:
      explicit PatchPartMapHolder(const TargetSetHolder<Shape_>& tsh) :
        _patch_map(tsh.template get_target_set<0>())
        {
        }
 
      bool build(const TargetSetHolder<Shape_>& tsh, const IndexSetHolder<Shape_>* DOXY(ish))
      {
        return _patch_map.build(tsh.template get_target_set<0>());
      }
 
      Index get_num_entities(int dim) const
      {
        XASSERT(dim == 0);
        return _patch_map.size();
      }
 
      template<typename AttributeSetContainer_>
      void fill_attribute_sets(
        AttributeSetContainer_& attribute_container_out,
        const AttributeSetContainer_& attribute_container_in) const
      {
        // Fill attribute set to shape dimension 0
        _patch_map.fill_attribute_set( attribute_container_out, attribute_container_in);
      }
 
      template<typename IndexSetHolder_>
      void fill_index_sets(IndexSetHolder_& DOXY(ish), const IndexSetHolder_& DOXY(ish_in)) const
      {
      }
 
      void fill_target_sets(TargetSetHolder<Shape_>& tsh) const
      {
        _patch_map.fill_target_set(tsh.template get_target_set<0>());
      }
 
 
      const std::map<Index, Index>& get_vertex_map() const
      {
        return _patch_map.get_io_map();
      }
    }; // class PatchPartMapHolder<Shape_, 0>
 
  } // namespace Geometry
} // namespace FEAT