ARComputer.h

Go to the documentation of this file.

#ifndef COMMA_PROJECT_ARCOMPUTER_H
#define COMMA_PROJECT_ARCOMPUTER_H
 
/*
 * CoMMA
 *
 * Copyright © 2024 ONERA
 *
 * Authors: Nicolas Lantos, Alberto Remigi, and Riccardo Milani
 * Contributors: Karim Anemiche
 *
 * This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at https://mozilla.org/MPL/2.0/.
 */
 
#include <cmath>
#include <memory>
#include <numeric>
#include <unordered_map>
#include <unordered_set>
#include <vector>
 
#include "CoMMA/Dual_Graph.h"
#include "CoMMA/Util.h"
 
namespace comma {
 
template<typename IndexT, typename RealT, typename IntT>
class CellFeatures {
public:
  RealT _measure;
 
  RealT _external_weights;
 
  RealT _internal_weights;
 
  IntT _n_internal_faces;
 
  RealT _sq_min_edge;
 
  RealT _sq_diam;
 
  std::vector<RealT> _sum_centers;
 
  std::unordered_map<IndexT, IntT> _external_facets;
 
  explicit CellFeatures(
    const RealT measure = std::numeric_limits<RealT>::min(),
    const RealT external_weights = 0.0,
    const RealT internal_weights = 0.0,
    const IntT n_internal_faces = 0,
    const RealT min_edge = std::numeric_limits<RealT>::max(),
    const RealT diam = std::numeric_limits<RealT>::min(),
    const std::vector<RealT> sum_centers = std::vector<RealT>{},
    const std::unordered_map<IndexT, IntT> external_facets =
      std::unordered_map<IndexT, IntT>{}
  ) :
    _measure(measure),
    _external_weights(external_weights),
    _internal_weights(internal_weights),
    _n_internal_faces(n_internal_faces),
    _sq_min_edge(min_edge * min_edge),
    _sq_diam(diam * diam),
    _sum_centers(sum_centers),
    _external_facets(external_facets){};
 
  CellFeatures(
    const IndexT index, std::shared_ptr<Dual_Graph<IndexT, RealT, IntT>> graph
  ) :
    _measure(graph->_volumes[index]),
    _external_weights(graph->estimated_total_weight(index)),
    _internal_weights(0.0),
    _n_internal_faces(0),
    _sq_min_edge(std::numeric_limits<RealT>::max()),
    _sq_diam(std::numeric_limits<RealT>::min()),
    _external_facets({{index, graph->get_total_n_faces(index)}}) {
    _sum_centers = graph->_centers[index];
    for (auto &xyz : _sum_centers) {
      xyz *= graph->_volumes[index];
    }
  }
 
  CellFeatures(const CellFeatures<IndexT, RealT, IntT> &other) = default;
 
  CellFeatures(CellFeatures<IndexT, RealT, IntT> &&other) noexcept = default;
 
  CellFeatures &operator=(CellFeatures<IndexT, RealT, IntT> &&other
  ) noexcept = default;
 
  CellFeatures &operator=(const CellFeatures<IndexT, RealT, IntT> &other
  ) = default;
 
  template<unsigned int dim>
  inline RealT constexpr get_radius() const {
    if constexpr (dim < 2) {
      return _measure;
    } else if constexpr (dim == 2) {
      return sqrt(_measure);
    } else if constexpr (dim == 3) {
      return cbrt(_measure);
    } else {
      return pow(_measure, 1.0 / dim);
    }
  }
};
 
template<typename IndexT, typename RealT, typename IntT>
class ARComputer {
public:
  explicit ARComputer(std::shared_ptr<Dual_Graph<IndexT, RealT, IntT>> graph) :
    _graph(graph){};
 
  virtual ~ARComputer() = default;
 
  virtual void compute_and_update_features(
    const IndexT i_fc,
    const CellFeatures<IndexT, RealT, IntT> &cc_feats,
    const std::unordered_set<IndexT> &fc_of_cc,
    // out
    IntT &shared_faces,
    RealT &aspect_ratio,
    CellFeatures<IndexT, RealT, IntT> &new_feats
  ) const = 0;
 
protected:
  std::shared_ptr<Dual_Graph<IndexT, RealT, IntT>> _graph;
 
