Program Listing for File element.hxx

Return to documentation for file (mesh_handle/element.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) 2025 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 "data_handler.hxx"
#include <t8_eclass/t8_eclass.h>
#include <t8_element/t8_element.h>
#include <t8_forest/t8_forest_general.h>
#include <t8_forest/t8_forest_geometrical.h>
#include <t8_schemes/t8_scheme.hxx>
#include <t8_types/t8_vec.hxx>
#include <vector>
#include <optional>

namespace t8_mesh_handle
{
template <typename TMeshClass, template <typename> class... TCompetences>
class element: public TCompetences<element<TMeshClass, TCompetences...>>... {
 private:
  using SelfType = element<TMeshClass, TCompetences...>;
  friend TMeshClass;
  friend struct element_data_element_competence<
    SelfType>;
  element (TMeshClass* mesh, t8_locidx_t tree_id, t8_locidx_t element_id, bool is_ghost_element = false)
    : m_mesh (mesh), m_tree_id (tree_id), m_element_id (element_id), m_is_ghost_element (is_ghost_element)
  {
    // Cache the t8_element_t from the forest as it is often used.
    if (m_is_ghost_element) {
      // The local ghost tree id is per definition the local tree id - number of local (non-ghost) trees.
      m_element = t8_forest_ghost_get_leaf_element (
        m_mesh->m_forest, m_tree_id - t8_forest_get_num_local_trees (m_mesh->m_forest), m_element_id);
    }
    else {
      m_element = t8_forest_get_leaf_element_in_tree (m_mesh->m_forest, m_tree_id, m_element_id);
    }

    // Resize cache vectors for face properties to ensure clean access in case they may be only partially filled later.
    if constexpr (has_face_neighbor_cache ()) {
      const int num_faces = get_num_faces ();
      this->m_num_neighbors.resize (num_faces);
      this->m_dual_faces.resize (num_faces);
      this->m_neighbors.resize (num_faces);
    }
    if constexpr (has_face_areas_cache ()) {
      const int num_faces = get_num_faces ();
      this->m_face_areas.resize (num_faces);
    }
    if constexpr (has_face_centroids_cache ()) {
      const int num_faces = get_num_faces ();
      this->m_face_centroids.resize (num_faces);
    }
    if constexpr (has_face_normals_cache ()) {
      const int num_faces = get_num_faces ();
      this->m_face_normals.resize (num_faces);
    }
  }

 public:
  // --- Public functions to check if caches exist. ---
  static constexpr bool
  has_volume_cache ()
  {
    return requires (SelfType& element) { element.volume_cache_filled (); };
  }

  static constexpr bool
  has_diameter_cache ()
  {
    return requires (SelfType& element) { element.diameter_cache_filled (); };
  }

  static constexpr bool
  has_vertex_cache ()
  {
    return requires (SelfType& element) { element.vertex_cache_filled (); };
  }

  static constexpr bool
  has_centroid_cache ()
  {
    return requires (SelfType& element) { element.centroid_cache_filled (); };
  }

  static constexpr bool
  has_face_neighbor_cache ()
  {
    return requires (SelfType& element) { element.neighbor_cache_filled (0); };
  }
  static constexpr bool
  has_face_areas_cache ()
  {
    return requires (SelfType& element) { element.face_area_cache_filled (0); };
  }

  static constexpr bool
  has_face_centroids_cache ()
  {
    return requires (SelfType& element) { element.face_centroid_cache_filled (0); };
  }

  static constexpr bool
  has_face_normals_cache ()
  {
    return requires (SelfType& element) { element.face_normal_cache_filled (0); };
  }

  static constexpr bool
  has_element_data_handler_competence ()
  {
    return requires (SelfType& element) { element.get_element_data (); };
  }

  // --- Functionality of the element. In each function, it is checked if a cached version exists (and is used then). ---
  t8_element_level
  get_level () const
  {
    return t8_forest_get_scheme (m_mesh->m_forest)->element_get_level (get_tree_class (), m_element);
  }

  int
  get_num_faces () const
  {
    return t8_forest_get_scheme (m_mesh->m_forest)->element_get_num_faces (get_tree_class (), m_element);
  }

  int
  get_num_vertices () const
  {
    return t8_forest_get_scheme (m_mesh->m_forest)->element_get_num_corners (get_tree_class (), m_element);
  }

  t8_element_shape_t
  get_shape () const
  {
    return t8_forest_get_scheme (m_mesh->m_forest)->element_get_shape (get_tree_class (), m_element);
  }

