adaptive_parti_domain_control.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
 
#include <kernel/base_header.hpp>
#include <kernel/geometry/adaptive_mesh.hpp>
#include <kernel/geometry/adaptive_mesh_layer.hpp>
#include <kernel/geometry/mesh_node.hpp>
#include <kernel/geometry/mesh_permutation.hpp>
#include <kernel/geometry/subdivision_levels.hpp>
#include <kernel/util/dist.hpp>
 
#include <kernel/geometry/boundary_factory.hpp>
#include <kernel/geometry/common_factories.hpp>
#include <kernel/geometry/mesh_file_reader.hpp>
 
#include <control/domain/parti_domain_control.hpp>
#include <memory>
#include <optional>
#include <utility>
 
namespace FEAT::Control::Domain
{
  template<typename DomainLevel_, typename TemplateSet_>
  class AdaptiveLevelWrapper : public DomainLevel_
  {
  public:
    using BaseClass = DomainLevel_;
    using typename BaseClass::MeshType;
    using typename BaseClass::MeshNodeType;
    using AdaptiveMeshType = Geometry::AdaptiveMesh<TemplateSet_, typename MeshType::ShapeType>;
    using AdaptiveMeshLayerType = Geometry::AdaptiveMeshLayer<AdaptiveMeshType>;
 
    std::optional<AdaptiveMeshLayerType> mesh_layer;
 
    explicit AdaptiveLevelWrapper(int _lvl_idx, std::unique_ptr<MeshNodeType> node) :
      BaseClass(_lvl_idx, std::move(node)),
      mesh_layer(std::nullopt)
    {
    }
 
    explicit AdaptiveLevelWrapper(int _lvl_idx, std::unique_ptr<MeshNodeType> node, AdaptiveMeshLayerType&& layer) :
      BaseClass(_lvl_idx, std::move(node)),
      mesh_layer(std::make_optional<AdaptiveMeshLayerType>(std::move(layer)))
    {
    }
 
    bool is_partial_level()
    {
      return bool(mesh_layer);
    }
  }; // class AdaptiveLevelWrapper<...>
 
  template<typename DomainLevel_, typename TemplateSet_>
  class AdaptivePartiDomainControl : public PartiDomainControl<DomainLevel_>
  {
  public:
    using BaseClass = Control::Domain::PartiDomainControl<DomainLevel_>;
    using typename BaseClass::LevelType;
    using typename BaseClass::LayerType;
    using typename BaseClass::MeshType;
    using typename BaseClass::AtlasType;
    using typename BaseClass::MeshNodeType;
    using typename BaseClass::MeshPartType;
    using typename BaseClass::Ancestor;
    using typename BaseClass::WeightType;
 
    using AdaptiveMeshType = Geometry::AdaptiveMesh<TemplateSet_, typename MeshType::ShapeType>;
 
  protected:
    std::function<void(Geometry::SubdivisionLevels&, const MeshType&)> _refinement_strategy;
 
    std::shared_ptr<AdaptiveMeshType> _adaptive_mesh;
 
  public:
    explicit AdaptivePartiDomainControl(
      const Dist::Comm& comm_,
      bool support_multi_layered,
      std::function<void(Geometry::SubdivisionLevels&, const MeshType&)> refinement_strategy) :
      BaseClass(comm_, support_multi_layered),
      _refinement_strategy(std::move(refinement_strategy))
    {
      this->_allow_parti_2level = false;
    }
 
    virtual ~AdaptivePartiDomainControl() = default;
 
    const AdaptiveMeshType& get_adaptive_mesh() const
    {
      return *_adaptive_mesh;
    }
 
  protected:
 