  inline void compute_shared_faces(
    const IndexT i_fc,
    const std::unordered_set<IndexT> &fc_of_cc,
    RealT &shared_weights,
    std::unordered_set<IndexT> &shared_faces
  ) const {
    shared_weights = 0.;
    auto n_it = this->_graph->neighbours_cbegin(i_fc);
    auto w_it = this->_graph->weights_cbegin(i_fc);
    for (; n_it != this->_graph->neighbours_cend(i_fc); ++n_it, ++w_it) {
      if (*n_it != i_fc && (fc_of_cc.count(*n_it) != 0)) {
        shared_faces.insert(*n_it);
        shared_weights += *w_it;
      }
    }
  }
 
  inline void update_approximated_geometric_features(
    const IndexT i_fc,
    const CellFeatures<IndexT, RealT, IntT> &cc_feats,
    const std::unordered_set<IndexT> &fc_of_cc,
    CellFeatures<IndexT, RealT, IntT> &new_feats
  ) const {
    const std::vector<RealT> &cen_fc = this->_graph->_centers[i_fc];
    RealT max_diam = cc_feats._sq_diam, min_edge = cc_feats._sq_min_edge;
    for (const auto i_fc_cc : fc_of_cc) {
      const auto dist = squared_euclidean_distance<RealT>(
        cen_fc, this->_graph->_centers[i_fc_cc]
      );
      if (dist > max_diam) max_diam = dist;
      if (dist < min_edge) min_edge = dist;
    }  // for i_fc_cc
    new_feats._sq_diam = max_diam;
    new_feats._sq_min_edge = min_edge;
  }
 
  inline std::vector<RealT> compute_and_update_barycenter(
    const IndexT i_fc,
    const CellFeatures<IndexT, RealT, IntT> &cc_feats,
    const std::unordered_set<IndexT> &fc_of_cc,
    CellFeatures<IndexT, RealT, IntT> &new_feats
  ) const {
    new_feats._sum_centers = cc_feats._sum_centers;
    const auto &cur_c = this->_graph->_centers[i_fc];
    std::vector<RealT> bary(cur_c.size());
    const RealT ov_N = 1. / (fc_of_cc.size() + 1);
    for (auto i = decltype(cur_c.size()){0}; i < cur_c.size(); ++i) {
      new_feats._sum_centers[i] += this->_graph->_volumes[i_fc] * cur_c[i];
      bary[i] = new_feats._sum_centers[i] * ov_N;
    }
    return bary;
  }
 
  template<bool compute_weights = false, bool update_facets = false>
  inline void update_basic_features(
    const IndexT i_fc,
    const CellFeatures<IndexT, RealT, IntT> &cc_feats,
    const std::unordered_set<IndexT> &fc_of_cc,
    IntT &n_shared_faces,
    CellFeatures<IndexT, RealT, IntT> &new_feats
  ) const {
    new_feats._measure = cc_feats._measure + this->_graph->_volumes[i_fc];
    RealT shared_weights{};
    std::unordered_set<IndexT> shared_faces{};
    this->compute_shared_faces(i_fc, fc_of_cc, shared_weights, shared_faces);
    n_shared_faces = static_cast<IntT>(shared_faces.size());
    new_feats._n_internal_faces = cc_feats._n_internal_faces + n_shared_faces;
    if constexpr (compute_weights) {
      new_feats._external_weights = cc_feats._external_weights
        + this->_graph->estimated_total_weight(i_fc) - 2 * shared_weights;
      new_feats._internal_weights = cc_feats._internal_weights + shared_weights;
    }
    if constexpr (update_facets) {
      new_feats._external_facets = cc_feats._external_facets;  // OK copy
      for (const auto &face : shared_faces) {
        new_feats._external_facets.at(face)--;
      }
      new_feats._external_facets[i_fc] =
        this->_graph->get_total_n_faces(i_fc) - n_shared_faces;
    }
  }
 