  double
  get_volume () const
  {
    if constexpr (has_volume_cache ()) {
      if (!this->volume_cache_filled ()) {
        // Fill cache.
        this->m_volume = t8_forest_element_volume (m_mesh->m_forest, m_tree_id, m_element);
      }
      return this->m_volume.value ();
    }
    return t8_forest_element_volume (m_mesh->m_forest, m_tree_id, m_element);
  }

  double
  get_diameter () const
  {
    if constexpr (has_diameter_cache ()) {
      if (!this->diameter_cache_filled ()) {
        // Fill cache.
        this->m_diameter = t8_forest_element_diam (m_mesh->m_forest, m_tree_id, m_element);
      }
      return this->m_diameter.value ();
    }
    return t8_forest_element_diam (m_mesh->m_forest, m_tree_id, m_element);
  }

  std::vector<t8_3D_vec>
  get_vertex_coordinates () const
  {
    if constexpr (has_vertex_cache ()) {
      if (this->vertex_cache_filled ()) {
        return this->m_vertex_coordinates;
      }
    }
    // Calculate the vertex coordinates.
    const int num_corners = get_num_vertices ();
    std::vector<t8_3D_vec> vertex_coordinates (num_corners);
    for (int icorner = 0; icorner < num_corners; ++icorner) {
      t8_forest_element_coordinate (m_mesh->m_forest, m_tree_id, m_element, icorner,
                                    vertex_coordinates[icorner].data ());
    }
    // Fill the cache in the cached version.
    if constexpr (has_vertex_cache ()) {
      this->m_vertex_coordinates = std::move (vertex_coordinates);
      return this->m_vertex_coordinates;
    }
    return vertex_coordinates;
  }

  t8_3D_vec
  get_vertex_coordinates (int vertex) const
  {
    T8_ASSERT (vertex < get_num_vertices ());
    // Check if we have a cached version and if the cache has already been filled.
    if constexpr (has_vertex_cache ()) {
      if (!this->vertex_cache_filled ()) {
        get_vertex_coordinates ();
      }
      return this->m_vertex_coordinates[vertex];
    }
    t8_3D_vec coordinates;
    t8_forest_element_coordinate (m_mesh->m_forest, m_tree_id, m_element, vertex, coordinates.data ());
    return coordinates;
  }

  t8_3D_vec
  get_centroid () const
  {
    // Check if we have a cached version and if the cache has already been filled.
    if constexpr (has_centroid_cache ()) {
      if (this->centroid_cache_filled ()) {
        return this->m_centroid.value ();
      }
    }
    t8_3D_vec coordinates;
    t8_forest_element_centroid (m_mesh->m_forest, m_tree_id, m_element, coordinates.data ());
    // Fill the cache in the cached version.
    if constexpr (has_centroid_cache ()) {
      this->m_centroid = coordinates;
    }
    return coordinates;
  }

  std::vector<const SelfType*>
  get_face_neighbors (int face, std::optional<std::reference_wrapper<std::vector<int>>> dual_faces = std::nullopt) const
  {
    T8_ASSERT ((face >= 0) && (face < get_num_faces ()));
    SC_CHECK_ABORT (!m_is_ghost_element, "get_face_neighbors is not implemented for ghost elements.\n");
    if constexpr (has_face_neighbor_cache ()) {
      if (this->neighbor_cache_filled (face)) {
        if (dual_faces) {
          dual_faces->get () = this->m_dual_faces[face];
        }
        return this->m_neighbors[face];
      }
    }
    const t8_element_t** neighbors;
    int* dual_faces_internal;
    int num_neighbors;
    t8_locidx_t* neighids;
    t8_eclass_t neigh_class;
    t8_forest_leaf_face_neighbors (m_mesh->m_forest, m_tree_id, m_element, &neighbors, face, &dual_faces_internal,
                                   &num_neighbors, &neighids, &neigh_class);
    if (dual_faces) {
      dual_faces->get ().assign (dual_faces_internal, dual_faces_internal + num_neighbors);
    }
    std::vector<const SelfType*> neighbors_handle;
    for (int ineighs = 0; ineighs < num_neighbors; ineighs++) {
      neighbors_handle.push_back (&((*m_mesh)[neighids[ineighs]]));
    }
    if constexpr (has_face_neighbor_cache ()) {
      // Also store num_neighbors in cache to indicate that the cache is filled if a face does not have any neighbor.
      this->m_num_neighbors[face] = num_neighbors;
      this->m_dual_faces[face].assign (dual_faces_internal, dual_faces_internal + num_neighbors);
      this->m_neighbors[face] = std::move (neighbors_handle);
    }
    if (num_neighbors > 0) {
      // Free allocated memory.
      T8_FREE (neighbors);
      T8_FREE (dual_faces_internal);
      T8_FREE (neighids);
    }
    if constexpr (has_face_neighbor_cache ()) {
      return this->m_neighbors[face];
    }
    return neighbors_handle;
  }