    std::vector<WeightType> _compute_weights(Ancestor& DOXY(ancestor), const MeshNodeType& base_mesh_node) override
    {
      static constexpr int dim = MeshType::ShapeType::dimension;
      const WeightType average_elements_per_marking =
        WeightType(TemplateSet_::template average_elements_per_marking<dim>());
 
      const MeshType& mesh = *base_mesh_node.get_mesh();
      const auto& v_at_c = mesh.template get_index_set<dim, 0>();
 
      Geometry::SubdivisionLevels sdls(mesh.get_num_vertices());
      std::vector<WeightType> weights(mesh.get_num_elements());
 
      // Get subdivision levels for mesh
      _refinement_strategy(sdls, mesh);
 
      // Determine weights from subdivision levels
      for(Index cell(0); cell < weights.size(); cell++)
      {
        std::uint64_t markings(0);
        for(auto i = v_at_c.image_begin(cell); i != v_at_c.image_end(cell); ++i)
        {
          markings += sdls[*i];
        }
        // Minimum weight is 1
        weights[cell] = Math::max(WeightType(1), WeightType(markings) * average_elements_per_marking);
      }
 
      return weights;
    }
 
    void _create_single_process(std::unique_ptr<MeshNodeType> base_mesh_node) override
    {
      // Perform regular refinement of conformal mesh
      BaseClass::_create_single_process(std::move(base_mesh_node));
      refine_partially();
    }
 
    // Note: all following member functions are only required for parallel builds,
    // so we enclose them in the following #if-block to reduce compile times.
 
#if defined(FEAT_HAVE_MPI) || defined(DOXYGEN)
 
    void _create_single_layered(std::unique_ptr<MeshNodeType> base_mesh_node) override
    {
      BaseClass::_create_single_layered(std::move(base_mesh_node));
      refine_partially();
    }
 
    void _create_multi_layered(std::unique_ptr<MeshNodeType> base_mesh_node) override
    {
      BaseClass::_create_multi_layered(std::move(base_mesh_node));
      refine_partially();
    }
 
#endif // defined(FEAT_HAVE_MPI) || defined(DOXYGEN)
    void refine_partially()
    {
      const MeshNodeType& mesh_node = *this->_layer_levels.at(0).front()->get_mesh_node();
      Geometry::SubdivisionLevels sdls(mesh_node.get_mesh()->get_num_vertices());
 
      // Get requested subdivision levels
      _refinement_strategy(sdls, *mesh_node.get_mesh());
 
      // Sync subdivision levels across processes
      sync_subdivision_levels(sdls, mesh_node);
 
      // Sync maximum subdivision level across ranks
      // This value will determine how many copy-layers we need to add
      const std::uint64_t sdl_max_local = sdls.maximum_level();
      std::uint64_t sdl_max;
      this->comm().allreduce(&sdl_max_local, &sdl_max, std::size_t(1), Dist::op_max);
 
      // Create adaptive mesh
      _adaptive_mesh = std::make_shared<AdaptiveMeshType>(*mesh_node.get_mesh());
      _adaptive_mesh->adapt(sdls, AdaptiveMeshType::ImportBehaviour::All);
 
      int lvl = this->_chosen_levels.front().first;
      for(Index i(1); i <= sdl_max; i++)
      {
        // Either push the next layer of the adaptive mesh or a copy of the last layer.
        // We push a copy to ensure we correctly participate in global operations on all levels.
        // TODO: Use a more efficient method, rather than copying the previous level
        Geometry::Layer layer_idx{Math::min(i, _adaptive_mesh->num_layers() - 1)};
        Geometry::AdaptiveMeshLayer<AdaptiveMeshType> layer(_adaptive_mesh, layer_idx);
        this->push_level_front(
          0,
          std::make_shared<LevelType>(
            lvl++,
            project_mesh_node(*_adaptive_mesh, layer_idx, mesh_node),
            std::move(layer)));
      }
 
      this->_chosen_levels.emplace_front(lvl, this->_comm.size());
    }
 
