Util.h

Go to the documentation of this file.

#ifndef COMMA_PROJECT_UTIL_H
#define COMMA_PROJECT_UTIL_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 <algorithm>
#include <cassert>
#include <cmath>
#include <functional>
#include <iterator>
#include <map>
#include <numeric>
#include <queue>
#include <set>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
 
#include "CoMMA/Args.h"
 
namespace comma {
 
#define CoMMAUnused(var) (void)(var)
 
constexpr double deviate_thresh = 0.9396926207859084;
 
template<typename T>
inline T constexpr _sq(const T x) {
  return x * x;
}
 
template<typename T>
inline T constexpr _cb(const T x) {
  return x * x * x;
}
 
template<unsigned int p, typename T>
inline T constexpr int_power(const T x) {
  // See https://stackoverflow.com/a/1505791/12152457
  if constexpr (p == 0) return 1;
  if constexpr (p == 1) return x;
  if constexpr (p == 2) return x * x;
  if constexpr (p == 3) return x * x * x;
 
  constexpr int tmp = int_power<p / 2>(x);
  if constexpr ((p % 2) == 0) {
    return tmp * tmp;
  } else {
    return x * tmp * tmp;
  }
}
 
template<typename T>
inline bool constexpr dot_deviate(const T dot) {
  return fabs(dot) < decltype(fabs(dot)){deviate_thresh};
}
 
template<typename T>
inline T dot_product(const std::vector<T> &a, const std::vector<T> &b) {
#if 0
  const CoMMAWeightType dot = inner_product(
      prev_dir.begin(), prev_dir.end(), cur_dir.begin(),
      CoMMAWeightType{0.});
#endif
  T dot{0.};
  for (auto i = decltype(a.size()){0}; i < a.size(); ++i)
    dot += a[i] * b[i];
  return dot;
}
 
template<typename T>
inline T get_direction(
  const std::vector<T> &a, const std::vector<T> &b, std::vector<T> &dir
) {
  T norm{0.};
  for (auto i = decltype(a.size()){0}; i < a.size(); ++i) {
    const T di = b[i] - a[i];
    dir[i] = di;
    norm += di * di;
  }
  norm = sqrt(norm);
  const T ov_norm = T{1.} / norm;
  for (auto &di : dir)
    di *= ov_norm;
  return norm;
}
 
template<typename T>
inline T squared_euclidean_distance(
  const std::vector<T> &a, const std::vector<T> &b
) {
#if 0
  return
      transform_reduce(a.cbegin(), a.cend(), b.cbegin(), T{0.},
                       [](const auto sum, const auto val){return sum + val*val;},
                       minus<T>());
#endif
  T dist{0.};
  for (auto i = decltype(a.size()){0}; i < a.size(); ++i) {
    const auto diff = a[i] - b[i];
    dist += diff * diff;
  }
  return dist;
}
 
template<class PairT>
struct CustomPairLessFunctor {
  inline bool operator()(const PairT &left, const PairT &right) const {
    if (left.second < right.second) return true;
    if (left.second > right.second) return false;
    /* left.second == right.second */
    return left.first > right.first;
  }
};
 
template<typename PairT>
struct CustomPairGreaterFunctor {
  inline bool operator()(const PairT &left, const PairT &right) const {
    // I tried this: return !CustomPairLessFunctor<PairT>()(left,right);
    // but it didn't work well for equality
    if (left.second > right.second) return true;
    if (left.second < right.second) return false;
    /* left.second == right.second */
    return left.first < right.first;
  }
};
 
template<typename PairT>
struct PairSecondBasedLessFunctor {
  inline bool operator()(const PairT &left, const PairT &right) const {
    return left.second < right.second;
  }
};
 
template<typename CoMMAContainerPairType>
inline std::vector<typename CoMMAContainerPairType::value_type::first_type>
vector_of_first_elements(const CoMMAContainerPairType &cont) {
  std::vector<typename CoMMAContainerPairType::value_type::first_type> fe;
  fe.reserve(cont.size());
  for (const auto &pr : cont)
    fe.emplace_back(pr.first);
  return fe;
}
 
template<typename PairT>
class PairFindFirstBasedFunctor {
public:
  PairFindFirstBasedFunctor() :
    _target(){};
  // NOLINTNEXTLINE
  PairFindFirstBasedFunctor(const typename PairT::first_type &target) :
    _target(target){};
  ~PairFindFirstBasedFunctor() = default;
 
