boundary_factory.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/conformal_mesh.hpp>
#include <kernel/geometry/mesh_part.hpp>
#include <kernel/geometry/intern/boundary_computer.hpp>
#include <kernel/util/dist.hpp>
 
namespace FEAT
{
  namespace Geometry
  {
    template<typename ParentMesh_>
    class BoundaryFactory :
      public Factory<MeshPart<ParentMesh_>>
    {
    public:
      typedef Factory<MeshPart<ParentMesh_>> BaseClass;
      typedef ParentMesh_ ParentMeshType;
      typedef MeshPart<ParentMeshType> MeshType;
      typedef typename MeshType::ShapeType ShapeType;
      typedef typename MeshType::TargetSetHolderType TargetSetHolderType;
      static constexpr int shape_dim = ShapeType::dimension;
 
    private:
      const ParentMeshType& _mesh_in;
      Intern::BoundaryFaceComputer<ShapeType> _face_computer;
      std::array<std::vector<int>, std::size_t(shape_dim)> _bnd_masks;
      std::array<std::vector<Index>, std::size_t(shape_dim)> _face_idx;
 
    public:
      explicit BoundaryFactory(const ParentMeshType& mesh_in) :
        _mesh_in(mesh_in)
      {
        _face_computer.compute_all(mesh_in.get_index_set_holder(), _bnd_masks, _face_idx);
      }
 
      virtual Index get_num_entities(int dim) override
      {
        if(dim < shape_dim)
          return Index(_face_idx.at(std::size_t(dim)).size());
        else
          return Index(0);
      }
 
      virtual void fill_attribute_sets(typename BaseClass::AttributeSetContainer&) override
      {
      }
 
      virtual void fill_index_sets(std::unique_ptr<typename BaseClass::IndexSetHolderType>&) override
      {
      }
 
      virtual void fill_target_sets(TargetSetHolderType& target_set_holder) override
      {
        _face_computer.fill_target_sets(target_set_holder, _face_idx);
      }
 
    }; // class BoundaryFactory<...>
 
    template<typename Mesh_>
    MeshPart<Mesh_> make_boundary_meshpart(const Mesh_& mesh)
    {
      BoundaryFactory<Mesh_> factory(mesh);
      return factory.make();
    }
 
    template<typename Mesh_>
    std::unique_ptr<MeshPart<Mesh_>> make_unique_boundary_meshpart(const Mesh_& mesh)
    {
      BoundaryFactory<Mesh_> factory(mesh);
      return factory.make_unique();
    }
 
 
    template<typename ParentMesh_>
    class MaskedBoundaryFactory :
      public Factory<MeshPart<ParentMesh_>>
    {
    public:
      typedef Factory<MeshPart<ParentMesh_>> BaseClass;
      typedef ParentMesh_ ParentMeshType;
      typedef MeshPart<ParentMeshType> MeshType;
      typedef typename MeshType::ShapeType ShapeType;
      typedef typename MeshType::TargetSetHolderType TargetSetHolderType;
      static constexpr int shape_dim = ShapeType::dimension;
      static constexpr int facet_dim = shape_dim - 1;
 
    private:
      const ParentMeshType& _mesh_in;
      Intern::BoundaryFaceComputer<ShapeType> _face_computer;
      std::vector<int> _facet_mask;
      std::array<std::vector<int>, std::size_t(shape_dim)> _bnd_masks;
      std::array<std::vector<Index>, std::size_t(shape_dim)> _face_idx;
 
    public:
      explicit MaskedBoundaryFactory(const ParentMeshType& mesh_in) :
        _mesh_in(mesh_in),
        _facet_mask(mesh_in.get_num_entities(facet_dim), 0)
      {
      }
 
      // Adds a single facet to the boundary facet mask.
      void add_mask_facet(Index facet_idx)
      {
        _facet_mask.at(facet_idx) = 1;
      }
 
      void add_mask_meshpart(const MeshPart<ParentMeshType>& part)
      {
        const TargetSet& trg_set = part.template get_target_set<facet_dim>();
        for(Index i(0); i < trg_set.get_num_entities(); ++i)
          _facet_mask.at(trg_set[i]) = 1;
      }
 
      void compile()
      {
        _face_computer.compute_masks(_mesh_in.get_index_set_holder(), _bnd_masks, _face_idx, _facet_mask);
        _face_computer.compute_faces(_bnd_masks, _face_idx);
      }
 
      virtual Index get_num_entities(int dim) override
      {
        if(dim < shape_dim)
          return Index(_face_idx.at(std::size_t(dim)).size());
        else
          return Index(0);
      }
 
      virtual void fill_attribute_sets(typename BaseClass::AttributeSetContainer&) override
      {
      }
 
      virtual void fill_index_sets(std::unique_ptr<typename BaseClass::IndexSetHolderType>&) override
      {
      }
 
      virtual void fill_target_sets(TargetSetHolderType& target_set_holder) override
      {
        _face_computer.fill_target_sets(target_set_holder, _face_idx);
      }
    }; // class MaskedBoundaryFactory<...>
 