  template<unsigned int dim>
  inline RealT constexpr x_over_radius(
    const RealT sqx, const CellFeatures<IndexT, RealT, IntT> &feats
  ) const {
    if constexpr (dim < 2) {
      return sqrt(sqx) / feats.template get_radius<dim>();
    } else if constexpr (dim == 2) {
      return sqx / feats._measure;
    } else if constexpr (dim == 3) {
      return _cb(sqx) / _sq(feats._measure);
    } else {
      return int_power<dim>(sqx) / _sq(feats._measure);
    }
  }
};
 
template<typename IndexT, typename RealT, typename IntT, unsigned int dim>
class ARDiamOverRadius : public ARComputer<IndexT, RealT, IntT> {
public:
  explicit ARDiamOverRadius(
    std::shared_ptr<Dual_Graph<IndexT, RealT, IntT>> graph
  ) :
    ARComputer<IndexT, RealT, IntT>(graph){};
 
  ~ARDiamOverRadius() override = default;
 
  void compute_and_update_features(
    const IndexT i_fc,
    const CellFeatures<IndexT, RealT, IntT> &cc_feats,
    const std::unordered_set<IndexT> &fc_of_cc,
    // out
    IntT &shared_faces,
    RealT &aspect_ratio,
    CellFeatures<IndexT, RealT, IntT> &new_feats
  ) const override {
    // Update
    this->template update_basic_features<>(
      i_fc, cc_feats, fc_of_cc, shared_faces, new_feats
    );
    this->update_approximated_geometric_features(
      i_fc, cc_feats, fc_of_cc, new_feats
    );
    // Compute AR
    aspect_ratio =
      this->template x_over_radius<dim>(new_feats._sq_diam, new_feats);
  }
};
 
template<typename IndexT, typename RealT, typename IntT>
class AROverMeasure : public ARComputer<IndexT, RealT, IntT> {
public:
  explicit AROverMeasure(std::shared_ptr<Dual_Graph<IndexT, RealT, IntT>> graph
  ) :
    ARComputer<IndexT, RealT, IntT>(graph){};
 
  ~AROverMeasure() override = default;
 
  void compute_and_update_features(
    const IndexT i_fc,
    const CellFeatures<IndexT, RealT, IntT> &cc_feats,
    const std::unordered_set<IndexT> &fc_of_cc,
    // out
    IntT &shared_faces,
    RealT &aspect_ratio,
    CellFeatures<IndexT, RealT, IntT> &new_feats
  ) const override {
    // Update
    this->template update_basic_features<>(
      i_fc, cc_feats, fc_of_cc, shared_faces, new_feats
    );
    // Compute AR
    aspect_ratio = 1. / new_feats._measure;
  }
};
 
template<typename IndexT, typename RealT, typename IntT>
class ARDiamOverMinEdge : public ARComputer<IndexT, RealT, IntT> {
public:
  explicit ARDiamOverMinEdge(
    std::shared_ptr<Dual_Graph<IndexT, RealT, IntT>> graph
  ) :
    ARComputer<IndexT, RealT, IntT>(graph){};
 
  ~ARDiamOverMinEdge() override = default;
 
  void compute_and_update_features(
    const IndexT i_fc,
    const CellFeatures<IndexT, RealT, IntT> &cc_feats,
    const std::unordered_set<IndexT> &fc_of_cc,
    // out
    IntT &shared_faces,
    RealT &aspect_ratio,
    CellFeatures<IndexT, RealT, IntT> &new_feats
  ) const override {
    // Update
    this->template update_basic_features<>(
      i_fc, cc_feats, fc_of_cc, shared_faces, new_feats
    );
    this->update_approximated_geometric_features(
      i_fc, cc_feats, fc_of_cc, new_feats
    );
    // Compute AR
    // If only 2 cells (hence the reference coarse cell is only one cell), max
    // and min are the same, hence use only max
    aspect_ratio = fc_of_cc.size() == 1
      ? new_feats._sq_diam
      : new_feats._sq_diam / new_feats._sq_min_edge;
  }
};
 
template<typename IndexT, typename RealT, typename IntT, unsigned int dim>
class ARExternalWeightOverRadius : public ARComputer<IndexT, RealT, IntT> {
public:
  explicit ARExternalWeightOverRadius(
    std::shared_ptr<Dual_Graph<IndexT, RealT, IntT>> graph
  ) :
    ARComputer<IndexT, RealT, IntT>(graph){};
 
  ~ARExternalWeightOverRadius() override = default;
 