  void
  fill_face_neighbor_cache () const
    requires (has_face_neighbor_cache ())
  {
    for (int iface = 0; iface < get_num_faces (); iface++) {
      get_face_neighbors (iface);
    }
  }

  // --- Getter for face properties. ---
  double
  get_face_area (int face) const
  {
    T8_ASSERT ((face >= 0) && (face < get_num_faces ()));
    if constexpr (has_face_areas_cache ()) {
      if (!this->face_area_cache_filled (face)) {
        this->m_face_areas[face] = t8_forest_element_face_area (m_mesh->m_forest, m_tree_id, m_element, face);
      }
      return this->m_face_areas[face].value ();
    }
    return t8_forest_element_face_area (m_mesh->m_forest, m_tree_id, m_element, face);
  }

  t8_3D_vec
  get_face_centroid (int face) const
  {
    T8_ASSERT ((face >= 0) && (face < get_num_faces ()));
    if constexpr (has_face_centroids_cache ()) {
      if (this->face_centroid_cache_filled (face)) {
        return this->m_face_centroids[face].value ();
      }
    }
    t8_3D_vec coordinates;
    t8_forest_element_face_centroid (m_mesh->m_forest, m_tree_id, m_element, face, coordinates.data ());
    if constexpr (has_face_centroids_cache ()) {
      this->m_face_centroids[face] = coordinates;
    }
    return coordinates;
  }

  t8_3D_vec
  get_face_normal (int face) const
  {
    T8_ASSERT ((face >= 0) && (face < get_num_faces ()));
    if constexpr (has_face_normals_cache ()) {
      if (this->face_normal_cache_filled (face)) {
        return this->m_face_normals[face].value ();
      }
    }
    t8_3D_vec normal;
    t8_forest_element_face_normal (m_mesh->m_forest, m_tree_id, m_element, face, normal.data ());
    if constexpr (has_face_normals_cache ()) {
      this->m_face_normals[face] = normal;
    }
    return normal;
  }

  t8_element_shape_t
  get_face_shape (int face) const
  {
    T8_ASSERT ((face >= 0) && (face < get_num_faces ()));
    return t8_forest_get_scheme (m_mesh->m_forest)->element_get_face_shape (get_tree_class (), m_element, face);
  }

  // --- Print for the element for debugging purpose. ---
#if T8_ENABLE_DEBUG
  void
  print_element_debug () const
  {
    t8_forest_get_scheme (m_mesh->m_forest)->element_debug_print (get_tree_class (), m_element);
  }
#endif

  // --- Function to access mesh specific id. ---
  t8_locidx_t
  get_element_handle_id () const
  {
    if (m_is_ghost_element) {
      return m_mesh->get_num_local_elements ()
             + t8_forest_ghost_get_tree_element_offset (m_mesh->m_forest,
                                                        m_tree_id - t8_forest_get_num_local_trees (m_mesh->m_forest))
             + m_element_id;
    }
    return t8_forest_get_tree_element_offset (m_mesh->m_forest, m_tree_id) + m_element_id;
  }

  //--- Getter for the member variables. ---
  t8_locidx_t
  get_local_tree_id () const
  {
    return m_tree_id;
  }

  t8_locidx_t
  get_local_element_id () const
  {
    return m_element_id;
  }

  const TMeshClass*
  get_mesh () const
  {
    return m_mesh;
  }

  bool
  is_ghost_element () const
  {
    return m_is_ghost_element;
  }

 private:
  // --- Private member variables. ---
  TMeshClass* m_mesh;
  const t8_locidx_t m_tree_id;
  const t8_locidx_t m_element_id;
  const bool m_is_ghost_element;
  const t8_element_t* m_element;
  // --- Private helper function. ---
  t8_eclass_t
  get_tree_class () const
  {
    return t8_forest_get_tree_class (m_mesh->m_forest, m_tree_id);
  }
};

}  // namespace t8_mesh_handle