    std::unique_ptr<MeshNodeType> project_mesh_node(
      const AdaptiveMeshType& mesh,
      Geometry::Layer l,
      const MeshNodeType& root_node,
      Geometry::AdaptMode adapt_mode = Geometry::AdaptMode::chart)
    {
 
      auto cmesh = std::make_unique<MeshType>(mesh.to_conformal_mesh(l));
      auto result = std::make_unique<MeshNodeType>(std::move(cmesh));
 
      // Projects mesh parts
      for(const auto& name : root_node.get_mesh_part_names())
      {
        const auto& mesh_part_node = root_node.find_mesh_part_node(name);
        const auto* mesh_part = mesh_part_node->get_mesh();
 
        std::unique_ptr<MeshPartType> new_mesh_part;
        if(mesh_part)
        {
          new_mesh_part = std::make_unique<MeshPartType>(mesh.template project_meshpart<MeshType>(l, *mesh_part));
        }
        else
        {
          new_mesh_part = std::unique_ptr<MeshPartType>(nullptr);
        }
        auto new_mesh_part_node = std::make_unique<Geometry::MeshPartNode<MeshType>>(std::move(new_mesh_part));
 
        result->add_mesh_part_node(
          name,
          std::move(new_mesh_part_node),
          mesh_part_node->find_mesh_part_chart_name(name),
          mesh_part_node->find_mesh_part_chart(name));
      }
 
      // refine our halos
      for(const auto& [rank, part] : root_node.get_halo_map())
      {
        if(part)
        {
          XASSERT(bool(part));
          auto new_halo = std::make_unique<MeshPartType>(mesh.template project_meshpart<MeshType>(l, *part));
          result->add_halo(rank, std::move(new_halo));
        }
        else
        {
          result->add_halo(rank, nullptr);
        }
      }
 
      // refine our patch mesh-parts
      for(const auto& [rank, patch] : root_node.get_patch_map())
      {
        if(patch)
        {
          XASSERT(bool(patch));
          auto new_patch = std::make_unique<MeshPartType>(mesh.template project_meshpart<MeshType>(l, *patch));
          result->add_patch(rank, std::move(new_patch));
        }
        else
        {
          result->add_patch(rank, nullptr);
        }
      }
 
      // adapt by chart?
      if((adapt_mode & Geometry::AdaptMode::chart) != Geometry::AdaptMode::none)
      {
        result->adapt(true);
      }
 
      // TODO: Adapt dual
 
      return result;
    }
 
    void sync_subdivision_levels(Geometry::SubdivisionLevels& sdls, const MeshNodeType& mesh_node)
    {
      XASSERT(mesh_node.get_mesh()->get_num_vertices() == sdls.size());
 
      std::unordered_map<int, Dist::RequestVector> sends;
      std::unordered_map<int, Dist::RequestVector> receives;
 
      Geometry::SubdivisionLevels received_sdls(sdls.size(), 0);
 
      // Post send and receives
      for(const auto& [rank, halo] : mesh_node.get_halo_map())
      {
        const auto& vertices = halo->template get_target_set<0>();
 
        for(Index i(0); i < vertices.get_num_entities(); i++)
        {
          sends[rank].push_back(this->_comm.isend(&sdls[vertices[i]], 1, rank));
          receives[rank].push_back(this->_comm.irecv(&received_sdls[vertices[i]], 1, rank));
        }
      }
 
      // Wait for all receives and map sdls
      for(const auto& [rank, halo] : mesh_node.get_halo_map())
      {
        const auto& vertices = halo->template get_target_set<0>();
        for(Index idx(0); receives[rank].wait_any(idx);)
        {
          const Index vidx = vertices[idx];
          // Subdivision level is maximum of all processes values
          sdls[vidx] = Math::max(sdls[vidx], received_sdls[vidx]);
        }
      }
 
      // Wait for all sends to finish
      for(const auto& [rank, halo] : mesh_node.get_halo_map())
      {
        sends[rank].wait_all();
      }
    }
  };
} // namespace FEAT::Control::Domain