  void compute_and_update_features(
    const IndexT i_fc,
    const CellFeatures<IndexT, RealT, IntT> &cc_feats,
    const std::unordered_set<IndexT> &fc_of_cc,
    // out
    IntT &shared_faces,
    RealT &aspect_ratio,
    CellFeatures<IndexT, RealT, IntT> &new_feats
  ) const override {
    // Update
    this->template update_basic_features<true>(
      i_fc, cc_feats, fc_of_cc, shared_faces, new_feats
    );
    // Compute AR
    aspect_ratio = this->compute_AR(new_feats);
  }
 
protected:
  inline RealT constexpr compute_AR(
    const CellFeatures<IndexT, RealT, IntT> &feats
  ) const {
    if constexpr (dim < 2) {
      return feats._external_weights / feats.template get_radius<dim>();
    } else if constexpr (dim == 2) {
      // return feats._external_weights / feats.template get_radius<dim>();
      return _sq(feats._external_weights) / feats._measure;
    } else if constexpr (dim == 3) {
      // return sqrt(feats._external_weights) / feats.template
      // get_radius<dim>();
      return _cb(feats._external_weights) / _sq(feats._measure);
    } else {
      // return pow(feats._external_weights, 1.0 / (dim - 1))
      // / feats.template get_radius<dim>();
      return int_power<dim>(feats._external_weights)
        / int_power<dim - 1>(feats._measure);
    }
  }
};
 
template<typename IndexT, typename RealT, typename IntT>
class ARDiameter : public ARComputer<IndexT, RealT, IntT> {
public:
  explicit ARDiameter(std::shared_ptr<Dual_Graph<IndexT, RealT, IntT>> graph) :
    ARComputer<IndexT, RealT, IntT>(graph){};
 
  ~ARDiameter() override = default;
 
  void compute_and_update_features(
    const IndexT i_fc,
    const CellFeatures<IndexT, RealT, IntT> &cc_feats,
    const std::unordered_set<IndexT> &fc_of_cc,
    // out
    IntT &shared_faces,
    RealT &aspect_ratio,
    CellFeatures<IndexT, RealT, IntT> &new_feats
  ) const override {
    // Update
    this->template update_basic_features<>(
      i_fc, cc_feats, fc_of_cc, shared_faces, new_feats
    );
    this->update_approximated_geometric_features(
      i_fc, cc_feats, fc_of_cc, new_feats
    );
    // Compute AR
    aspect_ratio = new_feats._sq_diam;
  }
};
 
template<typename IndexT, typename RealT, typename IntT>
class ARExternalWeights : public ARComputer<IndexT, RealT, IntT> {
public:
  explicit ARExternalWeights(
    std::shared_ptr<Dual_Graph<IndexT, RealT, IntT>> graph
  ) :
    ARComputer<IndexT, RealT, IntT>(graph){};
 
  ~ARExternalWeights() override = default;
 
  void compute_and_update_features(
    const IndexT i_fc,
    const CellFeatures<IndexT, RealT, IntT> &cc_feats,
    const std::unordered_set<IndexT> &fc_of_cc,
    // out
    IntT &shared_faces,
    RealT &aspect_ratio,
    CellFeatures<IndexT, RealT, IntT> &new_feats
  ) const override {
    // Update
    this->template update_basic_features<true>(
      i_fc, cc_feats, fc_of_cc, shared_faces, new_feats
    );
    // Compute AR
    aspect_ratio = new_feats._external_weights;
  }
};
 
template<typename IndexT, typename RealT, typename IntT>
class AROverInternalWeights : public ARComputer<IndexT, RealT, IntT> {
public:
  explicit AROverInternalWeights(
    std::shared_ptr<Dual_Graph<IndexT, RealT, IntT>> graph
  ) :
    ARComputer<IndexT, RealT, IntT>(graph){};
 
  ~AROverInternalWeights() override = default;
 