    template<typename ParentMesh_>
    class GlobalMaskedBoundaryFactory :
      public Factory<MeshPart<ParentMesh_>>
    {
    public:
      typedef Factory<MeshPart<ParentMesh_>> BaseClass;
      typedef ParentMesh_ ParentMeshType;
      typedef MeshPart<ParentMeshType> MeshType;
      typedef typename MeshType::ShapeType ShapeType;
      typedef typename MeshType::TargetSetHolderType TargetSetHolderType;
      static constexpr int shape_dim = ShapeType::dimension;
      static constexpr int facet_dim = shape_dim - 1;
 
    private:
      const ParentMeshType& _mesh_in;
      Intern::BoundaryFaceComputer<ShapeType> _face_computer;
      std::vector<int> _facet_mask;
      std::array<std::vector<int>, std::size_t(shape_dim)> _bnd_masks;
      std::array<std::vector<Index>, std::size_t(shape_dim)> _face_idx;
      std::vector<std::pair<int, const MeshType*>> _halos;
 
    public:
      explicit GlobalMaskedBoundaryFactory(const ParentMeshType& mesh_in) :
        _mesh_in(mesh_in),
        _facet_mask(mesh_in.get_num_entities(facet_dim), 0)
      {
      }
 
      // Adds a single facet to the mask.
      void add_mask_facet(Index facet_idx)
      {
        _facet_mask.at(facet_idx) = 1;
      }
 
      void add_mask_meshpart(const MeshPart<ParentMeshType>& part)
      {
        const TargetSet& trg_set = part.template get_target_set<facet_dim>();
        for(Index i(0); i < trg_set.get_num_entities(); ++i)
          _facet_mask.at(trg_set[i]) = 1;
      }
 
      void add_halo(int rank, const MeshPart<ParentMeshType>& halo)
      {
        // add the halo as a mask meshpart
        add_mask_meshpart(halo);
        // and save the halo
        //_halos.emplace(rank, &halo);
        _halos.push_back(std::make_pair(rank, &halo));
      }
 
      void compile(const Dist::Comm& comm)
      {
        _face_computer.compute_masks(_mesh_in.get_index_set_holder(), _bnd_masks, _face_idx, _facet_mask);
        _sync_bnd_masks(comm);
        _face_computer.compute_faces(_bnd_masks, _face_idx);
      }
 
      virtual Index get_num_entities(int dim) override
      {
        if(dim < shape_dim)
          return Index(_face_idx.at(std::size_t(dim)).size());
        else
          return Index(0);
      }
 
      virtual void fill_attribute_sets(typename BaseClass::AttributeSetContainer&) override
      {
      }
 
      virtual void fill_index_sets(std::unique_ptr<typename BaseClass::IndexSetHolderType>&) override
      {
      }
 
      virtual void fill_target_sets(TargetSetHolderType& target_set_holder) override
      {
        _face_computer.fill_target_sets(target_set_holder, _face_idx);
      }
 
    protected:
      void _sync_bnd_masks(const Dist::Comm& comm)
      {
        if(_halos.empty())
          return;
 
        // allocate halo buffers
        std::vector<std::vector<int>> halo_send_bufs(_halos.size()), halo_recv_bufs(_halos.size());
 
        // create send buffers and allocate receive buffers
        for(std::size_t i(0); i < _halos.size(); ++i)
        {
          // count total number of entries in halo
          Index count = 0u;
          for(int dim = 0; dim < shape_dim; ++dim) // we don't need to sync dim=shape_dim here
            count += _halos[i].second->get_num_entities(dim);
 
          halo_recv_bufs[i].resize(count);
          halo_send_bufs[i].resize(count); // build_halo_buffer uses push_back
          _face_computer.build_halo_buffer(halo_send_bufs[i], _bnd_masks, _halos[i].second->get_target_set_holder());
        }
 
        // post receive requests
        Dist::RequestVector recv_reqs(_halos.size());
        for(std::size_t i(0); i < _halos.size(); ++i)
          recv_reqs.push_back(comm.irecv(halo_recv_bufs[i].data(), halo_recv_bufs[i].size(), _halos[i].first));
 
        // post send requests
        Dist::RequestVector send_reqs(_halos.size());
        for(std::size_t i(0); i < _halos.size(); ++i)
          send_reqs.push_back(comm.isend(halo_send_bufs[i].data(), halo_send_bufs[i].size(), _halos[i].first));
 
        // wait for all receives to finish
        recv_reqs.wait_all();
 
        // process receives halo buffers
        for(std::size_t i(0); i < _halos.size(); ++i)
        {
          _face_computer.mask_halo_buffer(_bnd_masks, halo_recv_bufs[i], _halos[i].second->get_target_set_holder());
        }
 
        // wait for all sends to finish
        send_reqs.wait_all();
      }
    }; // class GlobalMaskedBoundaryFactory<...>
 
    template<typename MeshNode_>
    MeshPart<typename MeshNode_::MeshType> make_boundary_meshpart(const MeshNode_& node, const Dist::Comm& comm)
    {
      GlobalMaskedBoundaryFactory<typename MeshNode_::MeshType> boundary_factory(*node.get_mesh());
 
      for(const auto& pair : node.get_halo_map())
      {
        boundary_factory.add_halo(pair.first, *pair.second);
      }
      boundary_factory.compile(comm);
      return boundary_factory.make();
    }
 
    template<typename MeshNode_>
    std::unique_ptr<MeshPart<typename MeshNode_::MeshType>> make_unique_boundary_meshpart(const MeshNode_& node, const Dist::Comm& comm)
    {
      GlobalMaskedBoundaryFactory<typename MeshNode_::MeshType> boundary_factory(*node.get_mesh());
 
      for(const auto& pair : node.get_halo_map())
      {
        boundary_factory.add_halo(pair.first, *pair.second);
      }
      boundary_factory.compile(comm);
      return boundary_factory.make_unique();
    }
  } // namespace Geometry
} // namespace FEAT