Program Listing for File mesh.hxx

Return to documentation for file (mesh_handle/mesh.hxx)

/*
  This file is part of t8code.
  t8code is a C library to manage a collection (a forest) of multiple
  connected adaptive space-trees of general element classes in parallel.

  Copyright (C) 2026 the developers

  t8code is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2 of the License, or
  (at your option) any later version.

  t8code is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with t8code; if not, write to the Free Software Foundation, Inc.,
  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/

#pragma once

#include <t8.h>
#include "element.hxx"
#include "competence_pack.hxx"
#include "internal/adapt.hxx"
#include "concepts.hxx"
#include "data_handler.hxx"
#include <t8_forest/t8_forest_balance.h>
#include <t8_forest/t8_forest_types.h>
#include <t8_forest/t8_forest_general.h>
#include <t8_forest/t8_forest_ghost.h>
#include <vector>
#include <functional>
#include <memory>
#include <span>

namespace t8_mesh_handle
{

template <typename TType>
concept ElementCompetencePack = requires { typename TType::is_element_competence_pack; };
template <typename TType>
concept MeshCompetencePack = requires { typename TType::is_mesh_competence_pack; };

template <ElementCompetencePack TElementCompetencePack = element_competence_pack<>,
          MeshCompetencePack TMeshCompetencePack = mesh_competence_pack<>>
class mesh: public TMeshCompetencePack::template apply<mesh<TElementCompetencePack, TMeshCompetencePack>> {
 public:
  using SelfType = mesh<TElementCompetencePack, TMeshCompetencePack>;
  using element_class = typename TElementCompetencePack::template apply<
    SelfType, element>;
  friend element_class;
  using mesh_tag = void;
  using mesh_const_iterator =
    typename std::vector<element_class>::const_iterator;
  using mesh_iterator =
    typename std::vector<element_class>::iterator;
  friend struct element_data_element_competence<element_class>;
  using adapt_callback_type = std::function<int (const SelfType& mesh, std::span<const element_class> elements)>;

  template <typename TUserDataType>
  using adapt_callback_type_with_userdata
    = std::function<int (const SelfType& mesh, std::span<const element_class> elements, TUserDataType user_data)>;

  mesh (t8_forest_t forest): m_forest (forest)
  {
    T8_ASSERT (t8_forest_is_committed (m_forest));
    update_elements ();
  }

  ~mesh ()
  {
    t8_forest_unref (&m_forest);
  }

  // --- Getter for mesh related information. ---
  t8_locidx_t
  get_num_local_elements () const
  {
    return t8_forest_get_local_num_leaf_elements (m_forest);
  }

  t8_locidx_t
  get_num_ghosts () const
  {
    return t8_forest_get_num_ghosts (m_forest);
  }

  int
  get_dimension () const
  {
    return t8_forest_get_dimension (m_forest);
  }

  t8_forest_t
  get_forest () const
  {
    return m_forest;
  }

  bool
  is_balanced ()
  {
    return t8_forest_is_balanced (m_forest);
  }

  // --- Methods to access elements. ---
  mesh_const_iterator
  cbegin () const
  {
    return m_elements.cbegin ();
  }

  mesh_const_iterator
  cend () const
  {
    return m_elements.cend ();
  }

  mesh_iterator
  begin ()
  {
    return m_elements.begin ();
  }

  mesh_iterator
  end ()
  {
    return m_elements.end ();
  }

  mesh_const_iterator
  begin () const
  {
    return this->cbegin ();
  }

  mesh_const_iterator
  end () const
  {
    return this->cend ();
  }

  const element_class&
  operator[] (t8_locidx_t local_index) const
  {
    T8_ASSERT (0 <= local_index && local_index < get_num_local_elements () + get_num_ghosts ());
    if (local_index < get_num_local_elements ()) {
      return m_elements[local_index];
    }
    else {
      return m_ghosts[local_index - get_num_local_elements ()];
    }
  }

  element_class&
  operator[] (t8_locidx_t local_index)
  {
    return const_cast<element_class&> (static_cast<const mesh*> (this)->operator[] (local_index));
  }

  // --- Methods to change the mesh, e.g. adapt, partition, balance, ... ---
  template <typename TUserDataType>
  static adapt_callback_type
  mesh_adapt_callback_wrapper (adapt_callback_type_with_userdata<TUserDataType> adapt_callback_with_userdata,
                               const TUserDataType& user_data)
  {
    return [=] (const SelfType& mesh, std::span<const element_class> elements) {
      return adapt_callback_with_userdata (mesh, elements, user_data);
    };
  }

  void
  set_adapt (adapt_callback_type adapt_callback)
  {
    if (!m_uncommitted_forest.has_value ()) {
      m_uncommitted_forest.emplace ();
      t8_forest_init (&*m_uncommitted_forest);
    }
    // Create and register adaptation context holding the mesh handle and the user defined callback.
    detail::adapt_registry::register_context (
      m_forest, std::make_unique<detail::mesh_adapt_context<SelfType>> (*this, std::move (adapt_callback)));

    // Set up the forest for adaptation using the wrapper callback.
    // Recursive adaptation is currently not supported.
    t8_forest_set_adapt (m_uncommitted_forest.value (), m_forest, detail::mesh_adapt_callback_wrapper, false);
  }

  void
  set_partition (bool set_for_coarsening = false)
  {
    if (!m_uncommitted_forest.has_value ()) {
      t8_forest_t new_forest;
      t8_forest_init (&new_forest);
      m_uncommitted_forest = new_forest;
    }
    t8_forest_set_partition (m_uncommitted_forest.value (), m_forest, set_for_coarsening);
  }

  void
  set_balance (bool no_repartition = false)
  {
    if (!m_uncommitted_forest.has_value ()) {
      t8_forest_t new_forest;
      t8_forest_init (&new_forest);
      m_uncommitted_forest = new_forest;
    }
    t8_forest_set_balance (m_uncommitted_forest.value (), m_forest, no_repartition);
  }

  void
  set_ghost (bool do_ghost = true, t8_ghost_type_t ghost_type = T8_GHOST_FACES)
  {
    if (!m_uncommitted_forest.has_value ()) {
      m_uncommitted_forest.emplace ();
      t8_forest_init (&*m_uncommitted_forest);
    }
    t8_forest_set_ghost (m_uncommitted_forest.value (), do_ghost, ghost_type);
  }

  void
  commit ()
  {
    if (!m_uncommitted_forest.has_value ()) {
      m_uncommitted_forest.emplace ();
      t8_forest_init (&*m_uncommitted_forest);
    }
    /* It can happen that the user only calls set_ghost before commit.
    This does not set the set_from member of the forest and we copy the current forest in this case. */
    if (m_uncommitted_forest.value ()->set_from == NULL) {
      t8_forest_set_copy (m_uncommitted_forest.value (), m_forest);
    }
    t8_forest_ref (m_forest);
    t8_forest_commit (m_uncommitted_forest.value ());
    // Check if we adapted and unregister the adapt context if so.
    if (detail::adapt_registry::get (m_forest) != nullptr) {
      detail::adapt_registry::unregister_context (m_forest);
      if constexpr (has_element_data_handler_competence ()) {
        t8_global_infof (
          "Please note that the element data is not interpolated automatically during adaptation. Use the "
          "function set_element_data() to provide new adapted element data.\n");
      }
    }
    t8_forest_unref (&m_forest);
    // Update underlying forest of the mesh.
    m_forest = m_uncommitted_forest.value ();
    m_uncommitted_forest.reset ();
    update_elements ();
  }

  // --- Methods to check for mesh competences. ---
  static constexpr bool
  has_element_data_handler_competence ()
  {
    return requires (SelfType& mesh) { mesh.get_element_data (); };
  }

 private:
  void
  update_elements ()
  {
    m_elements.clear ();
    m_elements.reserve (get_num_local_elements ());
    // Iterate through forest elements and fill the element vector with newly created mesh elements.
    for (t8_locidx_t itree = 0; itree < t8_forest_get_num_local_trees (m_forest); ++itree) {
      const t8_locidx_t num_elems = t8_forest_get_tree_num_leaf_elements (m_forest, itree);
      for (t8_locidx_t ielem = 0; ielem < num_elems; ++ielem) {
        m_elements.push_back (element_class (this, itree, ielem));
      }
    }
    update_ghost_elements ();
  }

  void
  update_ghost_elements ()
  {
    if (get_num_ghosts () == 0) {
      m_ghosts.clear ();
      return;
    }
    m_ghosts.clear ();
    m_ghosts.reserve (get_num_ghosts ());
    t8_locidx_t num_loc_trees = t8_forest_get_num_local_trees (m_forest);

    for (t8_locidx_t itree = 0; itree < t8_forest_get_num_ghost_trees (m_forest); ++itree) {
      const t8_locidx_t num_elems = t8_forest_ghost_tree_num_leaf_elements (m_forest, itree);
      for (t8_locidx_t ielem = 0; ielem < num_elems; ++ielem) {
        m_ghosts.push_back (element_class (this, num_loc_trees + itree, ielem, true));
      }
    }
  }

  t8_forest_t m_forest;
  std::vector<element_class> m_elements;
  std::vector<element_class> m_ghosts;
  std::optional<t8_forest_t>
    m_uncommitted_forest;
};

}  // namespace t8_mesh_handle