  void compute_and_update_features(
    const IndexT i_fc,
    const CellFeatures<IndexT, RealT, IntT> &cc_feats,
    const std::unordered_set<IndexT> &fc_of_cc,
    // out
    IntT &shared_faces,
    RealT &aspect_ratio,
    CellFeatures<IndexT, RealT, IntT> &new_feats
  ) const override {
    // Update
    this->template update_basic_features<true>(
      i_fc, cc_feats, fc_of_cc, shared_faces, new_feats
    );
    // Compute AR
    aspect_ratio = 1. / new_feats._internal_weights;
  }
};
 
template<typename IndexT, typename RealT, typename IntT, unsigned int dim>
class ARMaxBaryDistanceOverRadius : public ARComputer<IndexT, RealT, IntT> {
public:
  explicit ARMaxBaryDistanceOverRadius(
    std::shared_ptr<Dual_Graph<IndexT, RealT, IntT>> graph
  ) :
    ARComputer<IndexT, RealT, IntT>(graph){};
 
  ~ARMaxBaryDistanceOverRadius() override = default;
 
  void compute_and_update_features(
    const IndexT i_fc,
    const CellFeatures<IndexT, RealT, IntT> &cc_feats,
    const std::unordered_set<IndexT> &fc_of_cc,
    // out
    IntT &shared_faces,
    RealT &aspect_ratio,
    CellFeatures<IndexT, RealT, IntT> &new_feats
  ) const override {
    // Update
    this->template update_basic_features<>(
      i_fc, cc_feats, fc_of_cc, shared_faces, new_feats
    );
    const auto bary =
      this->compute_and_update_barycenter(i_fc, cc_feats, fc_of_cc, new_feats);
    // Compute AR
    const auto dist_fc =
      squared_euclidean_distance<RealT>(bary, this->_graph->_centers[i_fc]);
    RealT max_dist = dist_fc;
    for (const auto i_fc_cc : fc_of_cc) {
      const auto dist = squared_euclidean_distance<RealT>(
        bary, this->_graph->_centers[i_fc_cc]
      );
      if (dist > max_dist) max_dist = dist;
    }  // for i_fc_cc
    aspect_ratio = this->template x_over_radius<dim>(max_dist, new_feats);
  }
};
 
template<typename IndexT, typename RealT, typename IntT>
class ARMaxOverMinBaryDistance : public ARComputer<IndexT, RealT, IntT> {
public:
  explicit ARMaxOverMinBaryDistance(
    std::shared_ptr<Dual_Graph<IndexT, RealT, IntT>> graph,
    const RealT tolerance
  ) :
    ARComputer<IndexT, RealT, IntT>(graph), _tolerance(tolerance){};
 
  ~ARMaxOverMinBaryDistance() override = default;
 
  void compute_and_update_features(
    const IndexT i_fc,
    const CellFeatures<IndexT, RealT, IntT> &cc_feats,
    const std::unordered_set<IndexT> &fc_of_cc,
    // out
    IntT &shared_faces,
    RealT &aspect_ratio,
    CellFeatures<IndexT, RealT, IntT> &new_feats
  ) const override {
    // Update
    this->template update_basic_features<false, true>(
      i_fc, cc_feats, fc_of_cc, shared_faces, new_feats
    );
    const auto bary =
      this->compute_and_update_barycenter(i_fc, cc_feats, fc_of_cc, new_feats);
    // Compute AR
    RealT min_dist = std::numeric_limits<RealT>::max(),
          max_dist = std::numeric_limits<RealT>::min();
    // Skipping if totally internal cell
    if (new_feats._external_facets[i_fc] > 0) {
      // Since we are comparing to a threshold, keep sqrt here
      const auto dist = sqrt(
        squared_euclidean_distance<RealT>(bary, this->_graph->_centers[i_fc])
      );
      // Skipping if barycenter is the center of the cell
      if (dist > this->_tolerance) {
        min_dist = max_dist = dist;
      }
    }
    for (const auto i_fc_cc : fc_of_cc) {
      if (new_feats._external_facets[i_fc_cc] > 0) {
        const auto dist = sqrt(squared_euclidean_distance<RealT>(
          bary, this->_graph->_centers[i_fc_cc]
        ));
        if (dist > this->_tolerance) {
          if (dist > max_dist) max_dist = dist;
          if (dist < min_dist) min_dist = dist;
        }
      }
    }  // for i_fc_cc
    aspect_ratio = max_dist / min_dist;
  }
 
protected:
  RealT _tolerance;
};
 
}  // end namespace comma
 
#endif  // COMMA_PROJECT_ARCOMPUTER_H