  inline bool operator()(const PairT &pr) const { return pr.first == _target; }
 
private:
  typename PairT::first_type _target;
};
 
template<typename KeyT, typename ValueT>
inline std::unordered_set<KeyT> d_keys_to_set(
  const std::unordered_map<KeyT, ValueT> &dict
) {
  std::unordered_set<KeyT> s_neighbours_of_seed = {};
  for (const auto &i_k_v : dict) {
    s_neighbours_of_seed.insert(i_k_v.first);
  }
  return s_neighbours_of_seed;
}
 
template<typename IndexT, typename IntT>
inline void filter_cells_by_n_edges(
  const std::vector<IndexT> &indices,
  const std::vector<IntT> &n_bnd_faces,
  const std::unordered_set<IntT> &allowed,
  std::vector<IndexT> &filtered
) {
  filtered.clear();
  auto target = decltype(allowed){};
  std::copy_if(
    allowed.begin(),
    allowed.end(),
    std::inserter(target, target.begin()),
    [](const auto n) { return n > 0; }
  );
  if (!target.empty()) {
    const auto n_cells = static_cast<IndexT>(n_bnd_faces.size());
    filtered.reserve(n_cells);
    for (auto c = decltype(n_cells){0}; c < n_cells; ++c) {
      if (allowed.find(
            static_cast<IntT>(indices[c + 1] - indices[c] + n_bnd_faces[c])
          )
          != allowed.end())
        filtered.emplace_back(c);
    }
    filtered.shrink_to_fit();
  }
}
 
template<typename IndexT, typename DistT>
void compute_neighbourhood_based_wall_distance(
  const std::vector<IndexT> &neigh_idxs,
  const std::vector<IndexT> &neighs,
  const std::vector<IndexT> &wall,
  std::vector<DistT> &dist
) {
  static_assert(
    std::is_signed<DistT>::value,
    "The distance type should be signed to allow flags (negative values)"
  );
  dist.resize(neigh_idxs.size() - 1);
  std::fill(dist.begin(), dist.end(), DistT{-1});
  std::queue<IndexT> to_visit{};
  for (const auto &cell : wall) {
    dist[cell] = DistT{0};
    to_visit.emplace(cell);
  }
  while (!to_visit.empty()) {
    // If you get wrong results, this could might be to the fact that one should
    // consider not specially the first of the queue, but the one with lowest
    // distance visited so far. So, the definition of `cur` below should be
    // changed, and with that, another container should be used instead of a
    // queue, since one might pop from whenever and not only the front (a deque
    // should work).
    const auto cur = to_visit.front();
    to_visit.pop();
    const auto cur_dist = dist[cur] + DistT{1};
    for (auto neigh = neighs.cbegin() + neigh_idxs[cur];
         neigh < neighs.cbegin() + neigh_idxs[cur + 1];
         ++neigh) {
      if (dist[*neigh] < DistT{0}) {
        dist[*neigh] = cur_dist;
        to_visit.emplace(*neigh);
      } else if (dist[*neigh] > cur_dist) {
        dist[*neigh] = cur_dist;
      }
    }
  }
}
 
template<typename IndexT, typename RealT, typename IntT>
void build_coarse_graph(
  const std::vector<IndexT> &f2c,
  const std::vector<IndexT> &f_adj_idx,
  const std::vector<IndexT> &f_adj,
  const std::vector<RealT> &f_weights,
  const std::vector<RealT> &f_volumes,
  const std::vector<std::vector<RealT>> &f_centers,
  const std::vector<RealT> &f_priority,
  const std::vector<IntT> &f_n_bnd,
  // output
  std::vector<IndexT> &c_adj_idx,
  std::vector<IndexT> &c_adj,
  std::vector<RealT> &c_weights,
  std::vector<RealT> &c_volumes,
  std::vector<std::vector<RealT>> &c_centers,
  std::vector<RealT> &c_priority,
  std::vector<IntT> &c_n_bnd
) {
  const auto n_cells = *(std::max_element(f2c.begin(), f2c.end())) + 1;
  const auto dim = static_cast<IntT>(f_centers[0].size());
  // Preparing outputs
  c_adj_idx.clear();
  c_adj_idx.emplace_back(0);
  c_adj.clear();
  c_weights.clear();
  // For these, we have the final size
  c_volumes.clear();
  c_volumes.resize(n_cells);
  std::fill(c_volumes.begin(), c_volumes.end(), RealT{0});
  c_n_bnd.clear();
  c_n_bnd.resize(n_cells);
  c_priority.clear();
  c_priority.resize(n_cells);
  std::fill(c_n_bnd.begin(), c_n_bnd.end(), IntT{0});
  c_centers.clear();
  c_centers.reserve(n_cells);
  for (auto cc = decltype(n_cells){0}; cc < n_cells; ++cc) {
    // https://stackoverflow.com/questions/18189362/how-best-to-fill-a-vector-of-vectors-avoiding-wasting-memory-and-unnecessary-al
    auto tmp = std::vector<RealT>(dim, RealT{0});
    c_centers.emplace_back(std::move(tmp));
  }  // for cc
  // Now building
  // For each coarse cell, find its fine cells
  std::vector<std::set<IndexT>> c2f(n_cells);
  for (auto fc = decltype(f2c.size()){0}; fc < f2c.size(); ++fc) {
    const auto cc = f2c[fc];
    // Adding FC to CC
    c2f[cc].emplace(fc);
    // Adding volume
    c_volumes[cc] += f_volumes[fc];
    // Updating center - Weighted average
    for (auto xyz = decltype(dim){0}; xyz < dim; ++xyz)
      c_centers[cc][xyz] += f_volumes[fc] * f_centers[fc][xyz];
  }  // for fc
  for (auto cc = decltype(n_cells){0}; cc < n_cells; ++cc) {
    const auto &fcs = c2f[cc];
    const RealT ov_n_fc(RealT{1.} / fcs.size());
    // Finishing center
    for (auto xyz = decltype(dim){0}; xyz < dim; ++xyz)
      c_centers[cc][xyz] *= ov_n_fc;
    auto n_bnd = std::numeric_limits<IntT>::min();
    auto priority = std::numeric_limits<RealT>::min();
    std::map<IndexT, RealT> cc_neighs;
    for (const auto fc : fcs) {
      if (f_n_bnd[fc] > n_bnd) n_bnd = f_n_bnd[fc];
      if (f_priority[fc] > priority) priority = f_priority[fc];
      for (auto fc_neigh = f_adj.cbegin() + f_adj_idx[fc];
           fc_neigh != f_adj.cbegin() + f_adj_idx[fc + 1];
           ++fc_neigh) {
        // Update if not one of the FC of the CC
        if (fcs.find(*fc_neigh) == fcs.end()) {
          const auto cc_neigh = f2c[*fc_neigh];
          // Add if not yet present
          cc_neighs.try_emplace(cc_neigh, 0.0);
          cc_neighs.at(cc_neigh) += f_weights[*fc_neigh];
        }
      }  // for fc_neigh
    }  // for fc
    // Max of f_bnd
    c_n_bnd[cc] = n_bnd;
    c_priority[cc] = priority;
    c_adj_idx.emplace_back(c_adj_idx.back() + cc_neighs.size());
    for (const auto &[neigh, weight] : cc_neighs) {
      c_adj.emplace_back(neigh);
      c_weights.emplace_back(weight);
    }
  }  // for cc
}
 
}  // end namespace comma
 
#endif  // COMMA_PROJECT_UTIL_H